Coverage for ibllib/pipes/mesoscope_tasks.py: 90%

669 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-07-02 18:55 +0100

1"""The mesoscope data extraction pipeline. 

2 

3The mesoscope pipeline currently comprises raw image registration and timestamps extraction. In 

4the future there will be compression (and potential cropping), FOV metadata extraction, and ROI 

5extraction. 

6 

7Pipeline: 

8 1. Register reference images and upload snapshots and notes to Alyx 

9 2. Run ROI cell detection 

10 3. Calculate the pixel and ROI brain locations and register fields of view to Alyx 

11 4. Compress the raw imaging data 

12 5. Extract the imaging times from the main DAQ 

13""" 

14import json 

15import logging 

16import subprocess 

17from zipfile import ZipFile, ZIP_DEFLATED 

18import uuid 

19from pathlib import Path 

20from itertools import chain, groupby 

21from collections import Counter 

22from fnmatch import fnmatch 

23import enum 

24import re 

25import time 

26 

27import numba as nb 

28import numpy as np 

29import pandas as pd 

30import sparse 

31from scipy.io import loadmat 

32from scipy.interpolate import interpn 

33import one.alf.io as alfio 

34from one.alf.spec import to_alf 

35from one.alf.path import filename_parts, session_path_parts 

36import one.alf.exceptions as alferr 

37from iblutil.util import flatten, ensure_list 

38from iblatlas.atlas import ALLEN_CCF_LANDMARKS_MLAPDV_UM, MRITorontoAtlas 

39 

40from ibllib.pipes import base_tasks 

41from ibllib.oneibl.data_handlers import ExpectedDataset, dataset_from_name 

42from ibllib.io.extractors import mesoscope 

43 

44 

45_logger = logging.getLogger(__name__) 

46Provenance = enum.Enum('Provenance', ['ESTIMATE', 'FUNCTIONAL', 'LANDMARK', 'HISTOLOGY']) # py3.11 make StrEnum 

47 

48 

49class MesoscopeRegisterSnapshots(base_tasks.MesoscopeTask, base_tasks.RegisterRawDataTask): 

50 """Upload snapshots as Alyx notes and register the 2P reference image(s).""" 

51 priority = 100 

52 job_size = 'small' 

53 

54 @property 

55 def signature(self): 

56 I = ExpectedDataset.input # noqa 1tq

57 signature = { 1tq

58 'input_files': [I('referenceImage.raw.tif', f'{self.device_collection}/reference', False, register=True), 

59 I('referenceImage.stack.tif', f'{self.device_collection}/reference', False, register=True), 

60 I('referenceImage.meta.json', f'{self.device_collection}/reference', False, register=True)], 

61 'output_files': [] 

62 } 

63 return signature 1tq

64 

65 def __init__(self, session_path, **kwargs): 

66 super().__init__(session_path, **kwargs) 1utq

67 self.device_collection = self.get_device_collection('mesoscope', 1utq

68 kwargs.get('device_collection', 'raw_imaging_data_??')) 

69 

70 def _run(self): 

71 """ 

72 Assert one reference image per collection and rename it. Register snapshots. 

73 

74 Returns 

75 ------- 

76 list of pathlib.Path containing renamed reference image. 

77 """ 

78 # Assert that only one tif file exists per collection 

79 dsets = dataset_from_name('referenceImage.raw.tif', self.input_files) 1q

80 reference_images = list(chain.from_iterable(map(lambda x: x.find_files(self.session_path)[1], dsets))) 1q

81 assert len(set(x.parent for x in reference_images)) == len(reference_images) 1q

82 # Rename the reference images 

83 out_files = super()._run() 1q

84 # Register snapshots in base session folder and raw_imaging_data folders 

85 self.register_snapshots(collection=[self.device_collection, '']) 1q

86 return out_files 1q

87 

88 

89class MesoscopeCompress(base_tasks.MesoscopeTask): 

90 """ Tar compress raw 2p tif files, optionally remove uncompressed data.""" 

91 

92 priority = 90 

93 job_size = 'large' 

94 _log_level = None 

95 

96 @property 

97 def signature(self): 

98 signature = { 1f

99 'input_files': [('*.tif', self.device_collection, True)], 

100 'output_files': [('imaging.frames.tar.bz2', self.device_collection, True)] 

101 } 

102 return signature 1f

103 

104 def setUp(self, **kwargs): 

105 """Run at higher log level""" 

106 self._log_level = _logger.level 1f

107 _logger.setLevel(logging.DEBUG) 1f

108 return super().setUp(**kwargs) 1f

109 

110 def tearDown(self): 

111 _logger.setLevel(self._log_level or logging.INFO) 1f

112 return super().tearDown() 1f

113 

114 def _run(self, remove_uncompressed=False, verify_output=True, overwrite=False, **kwargs): 

115 """ 

116 Run tar compression on all tif files in the device collection. 

117 

118 Parameters 

119 ---------- 

120 remove_uncompressed: bool 

121 Whether to remove the original, uncompressed data. Default is False. 

122 verify_output: bool 

123 Whether to check that the compressed tar file can be uncompressed without errors. 

124 Default is True. 

125 

126 Returns 

127 ------- 

128 list of pathlib.Path 

129 Path to compressed tar file. 

130 """ 

131 outfiles = [] # should be one per raw_imaging_data folder 1f

132 _, all_tifs, _ = zip(*(x.find_files(self.session_path) for x in self.input_files)) 1f

133 if self.input_files[0].operator: # multiple device collections 1f

134 output_identifiers = self.output_files[0].identifiers 

135 # Check that the number of input ollections and output files match 

136 assert len(self.input_files[0].identifiers) == len(output_identifiers) 

137 else: 

138 output_identifiers = [self.output_files[0].identifiers] 1f

139 assert self.output_files[0].operator is None, 'only one output file expected' 1f

140 

141 # A list of tifs, grouped by raw imaging data collection 

142 input_files = groupby(chain.from_iterable(all_tifs), key=lambda x: x.parent) 1f

143 for (in_dir, infiles), out_id in zip(input_files, output_identifiers): 1f

144 infiles = list(infiles) 1f

145 outfile = self.session_path.joinpath(*filter(None, out_id)) 1f

146 if outfile.exists() and not overwrite: 1f

147 _logger.info('%s already exists; skipping...', outfile.relative_to(self.session_path)) 

148 outfiles.append(outfile) 

149 else: 

150 if not infiles: 1f

151 _logger.info('No image files found in %s', in_dir.relative_to(self.session_path)) 

152 continue 

153 

154 _logger.debug( 1f

155 'Input files:\n\t%s', '\n\t'.join(map(Path.as_posix, (x.relative_to(self.session_path) for x in infiles))) 

156 ) 

157 

158 uncompressed_size = sum(x.stat().st_size for x in infiles) 1f

159 _logger.info('Compressing %i file(s)', len(infiles)) 1f

160 cmd = 'tar -cjvf "{output}" "{input}"'.format( 1f

161 output=outfile.relative_to(in_dir), input='" "'.join(str(x.relative_to(in_dir)) for x in infiles)) 

162 _logger.debug(cmd) 1f

163 process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=in_dir) 1f

164 info, error = process.communicate() # b'2023-02-17_2_test_2P_00001_00001.tif\n' 1f

165 _logger.debug(info.decode()) 1f

166 assert process.returncode == 0, f'compression failed: {error.decode()}' 1f

167 

168 # Check the output 

169 assert outfile.exists(), 'output file missing' 1f

170 outfiles.append(outfile) 1f

171 compressed_size = outfile.stat().st_size 1f

172 min_size = kwargs.pop('verify_min_size', 1024) 1f

173 assert compressed_size > int(min_size), f'Compressed file < {min_size / 1024:.0f}KB' 1f

174 _logger.info('Compression ratio = %.3f, saving %.2f pct (%.2f MB)', 1f

175 uncompressed_size / compressed_size, 

176 round((1 - (compressed_size / uncompressed_size)) * 10000) / 100, 

177 (uncompressed_size - compressed_size) / 1024 / 1024) 

178 

179 if verify_output: 1f

180 # Test bzip 

181 cmd = f'bzip2 -tv {outfile.relative_to(in_dir)}' 1f

182 process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=in_dir) 1f

183 info, error = process.communicate() 1f

184 _logger.debug(info.decode()) 1f

185 assert process.returncode == 0, f'bzip compression test failed: {error}' 1f

186 # Check tar 

187 cmd = f'bunzip2 -dc {outfile.relative_to(in_dir)} | tar -tvf -' 1f

188 process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=in_dir) 1f

189 info, error = process.communicate() 1f

190 _logger.debug(info.decode()) 1f

191 assert process.returncode == 0, 'tarball decompression test failed' 1f

192 compressed_files = set(x.split()[-1] for x in filter(None, info.decode().split('\n'))) 1f

193 assert compressed_files == set(x.name for x in infiles) 1f

194 

195 if remove_uncompressed: 1f

196 _logger.info(f'Removing input files for {in_dir.relative_to(self.session_path)}') 1f

197 for file in infiles: 1f

198 file.unlink() 1f

199 

200 return outfiles 1f

201 

202 

203class MesoscopePreprocess(base_tasks.MesoscopeTask): 

204 """Run suite2p preprocessing on tif files.""" 

205 

206 priority = 80 

207 cpu = -1 

208 job_size = 'large' 

209 env = 'suite2p' 

210 

211 def __init__(self, *args, **kwargs): 

212 self._teardown_files = [] 1buaed

213 self.overwrite = False 1buaed

214 super().__init__(*args, **kwargs) 1buaed

215 

216 def setUp(self, **kwargs): 

217 """Set up task. 

218 

219 This will check the local filesystem for the raw tif files and if not present, will assume 

220 they have been compressed and deleted, in which case the signature will be replaced with 

221 the compressed input. 

222 

223 Note: this will not work correctly if only some collections have compressed tifs. 

224 """ 

225 self.overwrite = kwargs.get('overwrite', False) 1hm

226 all_files_present = super().setUp(**kwargs) # Ensure files present 1hm

227 bin_sig, = dataset_from_name('data.bin', self.input_files) 1hm

228 renamed_bin_sig, = dataset_from_name('imaging.frames_motionRegistered.bin', self.input_files) 1hm

229 if not self.overwrite and (bin_sig | renamed_bin_sig).find_files(self.session_path)[0]: 1hm

230 return all_files_present # We have local bin files; no need to extract tifs 

231 tif_sig = dataset_from_name('*.tif', self.input_files) 1hm

232 if not tif_sig: 1hm

233 return all_files_present # No tifs in the signature; just return 

234 tif_sig = tif_sig[0] 1hm

235 tifs_present, *_ = tif_sig.find_files(self.session_path) 1hm

236 if tifs_present or not all_files_present: 1hm

237 return all_files_present # Tifs present on disk; no need to decompress 1m

238 # Decompress imaging files 

239 tif_sigs = dataset_from_name('imaging.frames.tar.bz2', self.input_files) 1h

240 present, files, _ = zip(*(x.find_files(self.session_path) for x in tif_sigs)) 1h

241 if not all(present): 1h

242 return False # Compressed files missing; return 

243 files = flatten(files) 1h

244 _logger.info('Decompressing %i file(s)', len(files)) 1h

245 for file in files: 1h

246 cmd = 'tar -xvjf "{input}"'.format(input=file.name) 1h

247 _logger.debug(cmd) 1h

248 process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, cwd=file.parent) 1h

249 stdout, _ = process.communicate() # b'x 2023-02-17_2_test_2P_00001_00001.tif\n' 1h

250 _logger.debug(stdout.decode()) 1h

251 tifs = [file.parent.joinpath(x.split()[-1]) for x in stdout.decode().splitlines() if x.endswith('.tif')] 1h

252 assert process.returncode == 0 and len(tifs) > 0 1h

253 assert all(map(Path.exists, tifs)) 1h

254 self._teardown_files.extend(tifs) 1h

255 return all_files_present 1h

256 

257 def tearDown(self): 

258 """Tear down task. 

259 

260 This removes any decompressed tif files. 

261 """ 

262 for file in self._teardown_files: 1hm

263 _logger.debug('Removing %s', file) 1h

264 file.unlink() 1h

265 return super().tearDown() 1hm

266 

267 @property 

268 def signature(self): 

269 # The number of in and outputs will be dependent on the number of input raw imaging folders and output FOVs 

270 I = ExpectedDataset.input # noqa 1hmaed

271 signature = { 1hmaed

272 'input_files': [('_ibl_rawImagingData.meta.json', self.device_collection, True), 

273 I('*.tif', self.device_collection, True) | 

274 I('imaging.frames.tar.bz2', self.device_collection, True, unique=False), 

275 ('exptQC.mat', self.device_collection, False)], 

276 'output_files': [('mpci.ROIActivityF.npy', 'alf/FOV*', True), 

277 ('mpci.ROINeuropilActivityF.npy', 'alf/FOV*', True), 

278 ('mpci.ROIActivityDeconvolved.npy', 'alf/FOV*', True), 

279 ('mpci.badFrames.npy', 'alf/FOV*', True), 

280 ('mpci.mpciFrameQC.npy', 'alf/FOV*', True), 

281 ('mpciFrameQC.names.tsv', 'alf/FOV*', True), 

282 ('mpciMeanImage.images.npy', 'alf/FOV*', True), 

283 ('mpciROIs.stackPos.npy', 'alf/FOV*', True), 

284 ('mpciROIs.mpciROITypes.npy', 'alf/FOV*', True), 

285 ('mpciROIs.cellClassifier.npy', 'alf/FOV*', True), 

286 ('mpciROIs.uuids.csv', 'alf/FOV*', True), 

287 ('mpciROITypes.names.tsv', 'alf/FOV*', True), 

288 ('mpciROIs.masks.sparse_npz', 'alf/FOV*', True), 

289 ('mpciROIs.neuropilMasks.sparse_npz', 'alf/FOV*', True), 

290 ('_suite2p_ROIData.raw.zip', 'alf/FOV*', False), 

291 ('imaging.frames_motionRegistered.bin', 'suite2p/plane*', False)] 

292 } 

293 if not self.overwrite: # If not forcing re-registration, check whether bin files already exist on disk 1hmaed

294 # Including the data.bin in the expected signature ensures raw data files are not needlessly re-downloaded 

295 # and/or uncompressed during task setup as the local data.bin may be used instead 

296 # NB: The data.bin file is renamed to imaging.frames_motionRegistered.bin before registration to Alyx 

297 registered_bin = (I('data.bin', 'suite2p/plane*', True, unique=False) | 1hmaed

298 I('imaging.frames_motionRegistered.bin', 'suite2p/plane*', True, unique=False)) 

299 signature['input_files'][1] = registered_bin | signature['input_files'][1] 1hmaed

300 

301 return signature 1hmaed

302 

303 @staticmethod 

304 def _masks2sparse(stat, ops): 

305 """ 

306 Extract 3D sparse mask arrays from suit2p output. 

307 

308 Parameters 

309 ---------- 

310 stat : numpy.array 

311 The loaded stat.npy file. A structured array with fields ('lam', 'ypix', 'xpix', 'neuropil_mask'). 

312 ops : numpy.array 

313 The loaded ops.npy file. A structured array with fields ('Ly', 'Lx'). 

314 

315 Returns 

316 ------- 

317 sparse.GCXS 

318 A pydata sparse array of type float32, representing the ROI masks. 

319 sparse.GCXS 

320 A pydata sparse array of type float32, representing the neuropil ROI masks. 

321 

322 Notes 

323 ----- 

324 Save using sparse.save_npz. 

325 """ 

326 shape = (stat.shape[0], ops['Ly'], ops['Lx']) 1ed

327 npx = np.prod(shape[1:]) # Number of pixels per time point 1ed

328 coords = [[], [], []] 1ed

329 data = [] 1ed

330 pil_coords = [] 1ed

331 for i, s in enumerate(stat): 1ed

332 coords[0].append(np.full(s['ypix'].shape, i)) 1ed

333 coords[1].append(s['ypix']) 1ed

334 coords[2].append(s['xpix']) 1ed

335 data.append(s['lam']) 1ed

336 pil_coords.append(s['neuropil_mask'] + i * npx) 1ed

337 roi_mask_sp = sparse.COO(list(map(np.concatenate, coords)), np.concatenate(data), shape=shape) 1ed

338 pil_mask_sp = sparse.COO(np.unravel_index(np.concatenate(pil_coords), shape), True, shape=shape) 1ed

339 return sparse.GCXS.from_coo(roi_mask_sp), sparse.GCXS.from_coo(pil_mask_sp) 1ed

340 

341 def _rename_outputs(self, suite2p_dir, frameQC_names, frameQC, rename_dict=None): 

342 """ 

343 Convert suite2p output files to ALF datasets. 

344 

345 This also moves any data.bin and ops.npy files to raw_bin_files for quicker re-runs. 

346 

347 Parameters 

348 ---------- 

349 suite2p_dir : pathlib.Path 

350 The location of the suite2p output (typically session_path/suite2p). 

351 rename_dict : dict or None 

352 The suite2p output filenames and the corresponding ALF name. NB: These files are saved 

353 after transposition. Default is None, i.e. using the default mapping hardcoded in the function below. 

354 

355 Returns 

356 ------- 

357 list of pathlib.Path 

358 All paths found in FOV folders. 

359 """ 

360 if rename_dict is None: 1ed

361 rename_dict = { 1ed

362 'F.npy': 'mpci.ROIActivityF.npy', 

363 'spks.npy': 'mpci.ROIActivityDeconvolved.npy', 

364 'Fneu.npy': 'mpci.ROINeuropilActivityF.npy' 

365 } 

366 fov_dsets = [d[0] for d in self.signature['output_files'] if d[1].startswith('alf/FOV')] 1ed

367 for plane_dir in self._get_plane_paths(suite2p_dir): 1ed

368 # Rename the registered bin file 

369 if (bin_file := plane_dir.joinpath('data.bin')).exists(): 1ed

370 bin_file.rename(plane_dir.joinpath('imaging.frames_motionRegistered.bin')) 1d

371 # Archive the raw suite2p output before renaming 

372 n = int(plane_dir.name.split('plane')[1]) 1ed

373 fov_dir = self.session_path.joinpath('alf', f'FOV_{n:02}') 1ed

374 if fov_dir.exists(): 1ed

375 for f in filter(Path.exists, map(fov_dir.joinpath, fov_dsets)): 1d

376 _logger.debug('Removing old file %s', f.relative_to(self.session_path)) 1d

377 f.unlink() 1d

378 else: 

379 fov_dir.mkdir(parents=True) 1e

380 with ZipFile(fov_dir / '_suite2p_ROIData.raw.zip', 'w', compression=ZIP_DEFLATED) as zf: 1ed

381 for file in plane_dir.iterdir(): 1ed

382 if file.suffix != '.bin': 1ed

383 zf.write(file, arcname=file.name) 1ed

384 # save frameQC in each dir (for now, maybe there will be fov specific frame QC eventually) 

385 if frameQC is not None and len(frameQC) > 0: 1ed

386 np.save(fov_dir.joinpath('mpci.mpciFrameQC.npy'), frameQC) 1d

387 frameQC_names.to_csv(fov_dir.joinpath('mpciFrameQC.names.tsv'), sep='\t', index=False) 1d

388 # extract some other data from suite2p outputs 

389 ops = np.load(plane_dir.joinpath('ops.npy'), allow_pickle=True).item() 1ed

390 stat = np.load(plane_dir.joinpath('stat.npy'), allow_pickle=True) 1ed

391 iscell = np.load(plane_dir.joinpath('iscell.npy')) 1ed

392 # Save suite2p ROI activity outputs in transposed from (n_frames, n_ROI) 

393 for k, v in rename_dict.items(): 1ed

394 np.save(fov_dir.joinpath(v), np.load(plane_dir.joinpath(k)).T) 1ed

395 np.save(fov_dir.joinpath('mpci.badFrames.npy'), np.asarray(ops['badframes'], dtype=bool)) 1ed

396 np.save(fov_dir.joinpath('mpciMeanImage.images.npy'), np.asarray(ops['meanImg'], dtype=float)) 1ed

397 np.save(fov_dir.joinpath('mpciROIs.stackPos.npy'), np.asarray([(*s['med'], 0) for s in stat], dtype=int)) 1ed

398 np.save(fov_dir.joinpath('mpciROIs.cellClassifier.npy'), np.asarray(iscell[:, 1], dtype=float)) 1ed

399 np.save(fov_dir.joinpath('mpciROIs.mpciROITypes.npy'), np.asarray(iscell[:, 0], dtype=np.int16)) 1ed

400 # clusters uuids 

401 uuid_list = ['uuids'] + list(map(str, [uuid.uuid4() for _ in range(len(iscell))])) 1ed

402 with open(fov_dir.joinpath('mpciROIs.uuids.csv'), 'w+') as fid: 1ed

403 fid.write('\n'.join(uuid_list)) 1ed

404 (pd.DataFrame([(0, 'no cell'), (1, 'cell')], columns=['roi_values', 'roi_labels']) 1ed

405 .to_csv(fov_dir.joinpath('mpciROITypes.names.tsv'), sep='\t', index=False)) 

406 # ROI and neuropil masks 

407 roi_mask, pil_mask = self._masks2sparse(stat, ops) 1ed

408 with open(fov_dir.joinpath('mpciROIs.masks.sparse_npz'), 'wb') as fp: 1ed

409 sparse.save_npz(fp, roi_mask) 1ed

410 with open(fov_dir.joinpath('mpciROIs.neuropilMasks.sparse_npz'), 'wb') as fp: 1ed

411 sparse.save_npz(fp, pil_mask) 1ed

412 # Remove the suite2p output files, leaving only the bins and ops 

413 for path in sorted(plane_dir.iterdir()): 1ed

414 if path.name != 'ops.npy' and path.suffix != '.bin': 1ed

415 path.unlink() if path.is_file() else path.rmdir() 1ed

416 # Collect all files in those directories 

417 datasets = self.session_path.joinpath('alf').rglob('FOV_??/*.*.*') 1ed

418 return sorted(x for x in datasets if x.name in fov_dsets) 1ed

419 

420 def load_meta_files(self): 

421 """Load the extracted imaging metadata files. 

422 

423 Loads and consolidates the imaging data metadata from rawImagingData.meta.json files. 

424 These files contain ScanImage metadata extracted from the raw tiff headers by the 

425 function `mesoscopeMetadataExtraction.m` in iblscripts/deploy/mesoscope. 

426 

427 Returns 

428 ------- 

429 dict 

430 Single, consolidated dictionary containing metadata. 

431 list of dict 

432 The meta data for each individual imaging bout. 

433 """ 

434 # Load metadata and make sure all metadata is consistent across FOVs 

435 meta_files = sorted(self.session_path.glob(f'{self.device_collection}/*rawImagingData.meta.*')) 1ka

436 collections = sorted(set(f.parts[-2] for f in meta_files)) 1ka

437 # Check there is exactly 1 meta file per collection 

438 assert len(meta_files) == len(list(self.session_path.glob(self.device_collection))) == len(collections) 1ka

439 raw_meta = map(alfio.load_file_content, meta_files) 1ka

440 all_meta = list(map(mesoscope.patch_imaging_meta, raw_meta)) 1ka

441 return self._consolidate_metadata(all_meta) if len(all_meta) > 1 else all_meta[0], all_meta 1ka

442 

443 @staticmethod 

444 def _consolidate_metadata(meta_data_all: list) -> dict: 

445 """ 

446 Check that the metadata is consistent across all raw imaging folders. 

447 

448 Parameters 

449 ---------- 

450 meta_data_all: list of dicts 

451 List of metadata dictionaries to be checked for consistency. 

452 

453 Returns 

454 ------- 

455 dict 

456 Single, consolidated dictionary containing metadata. 

457 """ 

458 # Ignore the things we don't expect to match 

459 ignore = ('acquisitionStartTime', 'nFrames') 1a

460 ignore_sub = {'rawScanImageMeta': ('ImageDescription', 'Software')} 1a

461 

462 def equal_dicts(a, b, skip=None): 1a

463 ka = set(a).difference(skip or ()) 1a

464 kb = set(b).difference(skip or ()) 1a

465 return ka == kb and all(a[key] == b[key] for key in ka) 1a

466 

467 # Compare each dict with the first one in the list 

468 for meta in meta_data_all[1:]: 1a

469 if meta != meta_data_all[0]: # compare entire object first 1a

470 for k, v in meta_data_all[0].items(): # check key by key 1a

471 if not (equal_dicts(v, meta[k], ignore_sub[k]) # compare sub-dicts... 1a

472 if k in ignore_sub else # ... if we have keys to ignore in test 

473 not (k in ignore or v == meta[k])): 

474 _logger.warning(f'Mismatch in meta data between raw_imaging_data folders for key {k}. ' 1a

475 f'Using meta_data from first folder!') 

476 return meta_data_all[0] 1a

477 

478 @staticmethod 

479 def _consolidate_exptQC(exptQC): 

480 """ 

481 Consolidate exptQC.mat files into a single file. 

482 

483 Parameters 

484 ---------- 

485 exptQC : list of pandas.DataFrame 

486 The loaded 'exptQC.mat' files as squeezed and simplified data frames, with columns 

487 {'frameQC_frames', 'frameQC_names'}. 

488 

489 Returns 

490 ------- 

491 numpy.array 

492 An array of uint8 where 0 indicates good frames, and other values correspond to 

493 experimenter-defined issues (in 'qc_values' column of output data frame). 

494 pandas.DataFrame 

495 A data frame with columns {'qc_values', 'qc_labels'}, the former an unsigned int 

496 corresponding to a QC code; the latter a human-readable QC explanation. 

497 numpy.array 

498 An array of frame indices where QC code != 0. 

499 """ 

500 # Create a new enumeration combining all unique QC labels. 

501 # 'ok' will always have an enum of 0, the rest are determined by order alone 

502 qc_labels = ['ok'] 1ja

503 frame_qc = [] 1ja

504 for e in exptQC: 1ja

505 assert e.keys() >= set(['frameQC_names', 'frameQC_frames']) 1ja

506 # Initialize an NaN array the same size of frameQC_frames to fill with new enum values 

507 frames = np.full(e['frameQC_frames'].shape, fill_value=np.nan) 1ja

508 # May be numpy array of str or a single str, in both cases we cast to list of str 

509 names = list(ensure_list(e['frameQC_names'])) 1ja

510 # For each label for the old enum, populate initialized array with the new one 

511 for i_old, name in enumerate(names): 1ja

512 name = name if len(name) else 'unknown' # handle empty array and empty str 1ja

513 try: 1ja

514 i_new = qc_labels.index(name) 1ja

515 except ValueError: 1ja

516 i_new = len(qc_labels) 1ja

517 qc_labels.append(name) 1ja

518 frames[e['frameQC_frames'] == i_old] = i_new 1ja

519 frame_qc.append(frames) 1ja

520 # Concatenate frames 

521 frame_qc = np.concatenate(frame_qc) 1ja

522 # If any NaNs left over, assign 'unknown' label 

523 if (missing_name := np.isnan(frame_qc)).any(): 1ja

524 try: 1j

525 i = qc_labels.index('unknown') 1j

526 except ValueError: 

527 i = len(qc_labels) 

528 qc_labels.append('unknown') 

529 frame_qc[missing_name] = i 1j

530 frame_qc = frame_qc.astype(np.uint32) # case to uint 1ja

531 bad_frames, = np.where(frame_qc != 0) 1ja

532 # Convert labels to value -> label data frame 

533 frame_qc_names = pd.DataFrame(list(enumerate(qc_labels)), columns=['qc_values', 'qc_labels']) 1ja

534 return frame_qc, frame_qc_names, bad_frames 1ja

535 

536 def get_default_tau(self): 

537 """ 

538 Determine the tau (fluorescence decay) from the subject's genotype. 

539 

540 Returns 

541 ------- 

542 float 

543 The tau value to use. 

544 

545 See Also 

546 -------- 

547 https://suite2p.readthedocs.io/en/latest/settings.html 

548 """ 

549 # These settings are from the suite2P documentation 

550 TAU_MAP = {'G6s': 1.5, 'G6m': 1., 'G6f': .7, 'default': 1.5} 1s

551 _, subject, *_ = session_path_parts(self.session_path) 1s

552 genotype = self.one.alyx.rest('subjects', 'read', id=subject)['genotype'] 1s

553 match = next(filter(None, (re.match(r'.+-(G\d[fms])$', g['allele']) for g in genotype)), None) 1s

554 key = match.groups()[0] if match else 'default' 1s

555 return TAU_MAP.get(key, TAU_MAP['default']) 1s

556 

557 def _meta2ops(self, meta): 

558 """ 

559 Create the ops dictionary for suite2p based on metadata. 

560 

561 Parameters 

562 ---------- 

563 meta: dict 

564 Imaging metadata. 

565 

566 Returns 

567 ------- 

568 dict 

569 Inputs to suite2p run that deviate from default parameters. 

570 """ 

571 # Computing dx and dy 

572 cXY = np.array([fov['Deg']['topLeft'] for fov in meta['FOV']]) 1ka

573 cXY -= np.min(cXY, axis=0) 1ka

574 nXnYnZ = np.array([fov['nXnYnZ'] for fov in meta['FOV']]) 1ka

575 

576 # Currently supporting z-stacks but not supporting dual plane / volumetric imaging, assert that this is not the case 

577 if np.any(nXnYnZ[:, 2] > 1): 1ka

578 raise NotImplementedError('Dual-plane imaging not yet supported, data seems to more than one plane per FOV') 

579 

580 sW = np.sqrt(np.sum((np.array([fov['Deg']['topRight'] for fov in meta['FOV']]) - np.array( 1ka

581 [fov['Deg']['topLeft'] for fov in meta['FOV']])) ** 2, axis=1)) 

582 sH = np.sqrt(np.sum((np.array([fov['Deg']['bottomLeft'] for fov in meta['FOV']]) - np.array( 1ka

583 [fov['Deg']['topLeft'] for fov in meta['FOV']])) ** 2, axis=1)) 

584 pixSizeX = nXnYnZ[:, 0] / sW 1ka

585 pixSizeY = nXnYnZ[:, 1] / sH 1ka

586 dx = np.round(cXY[:, 0] * pixSizeX).astype(dtype=np.int32) 1ka

587 dy = np.round(cXY[:, 1] * pixSizeY).astype(dtype=np.int32) 1ka

588 nchannels = len(meta['channelSaved']) if isinstance(meta['channelSaved'], list) else 1 1ka

589 

590 # Computing number of unique z-planes (slices in tiff) 

591 # FIXME this should work if all FOVs are discrete or if all FOVs are continuous, but may not work for combination of both 

592 slice_ids = [fov['slice_id'] for fov in meta['FOV']] 1ka

593 nplanes = len(set(slice_ids)) 1ka

594 

595 # Figuring out how many SI Rois we have (one unique ROI may have several FOVs) 

596 # FIXME currently unused 

597 # roiUUIDs = np.array([fov['roiUUID'] for fov in meta['FOV']]) 

598 # nrois = len(np.unique(roiUUIDs)) 

599 

600 db = { 1ka

601 'data_path': sorted(map(str, self.session_path.glob(f'{self.device_collection}'))), 

602 'save_path0': str(self.session_path), 

603 'look_one_level_down': False, # don't look in the children folders as that is where the reference data is 

604 'num_workers': self.cpu, # this selects number of cores to parallelize over for the registration step 

605 'num_workers_roi': -1, # for parallelization over FOVs during cell detection, for now don't 

606 'keep_movie_raw': False, 

607 'delete_bin': False, # TODO: delete this on the long run 

608 'batch_size': 500, # SP reduced this from 1000 

609 'nimg_init': 400, 

610 'combined': False, 

611 'nonrigid': True, 

612 'maxregshift': 0.05, # default = 1 

613 'denoise': 1, # whether binned movie should be denoised before cell detection 

614 'block_size': [128, 128], 

615 'save_mat': True, # save the data to Fall.mat 

616 'move_bin': True, # move the binary file to save_path 

617 'mesoscan': True, 

618 'nplanes': nplanes, 

619 'nrois': len(meta['FOV']), 

620 'nchannels': nchannels, 

621 'fs': meta['scanImageParams']['hRoiManager']['scanVolumeRate'], 

622 'lines': [list(np.asarray(fov['lineIdx']) - 1) for fov in meta['FOV']], # subtracting 1 to make 0-based 

623 'slices': slice_ids, # this tells us which FOV corresponds to which tiff slices 

624 'tau': self.get_default_tau(), # deduce the GCamp used from Alyx mouse line (defaults to 1.5; that of GCaMP6s) 

625 'functional_chan': 1, # for now, eventually find(ismember(meta.channelSaved == meta.channelID.green)) 

626 'align_by_chan': 1, # for now, eventually find(ismember(meta.channelSaved == meta.channelID.red)) 

627 'dx': dx.tolist(), 

628 'dy': dy.tolist() 

629 } 

630 

631 return db 1ka

632 

633 @staticmethod 

634 def _get_plane_paths(path): 

635 """Return list of sorted suite2p plane folder paths. 

636 

637 Parameters 

638 ---------- 

639 path : pathlib.Path 

640 The path containing plane folders. 

641 

642 Returns 

643 ------- 

644 list of pathlib.Path 

645 The plane folder paths, ordered by number. 

646 """ 

647 pattern = re.compile(r'(?<=^plane)\d+$') 1xaed

648 return sorted(path.glob('plane?*'), key=lambda x: int(pattern.search(x.name).group())) 1xaed

649 

650 def bin_per_plane(self, metadata, **kwargs): 

651 """ 

652 Extracts a binary data file of imaging data per imaging plane. 

653 

654 Parameters 

655 ---------- 

656 metadata : dict 

657 A dictionary of extracted metadata. 

658 save_path0 : str, pathlib.Path 

659 The root path of the suite2p bin output. 

660 save_folder : str 

661 The subfolder within `save_path0` to save the suite2p bin output. 

662 kwargs 

663 Other optional arguments to overwrite the defaults for. 

664 

665 Returns 

666 ------- 

667 list of pathlib.Path 

668 Ordered list of output plane folders containing binary data and ops. 

669 dict 

670 Suite2p's modified options. 

671 """ 

672 import suite2p.io 1a

673 

674 options = ('nplanes', 'data_path', 'save_path0', 'save_folder', 'fast_disk', 'batch_size', 1a

675 'nchannels', 'keep_movie_raw', 'look_one_level_down', 'lines', 'dx', 'dy', 'force_sktiff', 

676 'do_registration', 'slices') 

677 ops = self._meta2ops(metadata) 1a

678 ops['force_sktiff'] = False 1a

679 ops['do_registration'] = True 1a

680 ops = {k: v for k, v in ops.items() if k in options} 1a

681 ops.update(kwargs) 1a

682 ops['save_path0'] = str(ops['save_path0']) # Path objs must be str for suite2p 1a

683 # Update the task kwargs attribute as it will be stored in the arguments json field in alyx 

684 self.kwargs = ops.copy() 1a

685 

686 ret = suite2p.io.mesoscan_to_binary(ops.copy()) 1a

687 

688 # Get ordered list of plane folders 

689 out_path = Path(ret['save_path0'], ret['save_folder']) 1a

690 assert out_path.exists() 1a

691 planes = self._get_plane_paths(out_path) 1a

692 assert len(planes) == ret['nplanes'] 1a

693 

694 return planes, ret 1a

695 

696 def image_motion_registration(self, ops): 

697 """Perform motion registration. 

698 

699 Parameters 

700 ---------- 

701 ops : dict 

702 A dict of suite2p options. 

703 

704 Returns 

705 ------- 

706 dict 

707 A dictionary of registration metrics: "regDX" is nPC x 3, where X[:,0] 

708 is rigid, X[:,1] is average nonrigid, X[:,2] is max nonrigid shifts; 

709 "regPC" is average of top and bottom frames for each PC; "tPC" is PC 

710 across time frames; "reg_metrics_avg" is the average of "regDX"; 

711 "reg_metrics_max" is the maximum of "regDX". 

712 

713 """ 

714 import suite2p 1r

715 ops['do_registration'] = True 1r

716 ops['do_regmetrics'] = True 1r

717 ops['roidetect'] = False 1r

718 ret = suite2p.run_plane(ops) 1r

719 metrics = {k: ret.get(k, None) for k in ('regDX', 'regPC', 'tPC')} 1r

720 has_metrics = ops['do_regmetrics'] and metrics['regDX'] is not None 1r

721 metrics['reg_metrics_avg'] = np.mean(metrics['regDX'], axis=0) if has_metrics else None 1r

722 metrics['reg_metrics_max'] = np.max(metrics['regDX'], axis=0) if has_metrics else None 1r

723 return metrics 1r

724 

725 def roi_detection(self, ops): 

726 """Perform ROI detection. 

727 

728 Parameters 

729 ---------- 

730 ops : dict 

731 A dict of suite2p options. 

732 

733 Returns 

734 ------- 

735 dict 

736 An updated copy of the ops after running ROI detection. 

737 """ 

738 import suite2p 1v

739 ops['do_registration'] = False 1v

740 ops['roidetect'] = True 1v

741 ret = suite2p.run_plane(ops) 1v

742 return ret 1v

743 

744 def _run(self, rename_files=True, use_badframes=True, **kwargs): 

745 """ 

746 Process inputs, run suite2p and make outputs alf compatible. 

747 

748 The suite2p processing takes place in a 'suite2p' folder within the session path. After running, 

749 the data.bin files are moved to 'raw_bin_files' and the rest of the folder is zipped up and moved 

750 to 'alf/ 

751 

752 Parameters 

753 ---------- 

754 rename_files: bool 

755 Whether to rename and reorganize the suite2p outputs to be alf compatible. Defaults is True. 

756 use_badframes: bool 

757 Whether to exclude bad frames indicated by the experimenter in badframes.mat. 

758 overwrite : bool 

759 Whether to re-perform extraction and motion registration. 

760 do_registration : bool 

761 Whether to perform motion registration. 

762 roidetect : bool 

763 Whether to perform ROI detection. 

764 

765 Returns 

766 ------- 

767 list of pathlib.Path 

768 All files created by the task. 

769 

770 """ 

771 

772 """ Metadata and parameters """ 1a

773 overwrite = kwargs.pop('overwrite', self.overwrite) 1a

774 # Load and consolidate the image metadata from JSON files 

775 metadata, all_meta = self.load_meta_files() 1a

776 

777 # Create suite2p output folder in root session path 

778 raw_image_collections = sorted(self.session_path.glob(f'{self.device_collection}')) 1a

779 save_path = self.session_path.joinpath(save_folder := 'suite2p') 1a

780 

781 # Check for previous intermediate files 

782 plane_folders = self._get_plane_paths(save_path) 1a

783 if len(plane_folders) > 0: 1a

784 for bin_file in save_path.glob('plane?*/imaging.frames_motionRegistered.bin'): 1a

785 # If there is a previously registered bin file, rename back to data.bin 

786 # This file name is hard-coded in suite2p but is renamed in _rename_outputs 

787 # If a data.bin file already exists, we don't want to overwrite it - 

788 # it may be from a manual registration or incomplete run 

789 if not bin_file.with_name('data.bin').exists(): 1a

790 bin_file.rename(bin_file.with_stem('data')) 1a

791 

792 if len(plane_folders) == 0 or overwrite: 1a

793 _logger.info('Extracting tif data per plane') 1a

794 # Ingest tiff files 

795 plane_folders, _ = self.bin_per_plane(metadata, save_folder=save_folder, save_path0=self.session_path) 1a

796 

797 """ Bad frames """ 1a

798 # exptQC.mat contains experimenter QC values that may not affect ROI detection (e.g. noises, pauses) 

799 qc_datasets = dataset_from_name('exptQC.mat', self.input_files) 1a

800 qc_paths = [next(self.session_path.glob(d.glob_pattern), None) for d in qc_datasets] 1a

801 qc_paths = sorted(map(str, filter(None, qc_paths))) 1a

802 exptQC = [loadmat(p, squeeze_me=True, simplify_cells=True) for p in qc_paths] 1a

803 if len(exptQC) > 0: 1a

804 frameQC, frameQC_names, _ = self._consolidate_exptQC(exptQC) 1a

805 else: 

806 _logger.warning('No frame QC (exptQC.mat) files found.') 

807 frameQC = np.array([], dtype='u1') 

808 frameQC_names = pd.DataFrame(columns=['qc_values', 'qc_labels']) 

809 

810 # If applicable, save as bad_frames.npy in first raw_imaging_folder for suite2p 

811 # badframes.mat contains QC values that do affect ROI detection (e.g. no PMT, lens artefacts) 

812 badframes = np.array([], dtype='uint32') 1a

813 total_frames = 0 1a

814 # Ensure all indices are relative to total cumulative frames 

815 for m, collection in zip(all_meta, raw_image_collections): 1a

816 badframes_path = self.session_path.joinpath(collection, 'badframes.mat') 1a

817 if badframes_path.exists(): 1a

818 raw_mat = loadmat(badframes_path, squeeze_me=True, simplify_cells=True) 1a

819 badframes = np.r_[badframes, raw_mat['badframes'].astype('uint32') + total_frames] 1a

820 total_frames += m['nFrames'] 1a

821 if len(badframes) > 0 and use_badframes is True: 1a

822 # The badframes array should always be a subset of the frameQC array 

823 assert np.max(badframes) < frameQC.size and np.all(frameQC[badframes]) 1a

824 np.save(raw_image_collections[0].joinpath('bad_frames.npy'), badframes) 1a

825 

826 """ Suite2p """ 1a

827 # Create alf if is doesn't exist 

828 self.session_path.joinpath('alf').mkdir(exist_ok=True) 1a

829 

830 # Perform registration 

831 if kwargs.get('do_registration', True): 1a

832 _logger.info('Performing registration') 1a

833 for plane in plane_folders: 1a

834 ops = np.load(plane.joinpath('ops.npy'), allow_pickle=True).item() 1a

835 ops.update(kwargs) 1a

836 # (ops['do_registration'], ops['reg_file'], ops['meanImg']) 

837 _ = self.image_motion_registration(ops) 1a

838 # TODO Handle metrics and QC here 

839 

840 # ROI detection 

841 if kwargs.get('roidetect', True): 1a

842 _logger.info('Performing ROI detection') 1a

843 for plane in plane_folders: 1a

844 ops = np.load(plane.joinpath('ops.npy'), allow_pickle=True).item() 1a

845 ops.update(kwargs) 1a

846 self.roi_detection(ops) 1a

847 

848 """ Outputs """ 1a

849 # Save and rename other outputs 

850 if rename_files: 1a

851 self._rename_outputs(save_path, frameQC_names, frameQC) 1a

852 # Only return output file that are in the signature (for registration) 

853 out_files = chain.from_iterable(map(lambda x: x.find_files(self.session_path)[1], self.output_files)) 1a

854 else: 

855 out_files = save_path.rglob('*.*') # Output all files in suite2p folder 1a

856 

857 return list(out_files) 1a

858 

859 

860class MesoscopeSync(base_tasks.MesoscopeTask): 

861 """Extract the frame times from the main DAQ.""" 

862 

863 priority = 40 

864 job_size = 'small' 

865 

866 @property 

867 def signature(self): 

868 I = ExpectedDataset.input # noqa 1nl

869 signature = { 1nl

870 'input_files': [I(f'_{self.sync_namespace}_DAQdata.raw.npy', self.sync_collection, True), 

871 I(f'_{self.sync_namespace}_DAQdata.timestamps.npy', self.sync_collection, True), 

872 I(f'_{self.sync_namespace}_DAQdata.meta.json', self.sync_collection, True), 

873 I('_ibl_rawImagingData.meta.json', self.device_collection, True, unique=False), 

874 I('rawImagingData.times_scanImage.npy', self.device_collection, True, True, unique=False), 

875 I(f'_{self.sync_namespace}_softwareEvents.log.htsv', self.sync_collection, False), ], 

876 'output_files': [('mpci.times.npy', 'alf/FOV*', True), 

877 ('mpciStack.timeshift.npy', 'alf/FOV*', True),] 

878 } 

879 return signature 1nl

880 

881 def _run(self, **kwargs): 

882 """ 

883 Extract the imaging times for all FOVs. 

884 

885 Returns 

886 ------- 

887 list of pathlib.Path 

888 Files containing frame timestamps for individual FOVs and time offsets for each line scan. 

889 

890 """ 

891 # TODO function to determine nFOVs 

892 try: 1nl

893 alf_path = self.session_path / self.sync_collection 1nl

894 events = alfio.load_object(alf_path, 'softwareEvents').get('log') 1nl

895 except alferr.ALFObjectNotFound: 1l

896 events = None 1l

897 if events is None or events.empty: 1nl

898 _logger.debug('No software events found for session %s', self.session_path) 1l

899 all_collections = flatten(map(lambda x: x.identifiers, self.input_files))[::3] 1nl

900 collections = set(filter(lambda x: fnmatch(x, self.device_collection), all_collections)) 1nl

901 # Load first meta data file to determine the number of FOVs 

902 # Changing FOV between imaging bouts is not supported currently! 

903 self.rawImagingData = alfio.load_object(self.session_path / next(iter(collections)), 'rawImagingData') 1nl

904 self.rawImagingData['meta'] = mesoscope.patch_imaging_meta(self.rawImagingData['meta']) 1nl

905 n_FOVs = len(self.rawImagingData['meta']['FOV']) 1nl

906 sync, chmap = self.load_sync() # Extract sync data from raw DAQ data 1nl

907 legacy = kwargs.get('legacy', False) # this option may be removed in the future once fully tested 1nl

908 mesosync = mesoscope.MesoscopeSyncTimeline(self.session_path, n_FOVs) 1nl

909 _, out_files = mesosync.extract( 1nl

910 save=True, sync=sync, chmap=chmap, device_collection=collections, events=events, use_volume_counter=legacy) 

911 return out_files 1nl

912 

913 

914class MesoscopeFOV(base_tasks.MesoscopeTask): 

915 """Create FOV and FOV location objects in Alyx from metadata.""" 

916 

917 priority = 40 

918 job_size = 'small' 

919 

920 @property 

921 def signature(self): 

922 signature = { 1g

923 'input_files': [('_ibl_rawImagingData.meta.json', self.device_collection, True), 

924 ('mpciROIs.stackPos.npy', 'alf/FOV*', True)], 

925 'output_files': [('mpciMeanImage.brainLocationIds*.npy', 'alf/FOV_*', True), 

926 ('mpciMeanImage.mlapdv*.npy', 'alf/FOV_*', True), 

927 ('mpciROIs.mlapdv*.npy', 'alf/FOV_*', True), 

928 ('mpciROIs.brainLocationIds*.npy', 'alf/FOV_*', True), 

929 ('_ibl_rawImagingData.meta.json', self.device_collection, True)] 

930 } 

931 return signature 1g

932 

933 def _run(self, *args, provenance=Provenance.ESTIMATE, **kwargs): 

934 """ 

935 Register fields of view (FOV) to Alyx and extract the coordinates and IDs of each ROI. 

936 

937 Steps: 

938 1. Save the mpciMeanImage.brainLocationIds_estimate and mlapdv datasets. 

939 2. Use mean image coordinates and ROI stack position datasets to extract brain location 

940 of each ROI. 

941 3. Register the location of each FOV in Alyx. 

942 

943 Parameters 

944 ---------- 

945 provenance : Provenance 

946 The provenance of the coordinates in the meta file. For all but 'HISTOLOGY', the 

947 provenance is added as a dataset suffix. Defaults to ESTIMATE. 

948 

949 Returns 

950 ------- 

951 dict 

952 The newly created FOV Alyx record. 

953 list 

954 The newly created FOV location Alyx records. 

955 

956 Notes 

957 ----- 

958 - Once the FOVs have been registered they cannot be updated with this task. Rerunning this 

959 task will result in an error. 

960 - This task modifies the first meta JSON file. All meta files are registered by this task. 

961 """ 

962 # Load necessary data 

963 (filename, collection, _), *_ = self.signature['input_files'] 1g

964 meta_files = sorted(self.session_path.glob(f'{collection}/{filename}')) 1g

965 meta = mesoscope.patch_imaging_meta(alfio.load_file_content(meta_files[0]) or {}) 1g

966 nFOV = len(meta.get('FOV', [])) 1g

967 

968 suffix = None if provenance is Provenance.HISTOLOGY else provenance.name.lower() 1g

969 _logger.info('Extracting %s MLAPDV datasets', suffix or 'final') 1g

970 

971 # Extract mean image MLAPDV coordinates and brain location IDs 

972 mean_image_mlapdv, mean_image_ids = self.project_mlapdv(meta) 1g

973 

974 # Save the meta data file with new coordinate fields 

975 with open(meta_files[0], 'w') as fp: 1g

976 json.dump(meta, fp) 1g

977 

978 # Save the mean image datasets 

979 mean_image_files = [] 1g

980 assert set(mean_image_mlapdv.keys()) == set(mean_image_ids.keys()) and len(mean_image_ids) == nFOV 1g

981 for i in range(nFOV): 1g

982 alf_path = self.session_path.joinpath('alf', f'FOV_{i:02}') 1g

983 for attr, arr, sfx in (('mlapdv', mean_image_mlapdv[i], suffix), 1g

984 ('brainLocationIds', mean_image_ids[i], ('ccf', '2017', suffix))): 

985 mean_image_files.append(alf_path / to_alf('mpciMeanImage', attr, 'npy', timescale=sfx)) 1g

986 mean_image_files[-1].parent.mkdir(parents=True, exist_ok=True) 1g

987 np.save(mean_image_files[-1], arr) 1g

988 

989 # Extract ROI MLAPDV coordinates and brain location IDs 

990 roi_mlapdv, roi_brain_ids = self.roi_mlapdv(nFOV, suffix=suffix) 1g

991 

992 # Write MLAPDV + brain location ID of ROIs to disk 

993 roi_files = [] 1g

994 assert set(roi_mlapdv.keys()) == set(roi_brain_ids.keys()) and len(roi_mlapdv) == nFOV 1g

995 for i in range(nFOV): 1g

996 alf_path = self.session_path.joinpath('alf', f'FOV_{i:02}') 1g

997 for attr, arr, sfx in (('mlapdv', roi_mlapdv[i], suffix), 1g

998 ('brainLocationIds', roi_brain_ids[i], ('ccf', '2017', suffix))): 

999 roi_files.append(alf_path / to_alf('mpciROIs', attr, 'npy', timescale=sfx)) 1g

1000 np.save(roi_files[-1], arr) 1g

1001 

1002 # Register FOVs in Alyx 

1003 self.register_fov(meta, suffix) 1g

1004 

1005 return sorted([*meta_files, *roi_files, *mean_image_files]) 1g

1006 

1007 def update_surgery_json(self, meta, normal_vector): 

1008 """ 

1009 Update surgery JSON with surface normal vector. 

1010 

1011 Adds the key 'surface_normal_unit_vector' to the most recent surgery JSON, containing the 

1012 provided three element vector. The recorded craniotomy center must match the coordinates 

1013 in the provided meta file. 

1014 

1015 Parameters 

1016 ---------- 

1017 meta : dict 

1018 The imaging meta data file containing the 'centerMM' key. 

1019 normal_vector : array_like 

1020 A three element unit vector normal to the surface of the craniotomy center. 

1021 

1022 Returns 

1023 ------- 

1024 dict 

1025 The updated surgery record, or None if no surgeries found. 

1026 """ 

1027 if not self.one or self.one.offline: 1oc

1028 _logger.warning('failed to update surgery JSON: ONE offline') 1oc

1029 return 1oc

1030 # Update subject JSON with unit normal vector of craniotomy centre (used in histology) 

1031 subject = self.one.path2ref(self.session_path, parse=False)['subject'] 1o

1032 surgeries = self.one.alyx.rest('surgeries', 'list', subject=subject, procedure='craniotomy') 1o

1033 if not surgeries: 1o

1034 _logger.error(f'Surgery not found for subject "{subject}"') 1o

1035 return 1o

1036 surgery = surgeries[0] # Check most recent surgery in list 1o

1037 center = (meta['centerMM']['ML'], meta['centerMM']['AP']) 1o

1038 match = (k for k, v in (surgery['json'] or {}).items() if 1o

1039 str(k).startswith('craniotomy') and np.allclose(v['center'], center)) 

1040 if (key := next(match, None)) is None: 1o

1041 _logger.error('Failed to update surgery JSON: no matching craniotomy found') 1o

1042 return surgery 1o

1043 data = {key: {**surgery['json'][key], 'surface_normal_unit_vector': tuple(normal_vector)}} 1o

1044 surgery['json'] = self.one.alyx.json_field_update('subjects', subject, data=data) 1o

1045 return surgery 1o

1046 

1047 def roi_mlapdv(self, nFOV: int, suffix=None): 

1048 """ 

1049 Extract ROI MLAPDV coordinates and brain location IDs. 

1050 

1051 MLAPDV coordinates are in um relative to bregma. Location IDs are from the 2017 Allen 

1052 common coordinate framework atlas. 

1053 

1054 Parameters 

1055 ---------- 

1056 nFOV : int 

1057 The number of fields of view acquired. 

1058 suffix : {None, 'estimate'} 

1059 The attribute suffix of the mpciMeanImage datasets to load. If generating from 

1060 estimates, the suffix should be 'estimate'. 

1061 

1062 Returns 

1063 ------- 

1064 dict of int : numpy.array 

1065 A map of field of view to ROI MLAPDV coordinates. 

1066 dict of int : numpy.array 

1067 A map of field of view to ROI brain location IDs. 

1068 """ 

1069 all_mlapdv = {} 1g

1070 all_brain_ids = {} 1g

1071 for n in range(nFOV): 1g

1072 alf_path = self.session_path.joinpath('alf', f'FOV_{n:02}') 1g

1073 

1074 # Load neuron centroids in pixel space 

1075 stack_pos_file = next(alf_path.glob('mpciROIs.stackPos*'), None) 1g

1076 if not stack_pos_file: 1g

1077 raise FileNotFoundError(alf_path / 'mpci.stackPos*') 1g

1078 stack_pos = alfio.load_file_content(stack_pos_file) 1g

1079 

1080 # Load MLAPDV + brain location ID maps of pixels 

1081 mpciMeanImage = alfio.load_object( 1g

1082 alf_path, 'mpciMeanImage', attribute=['mlapdv', 'brainLocationIds']) 

1083 

1084 # Get centroid MLAPDV + brainID by indexing pixel-map with centroid locations 

1085 mlapdv = np.full(stack_pos.shape, np.nan) 1g

1086 brain_ids = np.full(stack_pos.shape[0], np.nan) 1g

1087 for i in np.arange(stack_pos.shape[0]): 1g

1088 idx = (stack_pos[i, 0], stack_pos[i, 1]) 1g

1089 sfx = f'_{suffix}' if suffix else '' 1g

1090 mlapdv[i, :] = mpciMeanImage['mlapdv' + sfx][idx] 1g

1091 brain_ids[i] = mpciMeanImage['brainLocationIds_ccf_2017' + sfx][idx] 1g

1092 assert ~np.isnan(brain_ids).any() 1g

1093 all_brain_ids[n] = brain_ids.astype(int) 1g

1094 all_mlapdv[n] = mlapdv 1g

1095 

1096 return all_mlapdv, all_brain_ids 1g

1097 

1098 @staticmethod 

1099 def get_provenance(filename): 

1100 """ 

1101 Get the field of view provenance from a mpciMeanImage or mpciROIs dataset. 

1102 

1103 Parameters 

1104 ---------- 

1105 filename : str, pathlib.Path 

1106 A filename to get the provenance from. 

1107 

1108 Returns 

1109 ------- 

1110 Provenance 

1111 The provenance of the file. 

1112 """ 

1113 filename = Path(filename).name 1w

1114 timescale = (filename_parts(filename)[3] or '').split('_') 1w

1115 provenances = [i.name.lower() for i in Provenance] 1w

1116 provenance = (Provenance[x.upper()] for x in timescale if x in provenances) 1w

1117 return next(provenance, None) or Provenance.HISTOLOGY 1w

1118 

1119 def register_fov(self, meta: dict, suffix: str = None) -> (list, list): 

1120 """ 

1121 Create FOV on Alyx. 

1122 

1123 Assumes field of view recorded perpendicular to objective. 

1124 Assumes field of view is plane (negligible volume). 

1125 

1126 Required Alyx fixtures: 

1127 - experiments.ImagingType(name='mesoscope') 

1128 - experiments.CoordinateSystem(name='IBL-Allen') 

1129 

1130 Parameters 

1131 ---------- 

1132 meta : dict 

1133 The raw imaging meta data from _ibl_rawImagingData.meta.json. 

1134 suffix : str 

1135 The file attribute suffixes to load from the mpciMeanImage object. Either 'estimate' or 

1136 None. No suffix means the FOV location provenance will be L (Landmark). 

1137 

1138 Returns 

1139 ------- 

1140 list of dict 

1141 A list registered of field of view entries from Alyx. 

1142 

1143 TODO Determine dual plane ID for JSON field 

1144 """ 

1145 dry = self.one is None or self.one.offline 1i

1146 alyx_fovs = [] 1i

1147 # Count the number of slices per stack ID: only register stacks that contain more than one slice. 

1148 slice_counts = Counter(f['roiUUID'] for f in meta.get('FOV', [])) 1i

1149 # Create a new stack in Alyx for all stacks containing more than one slice. 

1150 # Map of ScanImage ROI UUID to Alyx ImageStack UUID. 

1151 if dry: 1i

1152 stack_ids = {i: uuid.uuid4() for i in slice_counts if slice_counts[i] > 1} 1i

1153 else: 

1154 stack_ids = {i: self.one.alyx.rest('imaging-stack', 'create', data={'name': i})['id'] 1i

1155 for i in slice_counts if slice_counts[i] > 1} 

1156 

1157 for i, fov in enumerate(meta.get('FOV', [])): 1i

1158 assert set(fov.keys()) >= {'MLAPDV', 'nXnYnZ', 'roiUUID'} 1i

1159 # Field of view 

1160 alyx_FOV = { 1i

1161 'session': self.session_path.as_posix() if dry else str(self.path2eid()), 

1162 'imaging_type': 'mesoscope', 'name': f'FOV_{i:02}', 

1163 'stack': stack_ids.get(fov['roiUUID']) 

1164 } 

1165 if dry: 1i

1166 print(alyx_FOV) 1i

1167 alyx_FOV['location'] = [] 1i

1168 alyx_fovs.append(alyx_FOV) 1i

1169 else: 

1170 alyx_fovs.append(self.one.alyx.rest('fields-of-view', 'create', data=alyx_FOV)) 1i

1171 

1172 # Field of view location 

1173 data = {'field_of_view': alyx_fovs[-1].get('id'), 1i

1174 'default_provenance': True, 

1175 'coordinate_system': 'IBL-Allen', 

1176 'n_xyz': fov['nXnYnZ']} 

1177 if suffix: 1i

1178 data['provenance'] = suffix[0].upper() 1i

1179 

1180 # Convert coordinates to 4 x 3 array (n corners by n dimensions) 

1181 # x1 = top left ml, y1 = top left ap, y2 = top right ap, etc. 

1182 coords = [fov['MLAPDV'][key] for key in ('topLeft', 'topRight', 'bottomLeft', 'bottomRight')] 1i

1183 coords = np.vstack(coords).T 1i

1184 data.update({k: arr.tolist() for k, arr in zip('xyz', coords)}) 1i

1185 

1186 # Load MLAPDV + brain location ID maps of pixels 

1187 filename = 'mpciMeanImage.brainLocationIds_ccf_2017' + (f'_{suffix}' if suffix else '') + '.npy' 1i

1188 filepath = self.session_path.joinpath('alf', f'FOV_{i:02}', filename) 1i

1189 mean_image_ids = alfio.load_file_content(filepath) 1i

1190 

1191 data['brain_region'] = np.unique(mean_image_ids).astype(int).tolist() 1i

1192 

1193 if dry: 1i

1194 print(data) 1i

1195 alyx_FOV['location'].append(data) 1i

1196 else: 

1197 alyx_fovs[-1]['location'].append(self.one.alyx.rest('fov-location', 'create', data=data)) 1i

1198 return alyx_fovs 1i

1199 

1200 def load_triangulation(self): 

1201 """ 

1202 Load the surface triangulation file. 

1203 

1204 A triangle mesh of the smoothed convex hull of the dorsal surface of the mouse brain, 

1205 generated from the 2017 Allen 10um annotation volume. This triangulation was generated in 

1206 MATLAB. 

1207 

1208 Returns 

1209 ------- 

1210 points : numpy.array 

1211 An N by 3 float array of x-y vertices, defining all points of the triangle mesh. These 

1212 are in um relative to the IBL bregma coordinates. 

1213 connectivity_list : numpy.array 

1214 An N by 3 integer array of vertex indices defining all points that form a triangle. 

1215 """ 

1216 fixture_path = Path(mesoscope.__file__).parent.joinpath('mesoscope') 1c

1217 surface_triangulation = np.load(fixture_path / 'surface_triangulation.npz') 1c

1218 points = surface_triangulation['points'].astype('f8') 1c

1219 connectivity_list = surface_triangulation['connectivity_list'] 1c

1220 surface_triangulation.close() 1c

1221 return points, connectivity_list 1c

1222 

1223 def project_mlapdv(self, meta, atlas=None): 

1224 """ 

1225 Calculate the mean image pixel locations in MLAPDV coordinates and determine the brain 

1226 location IDs. 

1227 

1228 MLAPDV coordinates are in um relative to bregma. Location IDs are from the 2017 Allen 

1229 common coordinate framework atlas. 

1230 

1231 Parameters 

1232 ---------- 

1233 meta : dict 

1234 The raw imaging data meta file, containing coordinates for the centre of each field of 

1235 view. 

1236 atlas : ibllib.atlas.Atlas 

1237 An atlas instance. 

1238 

1239 Returns 

1240 ------- 

1241 dict 

1242 A map of FOV number (int) to mean image MLAPDV coordinates as a 2D numpy array. 

1243 dict 

1244 A map of FOV number (int) to mean image brain location IDs as a 2D numpy int array. 

1245 """ 

1246 mlapdv = {} 1c

1247 location_id = {} 1c

1248 # Use the MRI atlas as this applies scaling, particularly along the DV axis to (hopefully) 

1249 # more accurately represent the living brain. 

1250 atlas = atlas or MRITorontoAtlas(res_um=10) 1c

1251 # The centre of the craniotomy / imaging window 

1252 coord_ml = meta['centerMM']['ML'] * 1e3 # mm -> um 1c

1253 coord_ap = meta['centerMM']['AP'] * 1e3 # mm -> um 1c

1254 pt = np.array([coord_ml, coord_ap]) 1c

1255 

1256 points, connectivity_list = self.load_triangulation() 1c

1257 # Only keep faces that have normals pointing up (positive DV value). 

1258 # Calculate the normal vector pointing out of the convex hull. 

1259 triangles = points[connectivity_list, :] 1c

1260 normals = surface_normal(triangles) 1c

1261 up_faces, = np.where(normals[:, -1] > 0) 1c

1262 # only keep triangles that have normal vector with positive DV component 

1263 dorsal_connectivity_list = connectivity_list[up_faces, :] 1c

1264 # Flatten triangulation by dropping the dorsal coordinates and find the location of the 

1265 # window center (we convert mm -> um here) 

1266 face_ind = find_triangle(pt * 1e-3, points[:, :2] * 1e-3, dorsal_connectivity_list.astype(np.intp)) 1c

1267 assert face_ind != -1 1c

1268 

1269 dorsal_triangle = points[dorsal_connectivity_list[face_ind, :], :] 1c

1270 

1271 # Get the surface normal unit vector of dorsal triangle 

1272 normal_vector = surface_normal(dorsal_triangle) 1c

1273 

1274 # Update the surgery JSON field with normal unit vector, for use in histology alignment 

1275 self.update_surgery_json(meta, normal_vector) 1c

1276 

1277 # find the coordDV that sits on the triangular face and had [coordML, coordAP] coordinates; 

1278 # the three vertices defining the triangle 

1279 face_vertices = points[dorsal_connectivity_list[face_ind, :], :] 1c

1280 

1281 # all the vertices should be on the plane ax + by + cz = 1, so we can find 

1282 # the abc coefficients by inverting the three equations for the three vertices 

1283 abc, *_ = np.linalg.lstsq(face_vertices, np.ones(3), rcond=None) 1c

1284 

1285 # and then find a point on that plane that corresponds to a given x-y 

1286 # coordinate (which is ML-AP coordinate) 

1287 coord_dv = (1 - pt @ abc[:2]) / abc[2] 1c

1288 

1289 # We should not use the actual surface of the brain for this, as it might be in one of the sulci 

1290 # DO NOT USE THIS: 

1291 # coordDV = interp2(axisMLmm, axisAPmm, surfaceDV, coordML, coordAP) 

1292 

1293 # Now we need to span the plane of the coverslip with two orthogonal unit vectors. 

1294 # We start with vY, because the order is important and we usually have less 

1295 # tilt along AP (pitch), which will cause less deviation in vX from pure ML. 

1296 vY = np.array([0, normal_vector[2], -normal_vector[1]]) # orthogonal to the normal of the plane 1c

1297 vX = np.cross(vY, normal_vector) # orthogonal to n and to vY 1c

1298 # normalize and flip the sign if necessary 

1299 vX = vX / np.sqrt(vX @ vX) * np.sign(vX[0]) # np.sqrt(vY @ vY) == LR norm of vX 1c

1300 vY = vY / np.sqrt(vY @ vY) * np.sign(vY[1]) 1c

1301 

1302 # what are the dimensions of the data arrays (ap, ml, dv) 

1303 (nAP, nML, nDV) = atlas.image.shape 1c

1304 # Let's shift the coordinates relative to bregma 

1305 voxel_size = atlas.res_um # [um] resolution of the atlas 1c

1306 bregma_coords = ALLEN_CCF_LANDMARKS_MLAPDV_UM['bregma'] / voxel_size # (ml, ap, dv) 1c

1307 axis_ml_um = (np.arange(nML) - bregma_coords[0]) * voxel_size 1c

1308 axis_ap_um = (np.arange(nAP) - bregma_coords[1]) * voxel_size * -1. 1c

1309 axis_dv_um = (np.arange(nDV) - bregma_coords[2]) * voxel_size * -1. 1c

1310 

1311 # projection of FOVs on the brain surface to get ML-AP-DV coordinates 

1312 _logger.info('Projecting in 3D') 1c

1313 for i, fov in enumerate(meta['FOV']): # i, fov = next(enumerate(meta['FOV'])) 1c

1314 start_time = time.time() 1c

1315 _logger.info(f'FOV {i + 1}/{len(meta["FOV"])}') 1c

1316 y_px_idx, x_px_idx = np.mgrid[0:fov['nXnYnZ'][0], 0:fov['nXnYnZ'][1]] 1c

1317 

1318 # xx and yy are in mm in coverslip space 

1319 points = ((0, fov['nXnYnZ'][0] - 1), (0, fov['nXnYnZ'][1] - 1)) 1c

1320 # The four corners of the FOV, determined by taking the center of the craniotomy in MM, 

1321 # the x-y coordinates of the imaging window center (from the tiled reference image) in 

1322 # galvanometer units, and the x-y coordinates of the FOV center in galvanometer units. 

1323 values = [[fov['MM']['topLeft'][0], fov['MM']['topRight'][0]], 1c

1324 [fov['MM']['bottomLeft'][0], fov['MM']['bottomRight'][0]]] 

1325 values = np.array(values) * 1e3 # mm -> um 1c

1326 xx = interpn(points, values, (y_px_idx, x_px_idx)) 1c

1327 

1328 values = [[fov['MM']['topLeft'][1], fov['MM']['topRight'][1]], 1c

1329 [fov['MM']['bottomLeft'][1], fov['MM']['bottomRight'][1]]] 

1330 values = np.array(values) * 1e3 # mm -> um 1c

1331 yy = interpn(points, values, (y_px_idx, x_px_idx)) 1c

1332 

1333 xx = xx.flatten() - coord_ml 1c

1334 yy = yy.flatten() - coord_ap 1c

1335 

1336 # rotate xx and yy in 3D 

1337 # the coords are still on the coverslip, but now have 3D values 

1338 coords = np.outer(xx, vX) + np.outer(yy, vY) # (vX * xx) + (vY * yy) 1c

1339 coords = coords + [coord_ml, coord_ap, coord_dv] 1c

1340 

1341 # for each point of the FOV create a line parametrization (trajectory normal to the coverslip plane). 

1342 # start just above the coverslip and go 3 mm down, should be enough to 'meet' the brain 

1343 t = np.arange(-voxel_size, 3e3, voxel_size) 1c

1344 

1345 # Find the MLAPDV atlas coordinate and brain location of each pixel. 

1346 MLAPDV, annotation = _update_points( 1c

1347 t, normal_vector, coords, axis_ml_um, axis_ap_um, axis_dv_um, atlas.label) 

1348 annotation = atlas.regions.index2id(annotation) # convert annotation indices to IDs 1c

1349 

1350 if np.any(np.isnan(MLAPDV)): 1c

1351 _logger.warning('Areas of FOV lie outside the brain') 1c

1352 _logger.info(f'done ({time.time() - start_time:3.1f} seconds)\n') 1c

1353 MLAPDV = np.reshape(MLAPDV, [*x_px_idx.shape, 3]) 1c

1354 annotation = np.reshape(annotation, x_px_idx.shape) 1c

1355 

1356 fov['MLAPDV'] = { 1c

1357 'topLeft': MLAPDV[0, 0, :].tolist(), 

1358 'topRight': MLAPDV[0, -1, :].tolist(), 

1359 'bottomLeft': MLAPDV[-1, 0, :].tolist(), 

1360 'bottomRight': MLAPDV[-1, -1, :].tolist(), 

1361 'center': MLAPDV[round(x_px_idx.shape[0] / 2) - 1, round(x_px_idx.shape[1] / 2) - 1, :].tolist() 

1362 } 

1363 

1364 # Save the brain regions of the corners/centers of FOV (annotation field) 

1365 fov['brainLocationIds'] = { 1c

1366 'topLeft': int(annotation[0, 0]), 

1367 'topRight': int(annotation[0, -1]), 

1368 'bottomLeft': int(annotation[-1, 0]), 

1369 'bottomRight': int(annotation[-1, -1]), 

1370 'center': int(annotation[round(x_px_idx.shape[0] / 2) - 1, round(x_px_idx.shape[1] / 2) - 1]) 

1371 } 

1372 

1373 mlapdv[i] = MLAPDV 1c

1374 location_id[i] = annotation 1c

1375 return mlapdv, location_id 1c

1376 

1377 

1378def surface_normal(triangle): 

1379 """ 

1380 Calculate the surface normal unit vector of one or more triangles. 

1381 

1382 Parameters 

1383 ---------- 

1384 triangle : numpy.array 

1385 An array of shape (n_triangles, 3, 3) representing (Px Py Pz). 

1386 

1387 Returns 

1388 ------- 

1389 numpy.array 

1390 The surface normal unit vector(s). 

1391 """ 

1392 if triangle.shape == (3, 3): 1pc

1393 triangle = triangle[np.newaxis, :, :] 1pc

1394 if triangle.shape[1:] != (3, 3): 1pc

1395 raise ValueError('expected array of shape (3, 3); 3 coordinates in x, y, and z') 1p

1396 V = triangle[:, 1, :] - triangle[:, 0, :] # V = P2 - P1 1pc

1397 W = triangle[:, 2, :] - triangle[:, 0, :] # W = P3 - P1 1pc

1398 

1399 Nx = (V[:, 1] * W[:, 2]) - (V[:, 2] * W[:, 1]) # Nx = (Vy * Wz) - (Vz * Wy) 1pc

1400 Ny = (V[:, 2] * W[:, 0]) - (V[:, 0] * W[:, 2]) # Ny = (Vz * Wx) - (Vx * Wz) 1pc

1401 Nz = (V[:, 0] * W[:, 1]) - (V[:, 1] * W[:, 0]) # Nz = (Vx * Wy) - (Vy * Wx) 1pc

1402 N = np.c_[Nx, Ny, Nz] 1pc

1403 # Calculate unit vector. Transpose allows vectorized operation. 

1404 A = N / np.sqrt((Nx ** 2) + (Ny ** 2) + (Nz ** 2))[np.newaxis].T 1pc

1405 return A.squeeze() 1pc

1406 

1407 

1408@nb.njit('b1(f8[:,:], f8[:])') 

1409def in_triangle(triangle, point): 

1410 """ 

1411 Check whether `point` lies within `triangle`. 

1412 

1413 Parameters 

1414 ---------- 

1415 triangle : numpy.array 

1416 A (2 x 3) array of x-y coordinates; A(x1, y1), B(x2, y2) and C(x3, y3). 

1417 point : numpy.array 

1418 A point, P(x, y). 

1419 

1420 Returns 

1421 ------- 

1422 bool 

1423 True if coordinate lies within triangle. 

1424 """ 

1425 def area(x1, y1, x2, y2, x3, y3): 

1426 """Calculate the area of a triangle, given its vertices.""" 

1427 return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2.) 

1428 

1429 x1, y1, x2, y2, x3, y3 = triangle.flat 

1430 x, y = point 

1431 A = area(x1, y1, x2, y2, x3, y3) # area of triangle ABC 

1432 A1 = area(x, y, x2, y2, x3, y3) # area of triangle PBC 

1433 A2 = area(x1, y1, x, y, x3, y3) # area of triangle PAC 

1434 A3 = area(x1, y1, x2, y2, x, y) # area of triangle PAB 

1435 # Check if sum of A1, A2 and A3 equals that of A 

1436 diff = np.abs((A1 + A2 + A3) - A) 

1437 REL_TOL = 1e-9 

1438 return diff <= np.abs(REL_TOL * A) # isclose not yet implemented in numba 0.57 

1439 

1440 

1441@nb.njit('i8(f8[:], f8[:,:], intp[:,:])', nogil=True) 

1442def find_triangle(point, vertices, connectivity_list): 

1443 """ 

1444 Find which vertices contain a given point. 

1445 

1446 Currently O(n) but could take advantage of connectivity order to be quicker. 

1447 

1448 Parameters 

1449 ---------- 

1450 point : numpy.array 

1451 The (x, y) coordinate of a point to locate within one of the triangles. 

1452 vertices : numpy.array 

1453 An N x 3 array of vertices representing a triangle mesh. 

1454 connectivity_list : numpy.array 

1455 An N x 3 array of indices representing the connectivity of `points`. 

1456 

1457 Returns 

1458 ------- 

1459 int 

1460 The index of the vertices containing `point`, or -1 if not within any triangle. 

1461 """ 

1462 face_ind = -1 

1463 for i in nb.prange(connectivity_list.shape[0]): 

1464 triangle = vertices[connectivity_list[i, :], :] 

1465 if in_triangle(triangle, point): 

1466 face_ind = i 

1467 break 

1468 return face_ind 

1469 

1470 

1471@nb.njit('Tuple((f8[:], intp[:]))(f8[:], f8[:])', nogil=True) 

1472def _nearest_neighbour_1d(x, x_new): 

1473 """ 

1474 Nearest neighbour interpolation with extrapolation. 

1475 

1476 This was adapted from scipy.interpolate.interp1d but returns the indices of each nearest x 

1477 value. Assumes x is not sorted. 

1478 

1479 Parameters 

1480 ---------- 

1481 x : (N,) array_like 

1482 A 1-D array of real values. 

1483 x_new : (N,) array_like 

1484 A 1D array of values to apply function to. 

1485 

1486 Returns 

1487 ------- 

1488 numpy.array 

1489 A 1D array of interpolated values. 

1490 numpy.array 

1491 A 1D array of indices. 

1492 """ 

1493 SIDE = 'left' # use 'right' to round up to nearest int instead of rounding down 

1494 # Sort values 

1495 ind = np.argsort(x, kind='mergesort') 

1496 x = x[ind] 

1497 x_bds = x / 2.0 # Do division before addition to prevent possible integer overflow 

1498 x_bds = x_bds[1:] + x_bds[:-1] 

1499 # Find where in the averaged data the values to interpolate would be inserted. 

1500 x_new_indices = np.searchsorted(x_bds, x_new, side=SIDE) 

1501 # Clip x_new_indices so that they are within the range of x indices. 

1502 x_new_indices = x_new_indices.clip(0, len(x) - 1).astype(np.intp) 

1503 # Calculate the actual value for each entry in x_new. 

1504 y_new = x[x_new_indices] 

1505 return y_new, ind[x_new_indices] 

1506 

1507 

1508@nb.njit('Tuple((f8[:,:], u2[:]))(f8[:], f8[:], f8[:,:], f8[:], f8[:], f8[:], u2[:,:,:])', nogil=True) 

1509def _update_points(t, normal_vector, coords, axis_ml_um, axis_ap_um, axis_dv_um, atlas_labels): 

1510 """ 

1511 Determine the MLAPDV coordinate and brain location index for each of the given coordinates. 

1512 

1513 This has been optimized in numba. The majority of the time savings come from replacing iterp1d 

1514 and ismember with _nearest_neighbour_1d which were extremely slow. Parallel iteration further 

1515 halved the time it took per 512x512 FOV. 

1516 

1517 Parameters 

1518 ---------- 

1519 t : numpy.array 

1520 An N x 3 evenly spaced set of coordinates representing points going down from the coverslip 

1521 towards the brain. 

1522 normal_vector : numpy.array 

1523 The unit vector of the face normal to the center of the window. 

1524 coords : numpy.array 

1525 A set of N x 3 coordinates representing the MLAPDV coordinates of each pixel relative to 

1526 the center of the window, in micrometers (um). 

1527 axis_ml_um : numpy.array 

1528 An evenly spaced array of medio-lateral brain coordinates relative to bregma in um, at the 

1529 resolution of the atlas image used. 

1530 axis_ap_um : numpy.array 

1531 An evenly spaced array of anterio-posterior brain coordinates relative to bregma in um, at 

1532 the resolution of the atlas image used. 

1533 axis_dv_um : numpy.array 

1534 An evenly spaced array of dorso-ventral brain coordinates relative to bregma in um, at 

1535 the resolution of the atlas image used. 

1536 atlas_labels : numpy.array 

1537 A 3D array of integers representing the brain location index of each voxel of a given 

1538 atlas. The shape is expected to be (nAP, nML, nDV). 

1539 

1540 Returns 

1541 ------- 

1542 numpy.array 

1543 An N by 3 array containing the MLAPDV coordinates in um of each pixel coordinate. 

1544 Coordinates outside of the brain are NaN. 

1545 numpy.array 

1546 A 1D array of atlas label indices the length of `coordinates`. 

1547 """ 

1548 # passing through the center of the craniotomy/coverslip 

1549 traj_coords_centered = np.outer(t, -normal_vector) 

1550 MLAPDV = np.full_like(coords, np.nan) 

1551 annotation = np.zeros(coords.shape[0], dtype=np.uint16) 

1552 n_points = coords.shape[0] 

1553 for p in nb.prange(n_points): 

1554 # Shifted to the correct point on the coverslip, in true ML-AP-DV coords 

1555 traj_coords = traj_coords_centered + coords[p, :] 

1556 

1557 # Find intersection coordinate with the brain. 

1558 # Only use coordinates that exist in the atlas (kind of nearest neighbour interpolation) 

1559 ml, ml_idx = _nearest_neighbour_1d(axis_ml_um, traj_coords[:, 0]) 

1560 ap, ap_idx = _nearest_neighbour_1d(axis_ap_um, traj_coords[:, 1]) 

1561 dv, dv_idx = _nearest_neighbour_1d(axis_dv_um, traj_coords[:, 2]) 

1562 

1563 # Iterate over coordinates to find the first (if any) that is within the brain 

1564 ind = -1 

1565 area = 0 # 0 = void; 1 = root 

1566 for i in nb.prange(traj_coords.shape[0]): 

1567 anno = atlas_labels[ap_idx[i], ml_idx[i], dv_idx[i]] 

1568 if anno > 0: # first coordinate in the brain 

1569 ind = i 

1570 area = anno 

1571 if area > 1: # non-root brain area; we're done 

1572 break 

1573 if area > 1: 

1574 point = traj_coords[ind, :] 

1575 MLAPDV[p, :] = point # in um 

1576 annotation[p] = area 

1577 else: 

1578 MLAPDV[p, :] = np.nan 

1579 annotation[p] = area # root or void 

1580 

1581 return MLAPDV, annotation