Coverage for ibllib/io/extractors/camera.py: 96%

1""" Camera extractor functions. 

2 

3This module handles extraction of camera timestamps for both Bpod and DAQ. 

4""" 

5import logging (empty)

6from functools import partial (empty)

7 

8import cv2 (empty)

9import numpy as np (empty)

10import matplotlib.pyplot as plt (empty)

11from iblutil.util import range_str (empty)

12 

13import ibldsp.utils as dsp (empty)

14from ibllib.plots import squares, vertical_lines (empty)

15from ibllib.io.video import assert_valid_label, VideoStreamer (empty)

16from iblutil.numerical import within_ranges (empty)

17from ibllib.io.extractors.base import get_session_extractor_type (empty)

18from ibllib.io.extractors.ephys_fpga import get_sync_fronts, get_sync_and_chn_map (empty)

19import ibllib.io.raw_data_loaders as raw (empty)

20import ibllib.io.extractors.video_motion as vmotion (empty)

21from ibllib.io.extractors.base import ( (empty)

22 BaseBpodTrialsExtractor, 

23 BaseExtractor, 

24 run_extractor_classes, 

25 _get_task_types_json_config 

26) 

27 

28_logger = logging.getLogger(__name__) (empty)

29 

30 

31def extract_camera_sync(sync, chmap=None): (empty)

32 """ 

33 Extract camera timestamps from the sync matrix. 

34 

35 :param sync: dictionary 'times', 'polarities' of fronts detected on sync trace 

36 :param chmap: dictionary containing channel indices. Default to constant. 

37 :return: dictionary containing camera timestamps 

38 """ 

39 assert chmap 9 ctx1qtpgwrefd

40 times = {} 9 ctx1qtpgwrefd

41 for k in filter(lambda x: x.endswith('_camera'), chmap): 9 ctx1qtpgwrefd

42 label, _ = k.rsplit('_', 1) 9 ctx1qtpgwrefd

43 times[label] = get_sync_fronts(sync, chmap[k]).times[::2] 9 ctx1qtpgwrefd

44 return times 9 ctx1qtpgwrefd

45 

46 

47def get_video_length(video_path): (empty)

48 """ 

49 Returns video length. 

50 

51 :param video_path: A path to the video 

52 :return: 

53 """ 

54 is_url = isinstance(video_path, str) and video_path.startswith('http') 16 ctx1lmqpkgcrefdxoban

55 cap = VideoStreamer(video_path).cap if is_url else cv2.VideoCapture(str(video_path)) 16 ctx1lmqpkgcrefdxoban

56 assert cap.isOpened(), f'Failed to open video file {video_path}' 16 ctx1lmqpkgcrefdxoban

57 length = int(cap.get(cv2.CAP_PROP_FRAME_COUNT)) 16 ctx1lmqpkgcrefdxoban

58 cap.release() 16 ctx1lmqpkgcrefdxoban

59 return length 16 ctx1lmqpkgcrefdxoban

60 

61 

62class CameraTimestampsFPGA(BaseExtractor): (empty)

63 """Extractor for videos using DAQ sync and channel map.""" 

64 

65 def __init__(self, label, session_path=None): (empty)

66 super().__init__(session_path) 7 ctx1qpgrefd

67 self.label = assert_valid_label(label) 7 ctx1qpgrefd

68 self.save_names = f'_ibl_{label}Camera.times.npy' 7 ctx1qpgrefd

69 self.var_names = f'{label}_camera_timestamps' 7 ctx1qpgrefd

70 self._log_level = _logger.level 7 ctx1qpgrefd

71 _logger.setLevel(logging.DEBUG) 7 ctx1qpgrefd

72 

73 def __del__(self): (empty)

74 _logger.setLevel(self._log_level) 6 ctx1sqpgef

75 

76 def _extract(self, sync=None, chmap=None, video_path=None, sync_label='audio', (empty)

77 display=False, extrapolate_missing=True, **kwargs): 

78 """ 

79 The raw timestamps are taken from the DAQ. These are the times of the camera's frame TTLs. 

80 If the pin state file exists, these timestamps are aligned to the video frames using 

81 task TTLs (typically the audio TTLs). Frames missing from the embedded frame count are 

82 removed from the timestamps array. 

83 If the pin state file does not exist, the left and right camera timestamps may be aligned 

84 using the wheel data. 

85 

86 Parameters 

87 ---------- 

88 sync : dict 

89 Dictionary 'times', 'polarities' of fronts detected on sync trace. 

90 chmap : dict 

91 Dictionary containing channel indices. Default to constant. 

92 video_path : str, pathlib.Path, int 

93 An optional path for fetching the number of frames. If None, the video is loaded from 

94 the session path. If an int is provided this is taken to be the total number of frames. 

95 sync_label : str 

96 The sync label of the channel that's wired to the GPIO for synchronising the times. 

97 display : bool 

98 If true, the TTL and GPIO fronts are plotted. 

99 extrapolate_missing : bool 

100 If true, any missing timestamps at the beginning and end of the session are 

101 extrapolated based on the median frame rate, otherwise they will be NaNs. 

102 **kwargs 

103 Extra keyword arguments (unused). 

104 

105 Returns 

106 ------- 

107 numpy.array 

108 The extracted camera timestamps. 

109 """ 

110 fpga_times = extract_camera_sync(sync=sync, chmap=chmap) 7 ctx1qpgrefd

111 count, (*_, gpio) = raw.load_embedded_frame_data(self.session_path, self.label) 7 ctx1qpgrefd

112 raw_ts = fpga_times[self.label] 7 ctx1qpgrefd

113 

114 if video_path is None: 7 ctx1qpgrefd

115 filename = f'_iblrig_{self.label}Camera.raw.mp4' 6 ctx1qprefd

116 video_path = self.session_path.joinpath('raw_video_data', filename) 6 ctx1qprefd

117 # Permit the video path to be the length for development and debugging purposes 

118 length = (video_path if isinstance(video_path, int) else get_video_length(video_path)) 7 ctx1qpgrefd

119 _logger.debug(f'Number of video frames = {length}') 7 ctx1qpgrefd

120 

121 if gpio is not None and gpio['indices'].size > 1 and sync_label is not None: 7 ctx1qpgrefd

122 _logger.info(f'Aligning to {sync_label} TTLs') 5 ctx1grefd

123 # Extract sync TTLs 

124 ttl = get_sync_fronts(sync, chmap[sync_label]) 5 ctx1grefd

125 _, ts = raw.load_camera_ssv_times(self.session_path, self.label) 5 ctx1grefd

126 try: 5 ctx1grefd

127 """ 5 ctx1grefd

128 NB: Some of the sync TTLs occur very close together, and are therefore not 

129 reflected in the pin state. This function removes those. Also converts frame 

130 times to DAQ time. 

131 """ 

132 gpio, ttl, ts = groom_pin_state(gpio, ttl, ts, display=display) 5 ctx1grefd

133 """ 4 ctx1gefd

134 The length of the count and pin state are regularly longer than the length of 

135 the video file. Here we assert that the video is either shorter or the same 

136 length as the arrays, and we make an assumption that the missing frames are 

137 right at the end of the video. We therefore simply shorten the arrays to match 

138 the length of the video. 

139 """ 

140 if count.size > length: 4 ctx1gefd

141 count = count[:length] 4 ctx1gefd

142 else: 

143 assert length == count.size, 'fewer counts than frames' 2 ctx1ef

144 assert raw_ts.shape[0] > 0, 'no timestamps found in channel indicated for ' \ 4 ctx1gefd

145 f'{self.label} camera' 

146 return align_with_gpio(raw_ts, ttl, gpio, count, 4 ctx1gefd

147 display=display, 

148 extrapolate_missing=extrapolate_missing) 

149 except AssertionError as ex: 1 ctx1r

150 _logger.critical('Failed to extract using %s: %s', sync_label, ex) 1 ctx1r

151 

152 # If you reach here extracting using sync TTLs was not possible, we attempt to align using wheel motion energy 

153 _logger.warning('Attempting to align using wheel') 3 ctx1qpr

154 

155 try: 3 ctx1qpr

156 if self.label not in ['left', 'right']: 3 ctx1qpr

157 # Can only use wheel alignment for left and right cameras 

158 raise ValueError(f'Wheel alignment not supported for {self.label} camera') 2 ctx1qp

159 

160 motion_class = vmotion.MotionAlignmentFullSession(self.session_path, self.label, sync='nidq', upload=True) 3 ctx1qpr

161 new_times = motion_class.process() 

162 if not motion_class.qc_outcome: 

163 raise ValueError(f'Wheel alignment for {self.label} camera failed to pass qc: {motion_class.qc}') 

164 else: 

165 _logger.warning(f'Wheel alignment for {self.label} camera successful, qc: {motion_class.qc}') 

166 return new_times 

167 

168 except Exception as err: 3 ctx1qpr

169 _logger.critical(f'Failed to align with wheel for {self.label} camera: {err}') 3 ctx1qpr

170 

171 if length < raw_ts.size: 3 ctx1qpr

172 df = raw_ts.size - length 2 ctx1qp

173 _logger.info(f'Discarding first {df} pulses') 2 ctx1qp

174 raw_ts = raw_ts[df:] 2 ctx1qp

175 

176 return raw_ts 3 ctx1qpr

177 

178 

179class CameraTimestampsCamlog(BaseExtractor): (empty)

180 def __init__(self, label, session_path=None): (empty)

181 super().__init__(session_path) 1 ctx1t

182 self.label = assert_valid_label(label) 1 ctx1t

183 self.save_names = f'_ibl_{label}Camera.times.npy' 1 ctx1t

184 self.var_names = f'{label}_camera_timestamps' 1 ctx1t

185 self._log_level = _logger.level 1 ctx1t

186 _logger.setLevel(logging.DEBUG) 1 ctx1t

187 

188 def __del__(self): (empty)

189 _logger.setLevel(self._log_level) 1 ctx1t

190 

191 def _extract(self, sync=None, chmap=None, video_path=None, (empty)

192 display=False, extrapolate_missing=True, **kwargs): 

193 

194 fpga_times = extract_camera_sync(sync=sync, chmap=chmap) 1 ctx1t

195 video_frames = get_video_length(self.session_path.joinpath('raw_video_data', f'_iblrig_{self.label}Camera.raw.mp4')) 1 ctx1t

196 raw_ts = fpga_times[self.label] 1 ctx1t

197 

198 # For left camera sometimes we have one extra pulse than video frame 

199 if (raw_ts.size - video_frames) == 1: 1 ctx1t

200 _logger.warning(f'One extra sync pulse detected for {self.label} camera') 

201 raw_ts = raw_ts[:-1] 

202 elif (raw_ts.size - video_frames) == -1: 1 ctx1t

203 _logger.warning(f'One extra video frame detected for {self.label} camera') 

204 med_time = np.median(np.diff(raw_ts)) 

205 raw_ts = np.r_[raw_ts, np.array([raw_ts[-1] + med_time])] 

206 

207 assert video_frames == raw_ts.size, f'dimension mismatch between video frames and TTL pulses for {self.label} camera' \ 1 ctx1t

208 f'by {np.abs(video_frames - raw_ts.size)} frames' 

209 

210 return raw_ts 1 ctx1t

211 

212 

213class CameraTimestampsBpod(BaseBpodTrialsExtractor): (empty)

214 """ 

215 Get the camera timestamps from the Bpod 

216 

217 The camera events are logged only during the events not in between, so the times need 

218 to be interpolated 

219 """ 

220 save_names = '_ibl_leftCamera.times.npy' (empty)

221 var_names = 'left_camera_timestamps' (empty)

222 

223 def __init__(self, *args, **kwargs): (empty)

224 super().__init__(*args, **kwargs) 8 ctx1lmkcoban

225 self._log_level = _logger.level 8 ctx1lmkcoban

226 _logger.setLevel(logging.DEBUG) 8 ctx1lmkcoban

227 

228 def __del__(self): (empty)

229 _logger.setLevel(self._log_level) 7 ctx1slmkcbn

230 

231 def _extract(self, video_path=None, display=False, extrapolate_missing=True, **kwargs): (empty)

232 """ 

233 The raw timestamps are taken from the Bpod. These are the times of the camera's frame TTLs. 

234 If the pin state file exists, these timestamps are aligned to the video frames using the 

235 sync TTLs. Frames missing from the embedded frame count are removed from the timestamps 

236 array. 

237 If the pin state file does not exist, the left camera timestamps may be aligned using the 

238 wheel data. 

239 :param video_path: an optional path for fetching the number of frames. If None, 

240 the video is loaded from the session path. If an int is provided this is taken to be 

241 the total number of frames. 

242 :param display: if True, the TTL and GPIO fronts are plotted. 

243 :param extrapolate_missing: if True, any missing timestamps at the beginning and end of 

244 the session are extrapolated based on the median frame rate, otherwise they will be NaNs. 

245 :return: a numpy array of camera timestamps 

246 """ 

247 raw_ts = self._times_from_bpod() 8 ctx1lmkcoban

248 count, (*_, gpio) = raw.load_embedded_frame_data(self.session_path, 'left') 8 ctx1lmkcoban

249 if video_path is None: 8 ctx1lmkcoban

250 filename = '_iblrig_leftCamera.raw.mp4' 6 ctx1lmoban

251 video_path = self.session_path.joinpath('raw_video_data', filename) 6 ctx1lmoban

252 # Permit the video path to be the length for development and debugging purposes 

253 length = video_path if isinstance(video_path, int) else get_video_length(video_path) 8 ctx1lmkcoban

254 _logger.debug(f'Number of video frames = {length}') 8 ctx1lmkcoban

255 

256 # Check if the GPIO is usable for extraction. GPIO is None if the file does not exist, 

257 # is empty, or contains only one value (i.e. doesn't change) 

258 if gpio is not None and gpio['indices'].size > 1: 8 ctx1lmkcoban

259 _logger.info('Aligning to sync TTLs') 3 ctx1cba

260 # Extract audio TTLs 

261 _, audio = raw.load_bpod_fronts(self.session_path, data=self.bpod_trials, 3 ctx1cba

262 task_collection=self.task_collection) 

263 _, ts = raw.load_camera_ssv_times(self.session_path, 'left') 3 ctx1cba

264 """ 3 ctx1cba

265 There are many sync TTLs that are for some reason missed by the GPIO. Conversely 

266 the last GPIO doesn't often correspond to any audio TTL. These will be removed. 

267 The drift appears to be less severe than the DAQ, so when assigning TTLs we'll take 

268 the nearest TTL within 500ms. The go cue TTLs comprise two short pulses ~3ms apart. 

269 We will fuse any TTLs less than 5ms apart to make assignment more accurate. 

270 """ 

271 try: 3 ctx1cba

272 gpio, audio, ts = groom_pin_state(gpio, audio, ts, take='nearest', 3 ctx1cba

273 tolerance=.5, min_diff=5e-3, display=display) 

274 if count.size > length: 3 ctx1cba

275 count = count[:length] 1 ctx1a

276 else: 

277 assert length == count.size, 'fewer counts than frames' 3 ctx1cba

278 

279 return align_with_gpio(raw_ts, audio, gpio, count, 3 ctx1cba

280 extrapolate_missing, display=display) 

281 except AssertionError as ex: 1 ctx1a

282 _logger.critical('Failed to extract using audio: %s', ex) 1 ctx1a

283 

284 # If you reach here extracting using audio TTLs was not possible 

285 _logger.warning('Alignment by wheel data not yet implemented') 6 ctx1lmkoan

286 # Extrapolate at median frame rate 

287 n_missing = length - raw_ts.size 6 ctx1lmkoan

288 if n_missing > 0: 6 ctx1lmkoan

289 _logger.warning(f'{n_missing} fewer Bpod timestamps than frames; ' 2 ctx1oa

290 f'{"extrapolating" if extrapolate_missing else "appending nans"}') 

291 frate = np.median(np.diff(raw_ts)) 2 ctx1oa

292 to_app = ((np.arange(n_missing, ) + 1) / frate + raw_ts[-1] 2 ctx1oa

293 if extrapolate_missing 

294 else np.full(n_missing, np.nan)) 

295 raw_ts = np.r_[raw_ts, to_app] # Append the missing times 2 ctx1oa

296 elif n_missing < 0: 5 ctx1lmkon

297 _logger.warning(f'{abs(n_missing)} fewer frames than Bpod timestamps') 4 ctx1lmon

298 _logger.info(f'Discarding first {abs(n_missing)} pulses') 4 ctx1lmon

299 raw_ts = raw_ts[abs(n_missing):] 4 ctx1lmon

300 

301 return raw_ts 6 ctx1lmkoan

302 

303 def _times_from_bpod(self): (empty)

304 ntrials = len(self.bpod_trials) 8 ctx1lmkcoban

305 

306 cam_times = [] 8 ctx1lmkcoban

307 n_frames = 0 8 ctx1lmkcoban

308 n_out_of_sync = 0 8 ctx1lmkcoban

309 missed_trials = [] 8 ctx1lmkcoban

310 for ind in np.arange(ntrials): 8 ctx1lmkcoban

311 # get upgoing and downgoing fronts 

312 pin = np.array(self.bpod_trials[ind]['behavior_data'] 8 ctx1lmkcoban

313 ['Events timestamps'].get('Port1In')) 

314 pout = np.array(self.bpod_trials[ind]['behavior_data'] 8 ctx1lmkcoban

315 ['Events timestamps'].get('Port1Out')) 

316 # some trials at startup may not have the camera working, discard 

317 if np.all(pin) is None: 8 ctx1lmkcoban

318 missed_trials.append(ind) 2 ctx1ka

319 continue 2 ctx1ka

320 # if the trial starts in the middle of a square, discard the first downgoing front 

321 if pout[0] < pin[0]: 8 ctx1lmkcoban

322 pout = pout[1:] 8 ctx1lmkcoban

323 # same if the last sample is during an upgoing front, 

324 # always put size as it happens last 

325 pin = pin[:pout.size] 8 ctx1lmkcoban

326 frate = np.median(np.diff(pin)) 8 ctx1lmkcoban

327 if ind > 0: 8 ctx1lmkcoban

328 """ 8 ctx1lmkcoban

329 assert that the pulses have the same length and that we don't miss frames during 

330 the trial, the refresh rate of bpod is 100us 

331 """ 

332 test1 = np.all(np.abs(1 - (pin - pout) / np.median(pin - pout)) < 0.1) 8 ctx1lmkcoban

333 test2 = np.all(np.abs(np.diff(pin) - frate) <= 0.00011) 8 ctx1lmkcoban

334 if not all([test1, test2]): 8 ctx1lmkcoban

335 n_out_of_sync += 1 4 ctx1coba

336 # grow a list of cam times for ech trial 

337 cam_times.append(pin) 8 ctx1lmkcoban

338 n_frames += pin.size 8 ctx1lmkcoban

339 

340 if missed_trials: 8 ctx1lmkcoban

341 _logger.debug('trial(s) %s missing TTL events', range_str(missed_trials)) 2 ctx1ka

342 if n_out_of_sync > 0: 8 ctx1lmkcoban

343 _logger.warning(f"{n_out_of_sync} trials with bpod camera frame times not within" 4 ctx1coba

344 f" 10% of the expected sampling rate") 

345 

346 t_first_frame = np.array([c[0] for c in cam_times]) 8 ctx1lmkcoban

347 t_last_frame = np.array([c[-1] for c in cam_times]) 8 ctx1lmkcoban

348 frate = 1 / np.nanmedian(np.array([np.median(np.diff(c)) for c in cam_times])) 8 ctx1lmkcoban

349 intertrial_duration = t_first_frame[1:] - t_last_frame[:-1] 8 ctx1lmkcoban

350 intertrial_missed_frames = np.int32(np.round(intertrial_duration * frate)) - 1 8 ctx1lmkcoban

351 

352 # initialize the full times array 

353 frame_times = np.zeros(n_frames + int(np.sum(intertrial_missed_frames))) 8 ctx1lmkcoban

354 ii = 0 8 ctx1lmkcoban

355 for trial, cam_time in enumerate(cam_times): 8 ctx1lmkcoban

356 if cam_time is not None: 8 ctx1lmkcoban

357 # populate first the recovered times within the trials 

358 frame_times[ii: ii + cam_time.size] = cam_time 8 ctx1lmkcoban

359 ii += cam_time.size 8 ctx1lmkcoban

360 if trial == (len(cam_times) - 1): 8 ctx1lmkcoban

361 break 8 ctx1lmkcoban

362 # then extrapolate in-between 

363 nmiss = intertrial_missed_frames[trial] 8 ctx1lmkcoban

364 frame_times[ii: ii + nmiss] = (cam_time[-1] + intertrial_duration[trial] / 8 ctx1lmkcoban

365 (nmiss + 1) * (np.arange(nmiss) + 1)) 

366 ii += nmiss 8 ctx1lmkcoban

367 assert all(np.diff(frame_times) > 0) # negative diffs implies a big problem 8 ctx1lmkcoban

368 return frame_times 8 ctx1lmkcoban

369 

370 

371def align_with_gpio(timestamps, ttl, pin_state, count, extrapolate_missing=True, display=False): (empty)

372 """ 

373 Groom the raw DAQ or Bpod camera timestamps using the frame embedded GPIO and frame counter. 

374 

375 Parameters 

376 ---------- 

377 timestamps : numpy.array 

378 An array of raw DAQ or Bpod camera timestamps. 

379 ttl : dict 

380 A dictionary of DAQ sync TTLs, with keys {'times', 'polarities'}. 

381 pin_state : dict 

382 A dictionary containing GPIO pin state values, with keys {'indices', 'polarities'}. 

383 count : numpy.array 

384 An array of frame numbers. 

385 extrapolate_missing : bool 

386 If true and the number of timestamps is fewer than the number of frame counts, the 

387 remaining timestamps are extrapolated based on the frame rate, otherwise they are NaNs. 

388 display : bool 

389 Plot the resulting timestamps. 

390 

391 Returns 

392 ------- 

393 numpy.array 

394 The corrected frame timestamps. 

395 """ 

396 # Some assertions made on the raw data 

397 # assert count.size == pin_state.size, 'frame count and pin state size mismatch' 

398 assert all(np.diff(count) > 0), 'frame count not strictly increasing' 7 ctx1gcefdba

399 assert all(np.diff(timestamps) > 0), 'DAQ/Bpod camera times not strictly increasing' 7 ctx1gcefdba

400 same_n_ttl = pin_state['times'].size == ttl['times'].size 7 ctx1gcefdba

401 assert same_n_ttl, 'more ttl TTLs detected on camera than TTLs sent' 7 ctx1gcefdba

402 

403 """Here we will ensure that the DAQ camera times match the number of video frames in 7 ctx1gcefdba

404 length. We will make the following assumptions: 

405 

406 1. The number of DAQ camera times is equal to or greater than the number of video frames. 

407 2. No TTLs were missed between the camera and DAQ. 

408 3. No pin states were missed by Bonsai. 

409 4 No pixel count data was missed by Bonsai. 

410 

411 In other words the count and pin state arrays accurately reflect the number of frames 

412 sent by the camera and should therefore be the same length, and the length of the frame 

413 counter should match the number of saved video frames. 

414 

415 The missing frame timestamps are removed in three stages: 

416 

417 1. Remove any timestamps that occurred before video frame acquisition in Bonsai. 

418 2. Remove any timestamps where the frame counter reported missing frames, i.e. remove the 

419 dropped frames which occurred throughout the session. 

420 3. Remove the trailing timestamps at the end of the session if the camera was turned off 

421 in the wrong order. 

422 """ 

423 # Align on first pin state change 

424 first_uptick = pin_state['indices'][0] 7 ctx1gcefdba

425 first_ttl = np.searchsorted(timestamps, ttl['times'][0]) 7 ctx1gcefdba

426 """Here we find up to which index in the DAQ times we discard by taking the difference 7 ctx1gcefdba

427 between the index of the first pin state change (when the sync TTL was reported by the 

428 camera) and the index of the first sync TTL in DAQ time. We subtract the difference 

429 between the frame count at the first pin state change and the index to account for any 

430 video frames that were not saved during this period (we will remove those from the 

431 camera DAQ times later). 

432 """ 

433 # Minus any frames that were dropped between the start of frame acquisition and the first TTL 

434 start = first_ttl - first_uptick - (count[first_uptick] - first_uptick) 7 ctx1gcefdba

435 # Get approximate frame rate for extrapolating timestamps (if required) 

436 frate = round(1 / np.nanmedian(np.diff(timestamps))) 7 ctx1gcefdba

437 

438 if start < 0: 7 ctx1gcefdba

439 n_missing = abs(start) 3 ctx1cba

440 _logger.warning(f'{n_missing} missing DAQ/Bpod timestamp(s) at start; ' 3 ctx1cba

441 f'{"extrapolating" if extrapolate_missing else "prepending nans"}') 

442 to_app = (timestamps[0] - (np.arange(n_missing, 0, -1) + 1) / frate 3 ctx1cba

443 if extrapolate_missing 

444 else np.full(n_missing, np.nan)) 

445 timestamps = np.r_[to_app, timestamps] # Prepend the missing times 3 ctx1cba

446 start = 0 3 ctx1cba

447 

448 # Remove the extraneous timestamps from the beginning and end 

449 end = count[-1] + 1 + start 7 ctx1gcefdba

450 ts = timestamps[start:end] 7 ctx1gcefdba

451 if (n_missing := count[-1] - ts.size + 1) > 0: 7 ctx1gcefdba

452 """ 7 ctx1gcefdba

453 For ephys sessions there may be fewer DAQ times than frame counts if DAQ acquisition is 

454 turned off before the video acquisition workflow. For Bpod this always occurs because Bpod 

455 finishes before the camera workflow. For Bpod the times are already extrapolated for 

456 these late frames.""" 

457 _logger.warning(f'{n_missing} fewer DAQ/Bpod timestamps than frame counts; ' 7 ctx1gcefdba

458 f'{"extrapolating" if extrapolate_missing else "appending nans"}') 

459 to_app = ((np.arange(n_missing, ) + 1) / frate + ts[-1] 7 ctx1gcefdba

460 if extrapolate_missing 

461 else np.full(n_missing, np.nan)) 

462 ts = np.r_[ts, to_app] # Append the missing times 7 ctx1gcefdba

463 assert ts.size >= count.size, 'fewer timestamps than frame counts' 7 ctx1gcefdba

464 assert ts.size == count[-1] + 1, 'more frames recorded in frame count than timestamps ' 7 ctx1gcefdba

465 

466 # Remove the rest of the dropped frames 

467 ts = ts[count] 7 ctx1gcefdba

468 assert np.searchsorted(ts, ttl['times'][0]) == first_uptick, \ 7 ctx1gcefdba

469 'time of first sync TTL doesn\'t match after alignment' 

470 if ts.size != count.size: 7 ctx1gcefdba

471 _logger.error('number of timestamps and frames don\'t match after alignment') 

472 

473 if display: 7 ctx1gcefdba

474 # Plot to check 

475 fig, axes = plt.subplots(1, 1) 2 ctx1da

476 y = within_ranges(np.arange(ts.size), pin_state['indices'].reshape(-1, 2)).astype(float) 2 ctx1da

477 y *= 1e-5 # For scale when zoomed in 2 ctx1da

478 axes.plot(ts, y, marker='d', color='blue', drawstyle='steps-pre', label='GPIO') 2 ctx1da

479 axes.plot(ts, np.zeros_like(ts), 'kx', label='DAQ timestamps') 2 ctx1da

480 vertical_lines(ttl['times'], ymin=0, ymax=1e-5, 2 ctx1da

481 color='r', linestyle=':', ax=axes, label='sync TTL') 

482 plt.legend() 2 ctx1da

483 

484 return ts 7 ctx1gcefdba

485 

486 

487def attribute_times(arr, events, tol=.1, injective=True, take='first'): (empty)

488 """ 

489 Returns the values of the first array that correspond to those of the second. 

490 

491 Given two arrays of timestamps, the function will return the values of the first array 

492 that most likely correspond to the values of the second. For each of the values in the 

493 second array, the absolute difference is taken and the index of either the first sufficiently 

494 close value, or simply the closest one, is assigned. 

495 

496 If injective is True, once a value has been assigned to an event it can't be assigned to 

497 another. In other words there is a one-to-one mapping between the two arrays. 

498 

499 Parameters 

500 ---------- 

501 arr : numpy.array 

502 An array of event times to attribute to those in `events`. 

503 events : numpy.array 

504 An array of event times considered a subset of `arr`. 

505 tol : float 

506 The max absolute difference between values in order to be considered a match. 

507 injective : bool 

508 If true, once a value has been assigned it will not be assigned again. 

509 take : {'first', 'nearest', 'after'} 

510 If 'first' the first value within tolerance is assigned; if 'nearest' the 

511 closest value is assigned; if 'after' assign the first event after. 

512 

513 Returns 

514 ------- 

515 numpy.array 

516 An array the same length as `events` containing indices of `arr` corresponding to each 

517 event. 

518 """ 

519 if (take := take.lower()) not in ('first', 'nearest', 'after'): 12 ctx1uivgcefdjbah

520 raise ValueError('Parameter `take` must be either "first", "nearest", or "after"') 1 ctx1u

521 stack = np.ma.masked_invalid(arr, copy=False) 12 ctx1uivgcefdjbah

522 stack.fill_value = np.inf 12 ctx1uivgcefdjbah

523 # If there are no invalid values, the mask is False so let's ensure it's a bool array 

524 if stack.mask is np.bool_(0): 12 ctx1uivgcefdjbah

525 stack.mask = np.zeros(arr.shape, dtype=bool) 

526 assigned = np.full(events.shape, -1, dtype=int) # Initialize output array 12 ctx1uivgcefdjbah

527 min_tol = 0 if take == 'after' else -tol 12 ctx1uivgcefdjbah

528 for i, x in enumerate(events): 12 ctx1uivgcefdjbah

529 dx = stack.filled() - x 12 ctx1uivgcefdjbah

530 candidates = np.logical_and(min_tol < dx, dx < tol) 12 ctx1uivgcefdjbah

531 if any(candidates): # is any value within tolerance 12 ctx1uivgcefdjbah

532 idx = np.abs(dx).argmin() if take == 'nearest' else np.where(candidates)[0][0] 12 ctx1uivgcefdjbah

533 assigned[i] = idx 12 ctx1uivgcefdjbah

534 stack.mask[idx] = injective # If one-to-one, remove the assigned value 12 ctx1uivgcefdjbah

535 return assigned 12 ctx1uivgcefdjbah

536 

537 

538def groom_pin_state(gpio, ttl, ts, tolerance=2., display=False, take='first', min_diff=0.): (empty)

539 """ 

540 Align the GPIO pin state to the DAQ sync TTLs. Any sync TTLs not reflected in the pin 

541 state are removed from the dict and the times of the detected fronts are converted to DAQ 

542 time. At the end of this the number of GPIO fronts should equal the number of TTLs. 

543 

544 Note: 

545 - This function is ultra safe: we probably don't need assign all the ups and down fronts. 

546 separately and could potentially even align the timestamps without removing the missed fronts 

547 - The input gpio and TTL dicts may be modified by this function. 

548 - For training sessions the frame rate is only 30Hz and the TTLs tend to be broken up by 

549 small gaps. Setting the min_diff to 5ms helps the timestamp assignment accuracy. 

550 

551 Parameters 

552 ---------- 

553 gpio : dict 

554 A dictionary containing GPIO pin state values, with keys {'indices', 'polarities'}. 

555 ttl : dict 

556 A dictionary of DAQ sync TTLs, with keys {'times', 'polarities'}. 

557 ts : numpy.array 

558 The camera frame times (the camera frame TTLs acquired by the main DAQ). 

559 tolerance : float 

560 Two pulses need to be within this many seconds to be considered related. 

561 display : bool 

562 If true, the resulting timestamps are plotted against the raw audio signal. 

563 take : {'first', 'nearest'} 

564 If 'first' the first value within tolerance is assigned; if 'nearest' the 

565 closest value is assigned. 

566 min_diff : float 

567 Sync TTL fronts less than min_diff seconds apart will be removed. 

568 

569 Returns 

570 ------- 

571 dict 

572 Dictionary of GPIO DAQ front indices, polarities and DAQ aligned times. 

573 dict 

574 Sync TTL times and polarities sans the TTLs not detected in the frame data. 

575 numpy.array 

576 Frame times in DAQ time. 

577 

578 See Also 

579 -------- 

580 ibllib.io.extractors.ephys_fpga._get_sync_fronts 

581 """ 

582 # Check that the dimensions match 

583 if np.any(gpio['indices'] >= ts.size): 11 ctx1igcrefdjbah

584 _logger.warning('GPIO events occurring beyond timestamps array length') 1 ctx1j

585 keep = gpio['indices'] < ts.size 1 ctx1j

586 gpio = {k: gpio[k][keep] for k, v in gpio.items()} 1 ctx1j

587 assert ttl and ttl['times'].size > 0, 'no sync TTLs for session' 11 ctx1igcrefdjbah

588 assert ttl['times'].size == ttl['polarities'].size, 'sync TTL data dimension mismatch' 11 ctx1igcrefdjbah

589 # make sure that there are no 2 consecutive fall or consecutive rise events 

590 assert np.all(np.abs(np.diff(ttl['polarities'])) == 2), 'consecutive high/low sync TTL events' 11 ctx1igcrefdjbah

591 # make sure first TTL is high 

592 assert ttl['polarities'][0] == 1 11 ctx1igcrefdjbah

593 # make sure ttl times in order 

594 assert np.all(np.diff(ttl['times']) > 0) 11 ctx1igcrefdjbah

595 # make sure raw timestamps increase 

596 assert np.all(np.diff(ts) > 0), 'timestamps must strictly increase' 11 ctx1igcrefdjbah

597 # make sure there are state changes 

598 assert gpio['indices'].any(), 'no TTLs detected in GPIO' 11 ctx1igcrefdjbah

599 # # make sure first GPIO state is high 

600 assert gpio['polarities'][0] == 1 11 ctx1igcrefdjbah

601 """ 10 ctx1igcefdjbah

602 Some sync TTLs appear to be so short that they are not recorded by the camera. These can 

603 be as short as a few microseconds. Applying a cutoff based on framerate was unsuccessful. 

604 Assigning each sync TTL to each pin state change is not easy because some onsets occur very 

605 close together (sometimes < 70ms), on the order of the delay between TTL and frame time. 

606 Also, the two clocks have some degree of drift, so the delay between sync TTL and pin state 

607 change may be zero or even negative. 

608 

609 Here we split the events into sync TTL onsets (lo->hi) and TTL offsets (hi->lo). For each 

610 uptick in the GPIO pin state, we take the first TTL onset time that was within 100ms of it. 

611 We ensure that each sync TTL is assigned only once, so a TTL that is closer to frame 3 than 

612 frame 1 may still be assigned to frame 1. 

613 """ 

614 ifronts = gpio['indices'] # The pin state flips 10 ctx1igcefdjbah

615 sync_times = ttl['times'] 10 ctx1igcefdjbah

616 if ifronts.size != ttl['times'].size: 10 ctx1igcefdjbah

617 _logger.warning('more sync TTLs than GPIO state changes, assigning timestamps') 10 ctx1igcefdjbah

618 to_remove = np.zeros(ifronts.size, dtype=bool) # unassigned GPIO fronts to remove 10 ctx1igcefdjbah

619 low2high = ifronts[gpio['polarities'] == 1] 10 ctx1igcefdjbah

620 high2low = ifronts[gpio['polarities'] == -1] 10 ctx1igcefdjbah

621 assert low2high.size >= high2low.size 10 ctx1igcefdjbah

622 

623 # Remove and/or fuse short TTLs 

624 if min_diff > 0: 10 ctx1igcefdjbah

625 short, = np.where(np.diff(ttl['times']) < min_diff) 5 ctx1icbah

626 sync_times = np.delete(ttl['times'], np.r_[short, short + 1]) 5 ctx1icbah

627 _logger.debug(f'Removed {short.size * 2} fronts TLLs less than ' 5 ctx1icbah

628 f'{min_diff * 1e3:.0f}ms apart') 

629 assert sync_times.size > 0, f'all sync TTLs less than {min_diff}s' 5 ctx1icbah

630 

631 # Onsets 

632 ups = ts[low2high] - ts[low2high][0] # times relative to first GPIO high 10 ctx1igcefdjbah

633 onsets = sync_times[::2] - sync_times[0] # TTL times relative to first onset 10 ctx1igcefdjbah

634 # assign GPIO fronts to ttl onset 

635 assigned = attribute_times(onsets, ups, tol=tolerance, take=take) 10 ctx1igcefdjbah

636 unassigned = np.setdiff1d(np.arange(onsets.size), assigned[assigned > -1]) 10 ctx1igcefdjbah

637 if unassigned.size > 0: 10 ctx1igcefdjbah

638 _logger.debug(f'{unassigned.size} sync TTL rises were not detected by the camera') 10 ctx1igcefdjbah

639 # Check that all pin state upticks could be attributed to an onset TTL 

640 if np.any(missed := assigned == -1): 10 ctx1igcefdjbah

641 _logger.warning(f'{sum(missed)} pin state rises could not be attributed to a sync TTL') 6 ctx1icjbah

642 if display: 6 ctx1icjbah

643 ax = plt.subplot() 1 ctx1a

644 vertical_lines(ups[assigned > -1], 1 ctx1a

645 linestyle='-', color='g', ax=ax, 

646 label='assigned GPIO up state') 

647 vertical_lines(ups[missed], 1 ctx1a

648 linestyle='-', color='r', ax=ax, 

649 label='unassigned GPIO up state') 

650 vertical_lines(onsets[unassigned], 1 ctx1a

651 linestyle=':', color='k', ax=ax, 

652 alpha=0.3, label='sync TTL onset') 

653 vertical_lines(onsets[assigned], 1 ctx1a

654 linestyle=':', color='b', ax=ax, label='assigned TTL onset') 

655 plt.legend() 1 ctx1a

656 plt.show() 1 ctx1a

657 # Remove the missed fronts 

658 to_remove = np.in1d(gpio['indices'], low2high[missed]) 6 ctx1icjbah

659 assigned = assigned[~missed] 6 ctx1icjbah

660 onsets_ = sync_times[::2][assigned] 10 ctx1igcefdjbah

661 

662 # Offsets 

663 downs = ts[high2low] - ts[high2low][0] 10 ctx1igcefdjbah

664 offsets = sync_times[1::2] - sync_times[1] 10 ctx1igcefdjbah

665 assigned = attribute_times(offsets, downs, tol=tolerance, take=take) 10 ctx1igcefdjbah

666 unassigned = np.setdiff1d(np.arange(offsets.size), assigned[assigned > -1]) 10 ctx1igcefdjbah

667 if unassigned.size > 0: 10 ctx1igcefdjbah

668 _logger.debug(f'{unassigned.size} sync TTL falls were not detected by the camera') 9 ctx1gcefdjbah

669 # Check that all pin state downticks could be attributed to an offset TTL 

670 if np.any(missed := assigned == -1): 10 ctx1igcefdjbah

671 _logger.warning(f'{sum(missed)} pin state falls could not be attributed to a sync TTL') 5 ctx1cjbah

672 # Remove the missed fronts 

673 to_remove |= np.in1d(gpio['indices'], high2low[missed]) 5 ctx1cjbah

674 assigned = assigned[~missed] 5 ctx1cjbah

675 offsets_ = sync_times[1::2][assigned] 10 ctx1igcefdjbah

676 

677 # Sync TTLs groomed 

678 if np.any(to_remove): 10 ctx1igcefdjbah

679 # Check for any orphaned fronts (only one pin state edge was assigned) 

680 to_remove = np.pad(to_remove, (0, to_remove.size % 2), 'edge') # Ensure even size 6 ctx1icjbah

681 # Perform xor to find GPIOs where only onset or offset is marked for removal 

682 orphaned = to_remove.reshape(-1, 2).sum(axis=1) == 1 6 ctx1icjbah

683 if orphaned.any(): 6 ctx1icjbah

684 """If there are orphaned GPIO fronts (i.e. only one edge was assigned to a sync 2 ctx1ih

685 TTL front), remove the orphaned front its assigned sync TTL. In other words 

686 if both edges cannot be assigned to a sync TTL, we ignore the TTL entirely. 

687 This is a sign that the assignment was bad and extraction may fail.""" 

688 _logger.warning('Some onsets but not offsets (or vice versa) were not assigned; ' 2 ctx1ih

689 'this may be a sign of faulty wiring or clock drift') 

690 # Find indices of GPIO upticks where only the downtick was marked for removal 

691 orphaned_onsets, = np.where(~to_remove.reshape(-1, 2)[:, 0] & orphaned) 2 ctx1ih

692 # The onsets_ array already has the other TTLs removed (same size as to_remove == 

693 # False) so subtract the number of removed elements from index. 

694 for i, v in enumerate(orphaned_onsets): 2 ctx1ih

695 orphaned_onsets[i] -= to_remove.reshape(-1, 2)[:v, 0].sum() 1 ctx1h

696 # Same for offsets... 

697 orphaned_offsets, = np.where(~to_remove.reshape(-1, 2)[:, 1] & orphaned) 2 ctx1ih

698 for i, v in enumerate(orphaned_offsets): 2 ctx1ih

699 orphaned_offsets[i] -= to_remove.reshape(-1, 2)[:v, 1].sum() 2 ctx1ih