Skip to content

Commit 4b5d0d8

Browse files
authored
Merge pull request #484 from int-brain-lab/atlas_xyz2i
Handling index out of volume
2 parents ce9f04b + 585d9da commit 4b5d0d8

File tree

7 files changed

+92
-29
lines changed

7 files changed

+92
-29
lines changed

brainbox/io/one.py

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -334,7 +334,6 @@ def _load_channel_locations_traj(eid, probe=None, one=None, revision=None, align
334334
# get the channels from histology tracing
335335
xyz = xyz[np.argsort(xyz[:, 2]), :]
336336
chans = histology.interpolate_along_track(xyz, (depths + TIP_SIZE_UM) / 1e6)
337-
338337
channels[probe] = _channels_traj2bunch(chans, brain_atlas)
339338
source = 'traced'
340339
channels[probe]['axial_um'] = chn_coords[:, 1]
@@ -894,6 +893,7 @@ class SpikeSortingLoader:
894893
collection: str = ''
895894
histology: str = '' # 'alf', 'resolved', 'aligned' or 'traced'
896895
spike_sorting_path: Path = None
896+
_sync: dict = None
897897

898898
def __post_init__(self):
899899
# pid gets precedence
@@ -1039,3 +1039,20 @@ def url(self):
10391039
"""Gets flatiron URL for the session"""
10401040
webclient = getattr(self.one, '_web_client', None)
10411041
return webclient.rel_path2url(get_alf_path(self.session_path)) if webclient else None
1042+
1043+
def samples2times(self, values, direction='forward'):
1044+
"""
1045+
:param values: numpy array of times in seconds or samples to resync
1046+
:param direction: 'forward' (samples probe time to seconds main time) or 'reverse'
1047+
(seconds main time to samples probe time)
1048+
:return:
1049+
"""
1050+
if self._sync is None:
1051+
timestamps = self.one.load_dataset(
1052+
self.eid, dataset='_spikeglx_*.timestamps.npy', collection=f'raw_ephys_data/{self.pname}')
1053+
self._sync = {
1054+
'timestamps': timestamps,
1055+
'forward': interp1d(timestamps[:, 0], timestamps[:, 1], fill_value='extrapolate'),
1056+
'reverse': interp1d(timestamps[:, 1], timestamps[:, 0], fill_value='extrapolate'),
1057+
}
1058+
return self._sync[direction](values)

ibllib/atlas/atlas.py

Lines changed: 59 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -110,22 +110,55 @@ def _round(i, round=True):
110110
else:
111111
return i
112112

113-
def x2i(self, x, round=True):
114-
return self._round((x - self.x0) / self.dx, round=round)
115-
116-
def y2i(self, y, round=True):
117-
return self._round((y - self.y0) / self.dy, round=round)
118-
119-
def z2i(self, z, round=True):
120-
return self._round((z - self.z0) / self.dz, round=round)
113+
def x2i(self, x, round=True, mode='raise'):
114+
i = np.asarray(self._round((x - self.x0) / self.dx, round=round))
115+
if np.any(i < 0) or np.any(i >= self.nx):
116+
if mode == 'clip':
117+
i[i < 0] = 0
118+
i[i >= self.nx] = self.nx - 1
119+
elif mode == 'raise':
120+
raise ValueError("At least one x value lies outside of the atlas volume.")
121+
elif mode == 'wrap':
122+
pass
123+
return i
124+
125+
def y2i(self, y, round=True, mode='raise'):
126+
i = np.asarray(self._round((y - self.y0) / self.dy, round=round))
127+
if np.any(i < 0) or np.any(i >= self.ny):
128+
if mode == 'clip':
129+
i[i < 0] = 0
130+
i[i >= self.ny] = self.ny - 1
131+
elif mode == 'raise':
132+
raise ValueError("At least one y value lies outside of the atlas volume.")
133+
elif mode == 'wrap':
134+
pass
135+
return i
136+
137+
def z2i(self, z, round=True, mode='raise'):
138+
i = np.asarray(self._round((z - self.z0) / self.dz, round=round))
139+
if np.any(i < 0) or np.any(i >= self.nz):
140+
if mode == 'clip':
141+
i[i < 0] = 0
142+
i[i >= self.nz] = self.nz - 1
143+
elif mode == 'raise':
144+
raise ValueError("At least one z value lies outside of the atlas volume.")
145+
elif mode == 'wrap':
146+
pass
147+
return i
121148

122-
def xyz2i(self, xyz, round=True):
149+
def xyz2i(self, xyz, round=True, mode='raise'):
150+
"""
151+
:param mode: {‘raise’, 'clip', 'wrap'} determines what to do when determined index lies outside the atlas volume
152+
'raise' will raise a ValueError
153+
'clip' will replace the index with the closest index inside the volume
154+
'wrap' will wrap around to the other side of the volume. This is only here for legacy reasons
155+
"""
123156
xyz = np.array(xyz)
124157
dt = int if round else float
125158
out = np.zeros_like(xyz, dtype=dt)
126-
out[..., 0] = self.x2i(xyz[..., 0], round=round)
127-
out[..., 1] = self.y2i(xyz[..., 1], round=round)
128-
out[..., 2] = self.z2i(xyz[..., 2], round=round)
159+
out[..., 0] = self.x2i(xyz[..., 0], round=round, mode=mode)
160+
out[..., 1] = self.y2i(xyz[..., 1], round=round, mode=mode)
161+
out[..., 2] = self.z2i(xyz[..., 2], round=round, mode=mode)
129162
return out
130163

131164
"""Methods indices to distance"""
@@ -227,7 +260,10 @@ def _get_cache_dir():
227260
def compute_surface(self):
228261
"""
229262
Get the volume top, bottom, left and right surfaces, and from these the outer surface of
230-
the image volume. This is needed to compute probe insertions intersections
263+
the image volume. This is needed to compute probe insertions intersections.
264+
265+
NOTE: In places where the top or bottom surface touch the top or bottom of the atlas volume, the surface
266+
will be set to np.nan. If you encounter issues working with these surfaces check if this might be the cause.
231267
"""
232268
if self.surface is None: # only compute if it hasn't already been computed
233269
axz = self.xyz2dims[2] # this is the dv axis
@@ -439,7 +475,12 @@ def slice(self, coordinate, axis, volume='image', mode='raise', region_values=No
439475
:param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
440476
:return: 2d array or 3d RGB numpy int8 array
441477
"""
442-
index = self.bc.xyz2i(np.array([coordinate] * 3))[axis]
478+
if axis == 0:
479+
index = self.bc.x2i(np.array(coordinate), mode=mode)
480+
elif axis == 1:
481+
index = self.bc.y2i(np.array(coordinate), mode=mode)
482+
elif axis == 2:
483+
index = self.bc.z2i(np.array(coordinate), mode=mode)
443484

444485
# np.take is 50 thousand times slower than straight slicing !
445486
def _take(vol, ind, axis):
@@ -765,7 +806,10 @@ def from_dict(d, brain_atlas=None):
765806
if brain_atlas:
766807
iy = brain_atlas.bc.y2i(d['y'] / 1e6)
767808
ix = brain_atlas.bc.x2i(d['x'] / 1e6)
768-
z = brain_atlas.top[iy, ix]
809+
# Only use the brain surface value as z if it isn't NaN (this happens when the surface touches the edges
810+
# of the atlas volume
811+
if not np.isnan(brain_atlas.top[iy, ix]):
812+
z = brain_atlas.top[iy, ix]
769813
return Insertion(x=d['x'] / 1e6, y=d['y'] / 1e6, z=z,
770814
phi=d['phi'], theta=d['theta'], depth=d['depth'] / 1e6,
771815
beta=d.get('beta', 0), label=d.get('label', ''))

ibllib/ephys/neuropixel.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,3 +8,4 @@
88
', change your imports to neuropixel !', DeprecationWarning)
99

1010
from neuropixel import * # noqa
11+
from neuropixel import SITES_COORDINATES # noqa

ibllib/pipes/local_server.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -147,7 +147,7 @@ def task_queue(mode='all', lab=None, one=None):
147147
if one is None:
148148
one = ONE(cache_rest=None)
149149
if lab is None:
150-
_logger.info("Trying to infer lab from globus installation")
150+
_logger.debug("Trying to infer lab from globus installation")
151151
lab = _get_lab(one)
152152
if lab is None:
153153
_logger.error("No lab provided or found")

ibllib/pipes/tasks.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -526,7 +526,8 @@ def run_alyx_task(tdict=None, session_path=None, one=None, job_deck=None,
526526
_logger.warning(f"{tdict['name']} has unmet dependencies")
527527
# if parents are waiting or failed, set the current task status to Held
528528
# once the parents ran, the descendent tasks will be set from Held to Waiting (see below)
529-
if any(map(lambda s: s in ['Errored', 'Held', 'Empty', 'Waiting'], parent_statuses)):
529+
if any(map(lambda s: s in ['Errored', 'Held', 'Empty', 'Waiting', 'Started', 'Abandoned'],
530+
parent_statuses)):
530531
tdict = one.alyx.rest('tasks', 'partial_update', id=tdict['id'],
531532
data={'status': 'Held'})
532533
return tdict, registered_dsets

ibllib/plots/figures.py

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -478,16 +478,16 @@ def gain2level(gain):
478478
if plot_backend == 'matplotlib':
479479
_, axs = plt.subplots(1, 2, gridspec_kw={'width_ratios': [.95, .05]}, figsize=(16, 9))
480480
eqcs.append(Density(butt, fs=fs, taxis=1, ax=axs[0], title='highpass', vmin=eqc_levels[0], vmax=eqc_levels[1],
481-
cmap='Greys'))
481+
cmap='Greys-r'))
482482

483483
if destripe:
484484
dest = voltage.destripe(raw, fs=fs, channel_labels=channel_labels)
485485
_, axs = plt.subplots(1, 2, gridspec_kw={'width_ratios': [.95, .05]}, figsize=(16, 9))
486486
eqcs.append(Density(dest, fs=fs, taxis=1, ax=axs[0], title='destripe', vmin=eqc_levels[0], vmax=eqc_levels[1],
487-
cmap='Greys'))
487+
cmap='Greys-r'))
488488
_, axs = plt.subplots(1, 2, gridspec_kw={'width_ratios': [.95, .05]}, figsize=(16, 9))
489489
eqcs.append(Density((butt - dest), fs=fs, taxis=1, ax=axs[0], title='difference', vmin=eqc_levels[0],
490-
vmax=eqc_levels[1], cmap='Greys'))
490+
vmax=eqc_levels[1], cmap='Greys-r'))
491491

492492
for eqc in eqcs:
493493
y, x = np.meshgrid(ioutside, np.linspace(0, rl * 1e3, 500))
@@ -618,7 +618,7 @@ def raw_destripe(raw, fs, t0, i_plt, n_plt,
618618
Tplot = Xs.shape[1] / fs
619619

620620
# PLOT RAW DATA
621-
d = Density(-Xs, fs=fs, taxis=1, ax=axs[i_plt], vmin=MIN_X, vmax=MAX_X, cmap='Greys') # noqa
621+
d = Density(-Xs, fs=fs, taxis=1, ax=axs[i_plt], vmin=MIN_X, vmax=MAX_X, cmap='Greys-r') # noqa
622622
axs[i_plt].set_ylabel('')
623623
axs[i_plt].set_xlim((0, Tplot * 1e3))
624624
axs[i_plt].set_ylim((0, nc))

ibllib/tests/test_atlas.py

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -312,11 +312,11 @@ def test_sagittal_slice(self):
312312
ax.clear()
313313

314314
def test_horizontal_slice(self):
315-
ax = self.ba.plot_hslice(dv_coordinate=0.002)
315+
ax = self.ba.plot_hslice(dv_coordinate=-0.002)
316316
im = ax.get_images()[0]
317317
assert im.get_array().shape == (self.ba.bc.ny, self.ba.bc.nx)
318318
ax.clear()
319-
ax = self.ba.plot_hslice(dv_coordinate=0.002, volume='annotation')
319+
ax = self.ba.plot_hslice(dv_coordinate=-0.002, volume='annotation')
320320
im = ax.get_images()[0]
321321
assert im.get_array().shape == (self.ba.bc.ny, self.ba.bc.nx, 3)
322322
ax.clear()
@@ -353,9 +353,9 @@ def test_slice(self):
353353
# tests output shapes
354354
self.assertTrue(ba.slice(axis=0, coordinate=0).shape == (ny, nz)) # sagittal
355355
self.assertTrue(ba.slice(axis=1, coordinate=0).shape == (nx, nz)) # coronal
356-
self.assertTrue(ba.slice(axis=2, coordinate=.002).shape == (ny, nx)) # horizontal
356+
self.assertTrue(ba.slice(axis=2, coordinate=-.002).shape == (ny, nx)) # horizontal
357357
# tests out of bound
358-
with self.assertRaises(IndexError):
358+
with self.assertRaises(ValueError):
359359
ba.slice(axis=1, coordinate=123)
360360
self.assertTrue(ba.slice(axis=1, coordinate=21, mode='clip').shape == (nx, nz))
361361
"""
@@ -553,9 +553,9 @@ def test_brain_coordinates(self):
553553
self.assertTrue(bc.ny == 7)
554554
self.assertTrue(bc.nz == 8)
555555
# test array functions
556-
in_out = [([6, 7, 8], np.array([6, 7, 8])),
557-
(np.array([6, 7, 8]), np.array([6, 7, 8])),
558-
(np.array([[6, 7, 8], [6, 7, 8]]), np.array([[6, 7, 8], [6, 7, 8]])),
556+
in_out = [([3, 4, 5], np.array([3, 4, 5])),
557+
(np.array([3, 4, 5]), np.array([3, 4, 5])),
558+
(np.array([[3, 4, 5], [3, 4, 5]]), np.array([[3, 4, 5], [3, 4, 5]])),
559559
]
560560
for io in in_out:
561561
self.assertTrue(np.all(bc.xyz2i(io[0]) == io[1]))

0 commit comments

Comments
 (0)