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
« prev ^ index » next coverage.py v7.9.1, created at 2025-07-02 18:55 +0100
1"""The mesoscope data extraction pipeline.
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.
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
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
40from ibllib.pipes import base_tasks
41from ibllib.oneibl.data_handlers import ExpectedDataset, dataset_from_name
42from ibllib.io.extractors import mesoscope
45_logger = logging.getLogger(__name__)
46Provenance = enum.Enum('Provenance', ['ESTIMATE', 'FUNCTIONAL', 'LANDMARK', 'HISTOLOGY']) # py3.11 make StrEnum
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'
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
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_??'))
70 def _run(self):
71 """
72 Assert one reference image per collection and rename it. Register snapshots.
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
89class MesoscopeCompress(base_tasks.MesoscopeTask):
90 """ Tar compress raw 2p tif files, optionally remove uncompressed data."""
92 priority = 90
93 job_size = 'large'
94 _log_level = None
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
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
110 def tearDown(self):
111 _logger.setLevel(self._log_level or logging.INFO) 1f
112 return super().tearDown() 1f
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.
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.
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
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
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 )
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
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)
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
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
200 return outfiles 1f
203class MesoscopePreprocess(base_tasks.MesoscopeTask):
204 """Run suite2p preprocessing on tif files."""
206 priority = 80
207 cpu = -1
208 job_size = 'large'
209 env = 'suite2p'
211 def __init__(self, *args, **kwargs):
212 self._teardown_files = [] 1buaed
213 self.overwrite = False 1buaed
214 super().__init__(*args, **kwargs) 1buaed
216 def setUp(self, **kwargs):
217 """Set up task.
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.
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
257 def tearDown(self):
258 """Tear down task.
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
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
301 return signature 1hmaed
303 @staticmethod
304 def _masks2sparse(stat, ops):
305 """
306 Extract 3D sparse mask arrays from suit2p output.
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').
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.
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
341 def _rename_outputs(self, suite2p_dir, frameQC_names, frameQC, rename_dict=None):
342 """
343 Convert suite2p output files to ALF datasets.
345 This also moves any data.bin and ops.npy files to raw_bin_files for quicker re-runs.
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.
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
420 def load_meta_files(self):
421 """Load the extracted imaging metadata files.
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.
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
443 @staticmethod
444 def _consolidate_metadata(meta_data_all: list) -> dict:
445 """
446 Check that the metadata is consistent across all raw imaging folders.
448 Parameters
449 ----------
450 meta_data_all: list of dicts
451 List of metadata dictionaries to be checked for consistency.
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
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
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
478 @staticmethod
479 def _consolidate_exptQC(exptQC):
480 """
481 Consolidate exptQC.mat files into a single file.
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'}.
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
536 def get_default_tau(self):
537 """
538 Determine the tau (fluorescence decay) from the subject's genotype.
540 Returns
541 -------
542 float
543 The tau value to use.
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
557 def _meta2ops(self, meta):
558 """
559 Create the ops dictionary for suite2p based on metadata.
561 Parameters
562 ----------
563 meta: dict
564 Imaging metadata.
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
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')
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
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
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))
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 }
631 return db 1ka
633 @staticmethod
634 def _get_plane_paths(path):
635 """Return list of sorted suite2p plane folder paths.
637 Parameters
638 ----------
639 path : pathlib.Path
640 The path containing plane folders.
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
650 def bin_per_plane(self, metadata, **kwargs):
651 """
652 Extracts a binary data file of imaging data per imaging plane.
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.
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
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
686 ret = suite2p.io.mesoscan_to_binary(ops.copy()) 1a
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
694 return planes, ret 1a
696 def image_motion_registration(self, ops):
697 """Perform motion registration.
699 Parameters
700 ----------
701 ops : dict
702 A dict of suite2p options.
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".
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
725 def roi_detection(self, ops):
726 """Perform ROI detection.
728 Parameters
729 ----------
730 ops : dict
731 A dict of suite2p options.
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
744 def _run(self, rename_files=True, use_badframes=True, **kwargs):
745 """
746 Process inputs, run suite2p and make outputs alf compatible.
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/
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.
765 Returns
766 -------
767 list of pathlib.Path
768 All files created by the task.
770 """
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
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
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
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
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'])
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
826 """ Suite2p """ 1a
827 # Create alf if is doesn't exist
828 self.session_path.joinpath('alf').mkdir(exist_ok=True) 1a
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
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
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
857 return list(out_files) 1a
860class MesoscopeSync(base_tasks.MesoscopeTask):
861 """Extract the frame times from the main DAQ."""
863 priority = 40
864 job_size = 'small'
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
881 def _run(self, **kwargs):
882 """
883 Extract the imaging times for all FOVs.
885 Returns
886 -------
887 list of pathlib.Path
888 Files containing frame timestamps for individual FOVs and time offsets for each line scan.
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
914class MesoscopeFOV(base_tasks.MesoscopeTask):
915 """Create FOV and FOV location objects in Alyx from metadata."""
917 priority = 40
918 job_size = 'small'
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
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.
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.
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.
949 Returns
950 -------
951 dict
952 The newly created FOV Alyx record.
953 list
954 The newly created FOV location Alyx records.
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
968 suffix = None if provenance is Provenance.HISTOLOGY else provenance.name.lower() 1g
969 _logger.info('Extracting %s MLAPDV datasets', suffix or 'final') 1g
971 # Extract mean image MLAPDV coordinates and brain location IDs
972 mean_image_mlapdv, mean_image_ids = self.project_mlapdv(meta) 1g
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
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
989 # Extract ROI MLAPDV coordinates and brain location IDs
990 roi_mlapdv, roi_brain_ids = self.roi_mlapdv(nFOV, suffix=suffix) 1g
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
1002 # Register FOVs in Alyx
1003 self.register_fov(meta, suffix) 1g
1005 return sorted([*meta_files, *roi_files, *mean_image_files]) 1g
1007 def update_surgery_json(self, meta, normal_vector):
1008 """
1009 Update surgery JSON with surface normal vector.
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.
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.
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
1047 def roi_mlapdv(self, nFOV: int, suffix=None):
1048 """
1049 Extract ROI MLAPDV coordinates and brain location IDs.
1051 MLAPDV coordinates are in um relative to bregma. Location IDs are from the 2017 Allen
1052 common coordinate framework atlas.
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'.
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
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
1080 # Load MLAPDV + brain location ID maps of pixels
1081 mpciMeanImage = alfio.load_object( 1g
1082 alf_path, 'mpciMeanImage', attribute=['mlapdv', 'brainLocationIds'])
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
1096 return all_mlapdv, all_brain_ids 1g
1098 @staticmethod
1099 def get_provenance(filename):
1100 """
1101 Get the field of view provenance from a mpciMeanImage or mpciROIs dataset.
1103 Parameters
1104 ----------
1105 filename : str, pathlib.Path
1106 A filename to get the provenance from.
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
1119 def register_fov(self, meta: dict, suffix: str = None) -> (list, list):
1120 """
1121 Create FOV on Alyx.
1123 Assumes field of view recorded perpendicular to objective.
1124 Assumes field of view is plane (negligible volume).
1126 Required Alyx fixtures:
1127 - experiments.ImagingType(name='mesoscope')
1128 - experiments.CoordinateSystem(name='IBL-Allen')
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).
1138 Returns
1139 -------
1140 list of dict
1141 A list registered of field of view entries from Alyx.
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}
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
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
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
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
1191 data['brain_region'] = np.unique(mean_image_ids).astype(int).tolist() 1i
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
1200 def load_triangulation(self):
1201 """
1202 Load the surface triangulation file.
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.
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
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.
1228 MLAPDV coordinates are in um relative to bregma. Location IDs are from the 2017 Allen
1229 common coordinate framework atlas.
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.
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
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
1269 dorsal_triangle = points[dorsal_connectivity_list[face_ind, :], :] 1c
1271 # Get the surface normal unit vector of dorsal triangle
1272 normal_vector = surface_normal(dorsal_triangle) 1c
1274 # Update the surgery JSON field with normal unit vector, for use in histology alignment
1275 self.update_surgery_json(meta, normal_vector) 1c
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
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
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
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)
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
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
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
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
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
1333 xx = xx.flatten() - coord_ml 1c
1334 yy = yy.flatten() - coord_ap 1c
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
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
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
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
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 }
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 }
1373 mlapdv[i] = MLAPDV 1c
1374 location_id[i] = annotation 1c
1375 return mlapdv, location_id 1c
1378def surface_normal(triangle):
1379 """
1380 Calculate the surface normal unit vector of one or more triangles.
1382 Parameters
1383 ----------
1384 triangle : numpy.array
1385 An array of shape (n_triangles, 3, 3) representing (Px Py Pz).
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
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
1408@nb.njit('b1(f8[:,:], f8[:])')
1409def in_triangle(triangle, point):
1410 """
1411 Check whether `point` lies within `triangle`.
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).
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.)
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
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.
1446 Currently O(n) but could take advantage of connectivity order to be quicker.
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`.
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
1471@nb.njit('Tuple((f8[:], intp[:]))(f8[:], f8[:])', nogil=True)
1472def _nearest_neighbour_1d(x, x_new):
1473 """
1474 Nearest neighbour interpolation with extrapolation.
1476 This was adapted from scipy.interpolate.interp1d but returns the indices of each nearest x
1477 value. Assumes x is not sorted.
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.
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]
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.
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.
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).
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, :]
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])
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
1581 return MLAPDV, annotation