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

673 statements  

« prev     ^ index     » next       coverage.py v7.7.0, created at 2025-03-17 15:25 +0000

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 

17import shutil 

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 continue 

149 if not infiles: 1f

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

151 continue 

152 

153 _logger.debug( 1f

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

155 ) 

156 

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

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

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

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

161 _logger.debug(cmd) 1f

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

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

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

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

166 

167 # Check the output 

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

169 outfiles.append(outfile) 1f

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

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

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

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

174 uncompressed_size / compressed_size, 

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

176 (uncompressed_size - compressed_size) / 1024 / 1024) 

177 

178 if verify_output: 1f

179 # Test bzip 

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

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

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

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

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

185 # Check tar 

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

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

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

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

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

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

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

193 

194 if remove_uncompressed: 1f

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

196 for file in infiles: 1f

197 file.unlink() 1f

198 

199 return outfiles 1f

200 

201 

202class MesoscopePreprocess(base_tasks.MesoscopeTask): 

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

204 

205 priority = 80 

206 cpu = -1 

207 job_size = 'large' 

208 env = 'suite2p' 

209 

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

211 self._teardown_files = [] 1buaed

212 self.overwrite = False 1buaed

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

214 

215 def setUp(self, **kwargs): 

216 """Set up task. 

217 

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

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

220 the compressed input. 

221 

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

223 """ 

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

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

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

227 if not self.overwrite and all(x.find_files(self.session_path)[0] for x in bin_sig): 1hm

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

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

230 if not tif_sig: 1hm

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

232 tif_sig = tif_sig[0] 1hm

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

234 if tifs_present or not all_files_present: 1hm

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

236 # Decompress imaging files 

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

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

239 if not all(present): 1h

240 return False # Compressed files missing; return 

241 files = flatten(files) 1h

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

243 for file in files: 1h

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

245 _logger.debug(cmd) 1h

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

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

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

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

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

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

252 self._teardown_files.extend(tifs) 1h

253 return all_files_present 1h

254 

255 def tearDown(self): 

256 """Tear down task. 

257 

258 This removes any decompressed tif files. 

259 """ 

260 for file in self._teardown_files: 1hm

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

262 file.unlink() 1h

263 return super().tearDown() 1hm

264 

265 @property 

266 def signature(self): 

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

268 I = ExpectedDataset.input # noqa 1hmaed

269 signature = { 1hmaed

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

289 } 

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

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

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

293 registered_bin = I('data.bin', 'raw_bin_files/plane*', True, unique=False) 1hmaed

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

295 

296 return signature 1hmaed

297 

298 @staticmethod 

299 def _masks2sparse(stat, ops): 

300 """ 

301 Extract 3D sparse mask arrays from suit2p output. 

302 

303 Parameters 

304 ---------- 

305 stat : numpy.array 

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

307 ops : numpy.array 

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

309 

310 Returns 

311 ------- 

312 sparse.GCXS 

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

314 sparse.GCXS 

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

316 

317 Notes 

318 ----- 

319 Save using sparse.save_npz. 

320 """ 

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

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

323 coords = [[], [], []] 1ed

324 data = [] 1ed

325 pil_coords = [] 1ed

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

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

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

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

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

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

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

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

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

335 

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

337 """ 

338 Convert suite2p output files to ALF datasets. 

339 

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

341 

342 Parameters 

343 ---------- 

344 suite2p_dir : pathlib.Path 

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

346 rename_dict : dict or None 

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

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

349 

350 Returns 

351 ------- 

352 list of pathlib.Path 

353 All paths found in FOV folders. 

354 """ 

355 if rename_dict is None: 1ed

356 rename_dict = { 1ed

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

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

359 'Fneu.npy': 'mpci.ROINeuropilActivityF.npy' 

360 } 

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

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

363 # Move bin file(s) out of the way 

364 bin_files = list(plane_dir.glob('data*.bin')) # e.g. data.bin, data_raw.bin, data_chan2_raw.bin 1ed

365 if any(bin_files): 1ed

366 (bin_files_dir := self.session_path.joinpath('raw_bin_files', plane_dir.name)).mkdir(parents=True, exist_ok=True) 1d

367 _logger.debug('Moving bin file(s) to %s', bin_files_dir.relative_to(self.session_path)) 1d

368 for bin_file in bin_files: 1d

369 dst = bin_files_dir.joinpath(bin_file.name) 1d

370 bin_file.replace(dst) 1d

371 # copy ops file for lazy re-runs 

372 shutil.copy(plane_dir.joinpath('ops.npy'), bin_files_dir.joinpath('ops.npy')) 1d

373 # Archive the raw suite2p output before renaming 

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

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

376 if fov_dir.exists(): 1ed

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

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

379 f.unlink() 1d

380 if not any(fov_dir.iterdir()): 1d

381 _logger.debug('Removing old folder %s', fov_dir.relative_to(self.session_path)) 

382 fov_dir.rmdir() 

383 prev_level = _logger.level 1ed

384 _logger.setLevel(logging.WARNING) 1ed

385 shutil.make_archive(str(fov_dir / '_suite2p_ROIData.raw'), 'zip', plane_dir, logger=_logger) 1ed

386 _logger.setLevel(prev_level) 1ed

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

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

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

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

391 # extract some other data from suite2p outputs 

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

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

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

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

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

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

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

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

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

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

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

403 # clusters uuids 

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

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

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

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

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

409 # ROI and neuropil masks 

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

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

412 sparse.save_npz(fp, roi_mask) 1ed

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

414 sparse.save_npz(fp, pil_mask) 1ed

415 # Remove old suite2p files 

416 shutil.rmtree(str(suite2p_dir), ignore_errors=False, onerror=None) 1ed

417 # Collect all files in those directories 

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

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

420 

421 def load_meta_files(self): 

422 """Load the extracted imaging metadata files. 

423 

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

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

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

427 

428 Returns 

429 ------- 

430 dict 

431 Single, consolidated dictionary containing metadata. 

432 list of dict 

433 The meta data for each individual imaging bout. 

434 """ 

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

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

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

438 # Check there is exactly 1 meta file per collection 

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

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

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

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

443 

444 @staticmethod 

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

446 """ 

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

448 

449 Parameters 

450 ---------- 

451 meta_data_all: list of dicts 

452 List of metadata dictionaries to be checked for consistency. 

453 

454 Returns 

455 ------- 

456 dict 

457 Single, consolidated dictionary containing metadata. 

458 """ 

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

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

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

462 

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

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

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

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

467 

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

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

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

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

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

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

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

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

476 f'Using meta_data from first folder!') 

477 return meta_data_all[0] 1a

478 

479 @staticmethod 

480 def _consolidate_exptQC(exptQC): 

481 """ 

482 Consolidate exptQC.mat files into a single file. 

483 

484 Parameters 

485 ---------- 

486 exptQC : list of pandas.DataFrame 

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

488 {'frameQC_frames', 'frameQC_names'}. 

489 

490 Returns 

491 ------- 

492 numpy.array 

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

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

495 pandas.DataFrame 

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

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

498 numpy.array 

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

500 """ 

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

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

503 qc_labels = ['ok'] 1ja

504 frame_qc = [] 1ja

505 for e in exptQC: 1ja

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

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

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

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

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

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

512 for name in names: 1ja

513 i_old = names.index(name) # old enumeration 1ja

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

515 try: 1ja

516 i_new = qc_labels.index(name) 1ja

517 except ValueError: 1ja

518 i_new = len(qc_labels) 1ja

519 qc_labels.append(name) 1ja

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

521 frame_qc.append(frames) 1ja

522 # Concatenate frames 

523 frame_qc = np.concatenate(frame_qc) 1ja

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

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

526 try: 1j

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

528 except ValueError: 

529 i = len(qc_labels) 

530 qc_labels.append('unknown') 

531 frame_qc[missing_name] = i 1j

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

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

534 # Convert labels to value -> label data frame 

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

536 return frame_qc, frame_qc_names, bad_frames 1ja

537 

538 def get_default_tau(self): 

539 """ 

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

541 

542 Returns 

543 ------- 

544 float 

545 The tau value to use. 

546 

547 See Also 

548 -------- 

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

550 """ 

551 # These settings are from the suite2P documentation 

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

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

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

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

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

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

558 

559 def _meta2ops(self, meta): 

560 """ 

561 Create the ops dictionary for suite2p based on metadata. 

562 

563 Parameters 

564 ---------- 

565 meta: dict 

566 Imaging metadata. 

567 

568 Returns 

569 ------- 

570 dict 

571 Inputs to suite2p run that deviate from default parameters. 

572 """ 

573 # Computing dx and dy 

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

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

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

577 

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

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

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

581 

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

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

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

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

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

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

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

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

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

591 

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

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

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

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

596 

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

598 # FIXME currently unused 

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

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

601 

602 db = { 1ka

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

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

605 'fast_disk': '', # TODO 

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

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

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

609 'keep_movie_raw': False, 

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

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

612 'nimg_init': 400, 

613 'combined': False, 

614 'nonrigid': True, 

615 'maxregshift': 0.05, # default = 1 

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

617 'block_size': [128, 128], 

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

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

620 'mesoscan': True, 

621 'nplanes': nplanes, 

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

623 'nchannels': nchannels, 

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

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

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

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

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

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

630 'dx': dx.tolist(), 

631 'dy': dy.tolist() 

632 } 

633 

634 return db 1ka

635 

636 @staticmethod 

637 def _get_plane_paths(path): 

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

639 

640 Parameters 

641 ---------- 

642 path : pathlib.Path 

643 The path containing plane folders. 

644 

645 Returns 

646 ------- 

647 list of pathlib.Path 

648 The plane folder paths, ordered by number. 

649 """ 

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

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

652 

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

654 """ 

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

656 

657 Parameters 

658 ---------- 

659 metadata : dict 

660 A dictionary of extracted metadata. 

661 save_path0 : str, pathlib.Path 

662 The root path of the suite2p bin output. 

663 save_folder : str 

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

665 kwargs 

666 Other optional arguments to overwrite the defaults for. 

667 

668 Returns 

669 ------- 

670 list of pathlib.Path 

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

672 dict 

673 Suite2p's modified options. 

674 """ 

675 import suite2p.io 1a

676 

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

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

679 'do_registration', 'slices') 

680 ops = self._meta2ops(metadata) 1a

681 ops['force_sktiff'] = False 1a

682 ops['do_registration'] = True 1a

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

684 ops.update(kwargs) 1a

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

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

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

688 

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

690 

691 # Get ordered list of plane folders 

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

693 assert out_path.exists() 1a

694 planes = self._get_plane_paths(out_path) 1a

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

696 

697 return planes, ret 1a

698 

699 def image_motion_registration(self, ops): 

700 """Perform motion registration. 

701 

702 Parameters 

703 ---------- 

704 ops : dict 

705 A dict of suite2p options. 

706 

707 Returns 

708 ------- 

709 dict 

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

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

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

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

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

715 

716 """ 

717 import suite2p 1r

718 ops['do_registration'] = True 1r

719 ops['do_regmetrics'] = True 1r

720 ops['roidetect'] = False 1r

721 ret = suite2p.run_plane(ops) 1r

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

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

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

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

726 return metrics 1r

727 

728 def roi_detection(self, ops): 

729 """Perform ROI detection. 

730 

731 Parameters 

732 ---------- 

733 ops : dict 

734 A dict of suite2p options. 

735 

736 Returns 

737 ------- 

738 dict 

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

740 """ 

741 import suite2p 1v

742 ops['do_registration'] = False 1v

743 ops['roidetect'] = True 1v

744 ret = suite2p.run_plane(ops) 1v

745 return ret 1v

746 

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

748 """ 

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

750 

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

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

753 to 'alf/ 

754 

755 Parameters 

756 ---------- 

757 rename_files: bool 

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

759 use_badframes: bool 

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

761 overwrite : bool 

762 Whether to re-perform extraction and motion registration. 

763 do_registration : bool 

764 Whether to perform motion registration. 

765 roidetect : bool 

766 Whether to perform ROI detection. 

767 

768 Returns 

769 ------- 

770 list of pathlib.Path 

771 All files created by the task. 

772 

773 """ 

774 

775 """ Metadata and parameters """ 1a

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

777 # Load and consolidate the image metadata from JSON files 

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

779 

780 # Create suite2p output folder in root session path 

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

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

783 

784 # Check for previous intermediate files 

785 plane_folders = self._get_plane_paths(save_path) 1a

786 if len(plane_folders) == 0 and self.session_path.joinpath('raw_bin_files').exists(): 1a

787 self.session_path.joinpath('raw_bin_files').replace(save_path) 1a

788 plane_folders = self._get_plane_paths(save_path) 1a

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

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

791 # Ingest tiff files 

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

793 

794 """ Bad frames """ 1a

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

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

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

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

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

800 if len(exptQC) > 0: 1a

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

802 else: 

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

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

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

806 

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

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

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

810 total_frames = 0 1a

811 # Ensure all indices are relative to total cumulative frames 

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

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

814 if badframes_path.exists(): 1a

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

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

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

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

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

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

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

822 

823 """ Suite2p """ 1a

824 # Create alf if is doesn't exist 

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

826 

827 # Perform registration 

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

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

830 for plane in plane_folders: 1a

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

832 ops.update(kwargs) 1a

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

834 _ = self.image_motion_registration(ops) 1a

835 # TODO Handle metrics and QC here 

836 

837 # ROI detection 

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

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

840 for plane in plane_folders: 1a

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

842 ops.update(kwargs) 1a

843 self.roi_detection(ops) 1a

844 

845 """ Outputs """ 1a

846 # Save and rename other outputs 

847 if rename_files: 1a

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

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

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

851 else: 

852 out_files = save_path.rglob('*.*') 1a

853 

854 return list(out_files) 1a

855 

856 

857class MesoscopeSync(base_tasks.MesoscopeTask): 

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

859 

860 priority = 40 

861 job_size = 'small' 

862 

863 @property 

864 def signature(self): 

865 I = ExpectedDataset.input # noqa 1lo

866 signature = { 1lo

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

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

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

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

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

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

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

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

875 } 

876 return signature 1lo

877 

878 def _run(self): 

879 """ 

880 Extract the imaging times for all FOVs. 

881 

882 Returns 

883 ------- 

884 list of pathlib.Path 

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

886 

887 """ 

888 # TODO function to determine nFOVs 

889 try: 1lo

890 alf_path = self.session_path / self.sync_collection 1lo

891 events = alfio.load_object(alf_path, 'softwareEvents').get('log') 1lo

892 except alferr.ALFObjectNotFound: 1l

893 events = None 1l

894 if events is None or events.empty: 1lo

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

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

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

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

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

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

901 self.rawImagingData['meta'] = mesoscope.patch_imaging_meta(self.rawImagingData['meta']) 1lo

902 n_FOVs = len(self.rawImagingData['meta']['FOV']) 1lo

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

904 mesosync = mesoscope.MesoscopeSyncTimeline(self.session_path, n_FOVs) 1lo

905 _, out_files = mesosync.extract( 1lo

906 save=True, sync=sync, chmap=chmap, device_collection=collections, events=events) 

907 return out_files 1lo

908 

909 

910class MesoscopeFOV(base_tasks.MesoscopeTask): 

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

912 

913 priority = 40 

914 job_size = 'small' 

915 

916 @property 

917 def signature(self): 

918 signature = { 1g

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

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

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

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

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

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

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

926 } 

927 return signature 1g

928 

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

930 """ 

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

932 

933 Steps: 

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

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

936 of each ROI. 

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

938 

939 Parameters 

940 ---------- 

941 provenance : Provenance 

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

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

944 

945 Returns 

946 ------- 

947 dict 

948 The newly created FOV Alyx record. 

949 list 

950 The newly created FOV location Alyx records. 

951 

952 Notes 

953 ----- 

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

955 task will result in an error. 

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

957 """ 

958 # Load necessary data 

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

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

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

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

963 

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

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

966 

967 # Extract mean image MLAPDV coordinates and brain location IDs 

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

969 

970 # Save the meta data file with new coordinate fields 

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

972 json.dump(meta, fp) 1g

973 

974 # Save the mean image datasets 

975 mean_image_files = [] 1g

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

977 for i in range(nFOV): 1g

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

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

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

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

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

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

984 

985 # Extract ROI MLAPDV coordinates and brain location IDs 

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

987 

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

989 roi_files = [] 1g

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

991 for i in range(nFOV): 1g

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

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

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

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

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

997 

998 # Register FOVs in Alyx 

999 self.register_fov(meta, suffix) 1g

1000 

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

1002 

1003 def update_surgery_json(self, meta, normal_vector): 

1004 """ 

1005 Update surgery JSON with surface normal vector. 

1006 

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

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

1009 in the provided meta file. 

1010 

1011 Parameters 

1012 ---------- 

1013 meta : dict 

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

1015 normal_vector : array_like 

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

1017 

1018 Returns 

1019 ------- 

1020 dict 

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

1022 """ 

1023 if not self.one or self.one.offline: 1nc

1024 _logger.warning('failed to update surgery JSON: ONE offline') 1nc

1025 return 1nc

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

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

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

1029 if not surgeries: 1n

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

1031 return 1n

1032 surgery = surgeries[0] # Check most recent surgery in list 1n

1033 center = (meta['centerMM']['ML'], meta['centerMM']['AP']) 1n

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

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

1036 if (key := next(match, None)) is None: 1n

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

1038 return surgery 1n

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

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

1041 return surgery 1n

1042 

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

1044 """ 

1045 Extract ROI MLAPDV coordinates and brain location IDs. 

1046 

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

1048 common coordinate framework atlas. 

1049 

1050 Parameters 

1051 ---------- 

1052 nFOV : int 

1053 The number of fields of view acquired. 

1054 suffix : {None, 'estimate'} 

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

1056 estimates, the suffix should be 'estimate'. 

1057 

1058 Returns 

1059 ------- 

1060 dict of int : numpy.array 

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

1062 dict of int : numpy.array 

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

1064 """ 

1065 all_mlapdv = {} 1g

1066 all_brain_ids = {} 1g

1067 for n in range(nFOV): 1g

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

1069 

1070 # Load neuron centroids in pixel space 

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

1072 if not stack_pos_file: 1g

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

1074 stack_pos = alfio.load_file_content(stack_pos_file) 1g

1075 

1076 # Load MLAPDV + brain location ID maps of pixels 

1077 mpciMeanImage = alfio.load_object( 1g

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

1079 

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

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

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

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

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

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

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

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

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

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

1090 all_mlapdv[n] = mlapdv 1g

1091 

1092 return all_mlapdv, all_brain_ids 1g

1093 

1094 @staticmethod 

1095 def get_provenance(filename): 

1096 """ 

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

1098 

1099 Parameters 

1100 ---------- 

1101 filename : str, pathlib.Path 

1102 A filename to get the provenance from. 

1103 

1104 Returns 

1105 ------- 

1106 Provenance 

1107 The provenance of the file. 

1108 """ 

1109 filename = Path(filename).name 1w

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

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

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

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

1114 

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

1116 """ 

1117 Create FOV on Alyx. 

1118 

1119 Assumes field of view recorded perpendicular to objective. 

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

1121 

1122 Required Alyx fixtures: 

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

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

1125 

1126 Parameters 

1127 ---------- 

1128 meta : dict 

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

1130 suffix : str 

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

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

1133 

1134 Returns 

1135 ------- 

1136 list of dict 

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

1138 

1139 TODO Determine dual plane ID for JSON field 

1140 """ 

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

1142 alyx_fovs = [] 1i

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

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

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

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

1147 if dry: 1i

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

1149 else: 

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

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

1152 

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

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

1155 # Field of view 

1156 alyx_FOV = { 1i

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

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

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

1160 } 

1161 if dry: 1i

1162 print(alyx_FOV) 1i

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

1164 alyx_fovs.append(alyx_FOV) 1i

1165 else: 

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

1167 

1168 # Field of view location 

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

1170 'default_provenance': True, 

1171 'coordinate_system': 'IBL-Allen', 

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

1173 if suffix: 1i

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

1175 

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

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

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

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

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

1181 

1182 # Load MLAPDV + brain location ID maps of pixels 

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

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

1185 mean_image_ids = alfio.load_file_content(filepath) 1i

1186 

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

1188 

1189 if dry: 1i

1190 print(data) 1i

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

1192 else: 

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

1194 return alyx_fovs 1i

1195 

1196 def load_triangulation(self): 

1197 """ 

1198 Load the surface triangulation file. 

1199 

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

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

1202 MATLAB. 

1203 

1204 Returns 

1205 ------- 

1206 points : numpy.array 

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

1208 are in um relative to the IBL bregma coordinates. 

1209 connectivity_list : numpy.array 

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

1211 """ 

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

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

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

1215 connectivity_list = surface_triangulation['connectivity_list'] 1c

1216 surface_triangulation.close() 1c

1217 return points, connectivity_list 1c

1218 

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

1220 """ 

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

1222 location IDs. 

1223 

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

1225 common coordinate framework atlas. 

1226 

1227 Parameters 

1228 ---------- 

1229 meta : dict 

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

1231 view. 

1232 atlas : ibllib.atlas.Atlas 

1233 An atlas instance. 

1234 

1235 Returns 

1236 ------- 

1237 dict 

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

1239 dict 

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

1241 """ 

1242 mlapdv = {} 1c

1243 location_id = {} 1c

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

1245 # more accurately represent the living brain. 

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

1247 # The centre of the craniotomy / imaging window 

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

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

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

1251 

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

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

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

1255 triangles = points[connectivity_list, :] 1c

1256 normals = surface_normal(triangles) 1c

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

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

1259 dorsal_connectivity_list = connectivity_list[up_faces, :] 1c

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

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

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

1263 assert face_ind != -1 1c

1264 

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

1266 

1267 # Get the surface normal unit vector of dorsal triangle 

1268 normal_vector = surface_normal(dorsal_triangle) 1c

1269 

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

1271 self.update_surgery_json(meta, normal_vector) 1c

1272 

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

1274 # the three vertices defining the triangle 

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

1276 

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

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

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

1280 

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

1282 # coordinate (which is ML-AP coordinate) 

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

1284 

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

1286 # DO NOT USE THIS: 

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

1288 

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

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

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

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

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

1294 # normalize and flip the sign if necessary 

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

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

1297 

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

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

1300 # Let's shift the coordinates relative to bregma 

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

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

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

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

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

1306 

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

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

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

1310 start_time = time.time() 1c

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

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

1313 

1314 # xx and yy are in mm in coverslip space 

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

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

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

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

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

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

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

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

1323 

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

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

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

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

1328 

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

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

1331 

1332 # rotate xx and yy in 3D 

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

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

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

1336 

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

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

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

1340 

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

1342 MLAPDV, annotation = _update_points( 1c

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

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

1345 

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

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

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

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

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

1351 

1352 fov['MLAPDV'] = { 1c

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

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

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

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

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

1358 } 

1359 

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

1361 fov['brainLocationIds'] = { 1c

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

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

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

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

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

1367 } 

1368 

1369 mlapdv[i] = MLAPDV 1c

1370 location_id[i] = annotation 1c

1371 return mlapdv, location_id 1c

1372 

1373 

1374def surface_normal(triangle): 

1375 """ 

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

1377 

1378 Parameters 

1379 ---------- 

1380 triangle : numpy.array 

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

1382 

1383 Returns 

1384 ------- 

1385 numpy.array 

1386 The surface normal unit vector(s). 

1387 """ 

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

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

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

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

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

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

1394 

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

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

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

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

1399 # Calculate unit vector. Transpose allows vectorized operation. 

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

1401 return A.squeeze() 1pc

1402 

1403 

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

1405def in_triangle(triangle, point): 

1406 """ 

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

1408 

1409 Parameters 

1410 ---------- 

1411 triangle : numpy.array 

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

1413 point : numpy.array 

1414 A point, P(x, y). 

1415 

1416 Returns 

1417 ------- 

1418 bool 

1419 True if coordinate lies within triangle. 

1420 """ 

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

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

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

1424 

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

1426 x, y = point 

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

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

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

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

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

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

1433 REL_TOL = 1e-9 

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

1435 

1436 

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

1438def find_triangle(point, vertices, connectivity_list): 

1439 """ 

1440 Find which vertices contain a given point. 

1441 

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

1443 

1444 Parameters 

1445 ---------- 

1446 point : numpy.array 

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

1448 vertices : numpy.array 

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

1450 connectivity_list : numpy.array 

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

1452 

1453 Returns 

1454 ------- 

1455 int 

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

1457 """ 

1458 face_ind = -1 

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

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

1461 if in_triangle(triangle, point): 

1462 face_ind = i 

1463 break 

1464 return face_ind 

1465 

1466 

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

1468def _nearest_neighbour_1d(x, x_new): 

1469 """ 

1470 Nearest neighbour interpolation with extrapolation. 

1471 

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

1473 value. Assumes x is not sorted. 

1474 

1475 Parameters 

1476 ---------- 

1477 x : (N,) array_like 

1478 A 1-D array of real values. 

1479 x_new : (N,) array_like 

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

1481 

1482 Returns 

1483 ------- 

1484 numpy.array 

1485 A 1D array of interpolated values. 

1486 numpy.array 

1487 A 1D array of indices. 

1488 """ 

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

1490 # Sort values 

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

1492 x = x[ind] 

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

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

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

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

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

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

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

1500 y_new = x[x_new_indices] 

1501 return y_new, ind[x_new_indices] 

1502 

1503 

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

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

1506 """ 

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

1508 

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

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

1511 halved the time it took per 512x512 FOV. 

1512 

1513 Parameters 

1514 ---------- 

1515 t : numpy.array 

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

1517 towards the brain. 

1518 normal_vector : numpy.array 

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

1520 coords : numpy.array 

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

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

1523 axis_ml_um : numpy.array 

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

1525 resolution of the atlas image used. 

1526 axis_ap_um : numpy.array 

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

1528 the resolution of the atlas image used. 

1529 axis_dv_um : numpy.array 

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

1531 the resolution of the atlas image used. 

1532 atlas_labels : numpy.array 

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

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

1535 

1536 Returns 

1537 ------- 

1538 numpy.array 

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

1540 Coordinates outside of the brain are NaN. 

1541 numpy.array 

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

1543 """ 

1544 # passing through the center of the craniotomy/coverslip 

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

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

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

1548 n_points = coords.shape[0] 

1549 for p in nb.prange(n_points): 

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

1551 traj_coords = traj_coords_centered + coords[p, :] 

1552 

1553 # Find intersection coordinate with the brain. 

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

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

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

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

1558 

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

1560 ind = -1 

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

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

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

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

1565 ind = i 

1566 area = anno 

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

1568 break 

1569 if area > 1: 

1570 point = traj_coords[ind, :] 

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

1572 annotation[p] = area 

1573 else: 

1574 MLAPDV[p, :] = np.nan 

1575 annotation[p] = area # root or void 

1576 

1577 return MLAPDV, annotation