diff --git a/python/lib/db/queries/meg_ctf_head_shape.py b/python/lib/db/queries/meg_ctf_head_shape.py new file mode 100644 index 000000000..43e7a5292 --- /dev/null +++ b/python/lib/db/queries/meg_ctf_head_shape.py @@ -0,0 +1,17 @@ +from pathlib import Path + +from sqlalchemy import select +from sqlalchemy.orm import Session as Database + +from lib.db.models.meg_ctf_head_shape_file import DbMegCtfHeadShapeFile + + +def try_get_meg_ctf_head_shape_file_with_path(db: Database, path: Path) -> DbMegCtfHeadShapeFile | None: + """ + Get a MEG CTF head shape file from the database using its path, or return `None` if no file was + found. + """ + + return db.execute(select(DbMegCtfHeadShapeFile) + .where(DbMegCtfHeadShapeFile.path == path) + ).scalar_one_or_none() diff --git a/python/lib/physio/ctf/head_shape.py b/python/lib/physio/ctf/head_shape.py new file mode 100644 index 000000000..fe92012e3 --- /dev/null +++ b/python/lib/physio/ctf/head_shape.py @@ -0,0 +1,41 @@ +from decimal import Decimal +from pathlib import Path + +from loris_bids_utils.meg.ctf.head_shape import MegCtfHeadShapeFile +from loris_utils.crypto import compute_file_blake2b_hash + +from lib.db.models.meg_ctf_head_shape_file import DbMegCtfHeadShapeFile +from lib.db.models.meg_ctf_head_shape_point import DbMegCtfHeadShapePoint +from lib.env import Env + + +def insert_meg_ctf_head_shape_file( + env: Env, + head_shape_file: MegCtfHeadShapeFile, + loris_head_shape_file_path: Path, +) -> DbMegCtfHeadShapeFile: + """ + Insert a MEG CTF head shape file into the LORIS database. + """ + + blake2b_hash = compute_file_blake2b_hash(head_shape_file.path) + + db_head_shape_file = DbMegCtfHeadShapeFile( + path = loris_head_shape_file_path, + blake2b_hash = blake2b_hash, + ) + + env.db.add(db_head_shape_file) + env.db.flush() + + for name, point in head_shape_file.points.items(): + env.db.add(DbMegCtfHeadShapePoint( + file_id = db_head_shape_file.id, + name = name, + x = Decimal(point.x), + y = Decimal(point.y), + z = Decimal(point.z), + )) + + env.db.flush() + return db_head_shape_file diff --git a/python/lib/physio/file.py b/python/lib/physio/file.py index e1891da3f..ebf9b5c90 100644 --- a/python/lib/physio/file.py +++ b/python/lib/physio/file.py @@ -3,6 +3,7 @@ from pathlib import Path from lib.db.models.imaging_file_type import DbImagingFileType +from lib.db.models.meg_ctf_head_shape_file import DbMegCtfHeadShapeFile from lib.db.models.physio_file import DbPhysioFile from lib.db.models.physio_modality import DbPhysioModality from lib.db.models.physio_output_type import DbPhysioOutputType @@ -18,6 +19,7 @@ def insert_physio_file( modality: DbPhysioModality, output_type: DbPhysioOutputType, acquisition_time: datetime | None, + head_shape_file: DbMegCtfHeadShapeFile | None = None, ) -> DbPhysioFile: """ Insert a physiological file into the database. @@ -31,6 +33,7 @@ def insert_physio_file( output_type_id = output_type.id, acquisition_time = acquisition_time, inserted_by_user = getpass.getuser(), + head_shape_file_id = head_shape_file.id if head_shape_file is not None else None, ) env.db.add(file) diff --git a/python/loris_bids_importer/src/loris_bids_importer/main.py b/python/loris_bids_importer/src/loris_bids_importer/main.py index acf6dbafb..5ae87bd30 100644 --- a/python/loris_bids_importer/src/loris_bids_importer/main.py +++ b/python/loris_bids_importer/src/loris_bids_importer/main.py @@ -7,6 +7,7 @@ from lib.db.queries.session import try_get_session_with_cand_id_visit_label from lib.env import Env from lib.logging import log, log_error, log_error_exit, log_warning +from loris_bids_utils.meg.reader import BidsMegDataTypeReader from loris_bids_utils.mri.reader import BidsMriDataTypeReader from loris_bids_utils.reader import BidsDatasetReader, BidsDataTypeReader, BidsSessionReader @@ -20,6 +21,7 @@ from loris_bids_importer.eeg.main import Eeg from loris_bids_importer.env import BidsImportEnv from loris_bids_importer.events import import_bids_root_event_dict_file +from loris_bids_importer.meg.ctf import import_bids_meg_data_type from loris_bids_importer.mri.main import import_bids_mri_data_type from loris_bids_importer.print import print_bids_import_summary, print_bids_info from loris_bids_importer.validation.sessions import validate_bids_sessions @@ -177,11 +179,13 @@ def import_bids_data_type( match data_type: case BidsMriDataTypeReader(): import_bids_mri_data_type(env, import_env, session, data_type) + case BidsMegDataTypeReader(): + import_bids_meg_data_type(env, import_env, session, data_type) case BidsDataTypeReader(): - import_bids_eeg_data_type_files(env, import_env, args, session, data_type, dataset_tag_dict, legacy_db) + import_bids_eeg_data_type(env, import_env, args, session, data_type, dataset_tag_dict, legacy_db) -def import_bids_eeg_data_type_files( +def import_bids_eeg_data_type( env: Env, import_env: BidsImportEnv, args: Args, diff --git a/python/loris_bids_importer/src/loris_bids_importer/meg/ctf.py b/python/loris_bids_importer/src/loris_bids_importer/meg/ctf.py new file mode 100644 index 000000000..b95c80f39 --- /dev/null +++ b/python/loris_bids_importer/src/loris_bids_importer/meg/ctf.py @@ -0,0 +1,178 @@ +from pathlib import Path + +from lib.config import get_ephys_visualization_enabled_config +from lib.db.models.meg_ctf_head_shape_file import DbMegCtfHeadShapeFile +from lib.db.models.session import DbSession +from lib.db.queries.hed_schema_node import get_all_hed_schema_nodes +from lib.db.queries.meg_ctf_head_shape import try_get_meg_ctf_head_shape_file_with_path +from lib.db.queries.physio_file import try_get_physio_file_with_path +from lib.env import Env +from lib.logging import log, log_warning +from lib.physio.chunking import create_physio_channels_chunks +from lib.physio.ctf.head_shape import insert_meg_ctf_head_shape_file +from lib.physio.events import EventDictFileSource +from lib.physio.file import insert_physio_file +from lib.physio.parameters import insert_physio_file_parameter +from loris_bids_utils.info import BidsAcquisitionInfo +from loris_bids_utils.meg.ctf.acquisition import MegCtfAcquisition +from loris_bids_utils.meg.ctf.copy_ds import copy_bids_meg_ctf_ds +from loris_bids_utils.meg.reader import BidsMegDataTypeReader +from loris_utils.error import group_errors_tuple + +from loris_bids_importer.acquisitions import BidsImportFileResult, BidsImportFileStatus, import_bids_acquisitions +from loris_bids_importer.channels import insert_bids_channels_file +from loris_bids_importer.copy_files import copy_loris_bids_file, get_loris_bids_file_path +from loris_bids_importer.env import BidsImportEnv +from loris_bids_importer.events import insert_bids_event_dict_file, insert_bids_events_file +from loris_bids_importer.file_type import get_check_bids_imaging_file_type +from loris_bids_importer.physio import ( + get_check_bids_physio_file_hash, + get_check_bids_physio_modality, + get_check_bids_physio_output_type, +) + + +def import_bids_meg_data_type( + env: Env, + import_env: BidsImportEnv, + session: DbSession, + data_type: BidsMegDataTypeReader, +): + if data_type.head_shape_file is not None: + head_shape_file_path = get_loris_bids_file_path( + import_env, + session, + data_type.name, + data_type.head_shape_file.path, + ) + + head_shape_file = try_get_meg_ctf_head_shape_file_with_path(env.db, head_shape_file_path) + if head_shape_file is None: + head_shape_file = insert_meg_ctf_head_shape_file(env, data_type.head_shape_file, head_shape_file_path) + copy_loris_bids_file(import_env, data_type.head_shape_file.path, head_shape_file_path) + else: + head_shape_file = None + + import_bids_acquisitions( + env, + import_env, + session, + data_type.acquisitions, + lambda acquisition, bids_info: import_bids_meg_acquisition( + env, + import_env, + session, + acquisition, + bids_info, + head_shape_file, + ), + ) + + +def import_bids_meg_acquisition( + env: Env, + import_env: BidsImportEnv, + session: DbSession, + acquisition: MegCtfAcquisition, + bids_info: BidsAcquisitionInfo, + head_shape_file: DbMegCtfHeadShapeFile | None, +) -> BidsImportFileResult: + loris_file_path = get_loris_bids_file_path(import_env, session, bids_info.data_type, acquisition.ctf_path) + + loris_file = try_get_physio_file_with_path(env.db, loris_file_path) + if loris_file is not None: + return BidsImportFileResult(BidsImportFileStatus.IGNORE, loris_file_path) + + modality, output_type, file_type, file_hash = group_errors_tuple( + f"Error while checking database information for MEG acquisition '{bids_info.name}'.", + lambda: get_check_bids_physio_modality(env, bids_info.data_type), + # TODO: Use real output type after rebasing with #1443 + # Code: `importer.args.type or 'raw'` + lambda: get_check_bids_physio_output_type(env, 'raw'), + lambda: get_check_bids_imaging_file_type(env, 'ctf'), + lambda: get_check_bids_physio_file_hash(env, acquisition.ctf_path), + ) + + # The files to copy to LORIS, with the source path on the left and the LORIS path on the right. + files_to_copy: list[tuple[Path, Path]] = [] + + check_bids_meg_metadata_files(env, acquisition, bids_info) + + physio_file = insert_physio_file( + env, + session, + loris_file_path, + file_type, + modality, + output_type, + bids_info.scan_row.get_acquisition_time() if bids_info.scan_row is not None else None, + head_shape_file, + ) + + insert_physio_file_parameter(env, physio_file, 'physiological_file_blake2b_hash', file_hash) + for name, value in acquisition.sidecar_file.data.items(): + insert_physio_file_parameter(env, physio_file, name, value) + + if acquisition.events_file is not None: + hed_union = get_all_hed_schema_nodes(env.db) + + loris_events_file_path = get_loris_bids_file_path( + import_env, session, bids_info.data_type, acquisition.events_file.path + ) + + insert_bids_events_file(env, physio_file, acquisition.events_file, loris_events_file_path, {}, {}, hed_union) + files_to_copy.append((acquisition.events_file.path, loris_events_file_path)) + if acquisition.events_file.dictionary is not None: + loris_event_dict_file_path = get_loris_bids_file_path( + import_env, session, bids_info.data_type, acquisition.events_file.dictionary.path + ) + + insert_bids_event_dict_file( + env, + EventDictFileSource.from_file(physio_file), + acquisition.events_file.dictionary, + loris_event_dict_file_path, + ) + + files_to_copy.append((acquisition.events_file.dictionary.path, loris_event_dict_file_path)) + + if acquisition.channels_file is not None: + insert_bids_channels_file(env, import_env, physio_file, session, bids_info, acquisition.channels_file) + loris_channels_file_path = get_loris_bids_file_path( + import_env, session, bids_info.data_type, acquisition.channels_file.path + ) + files_to_copy.append((acquisition.channels_file.path, loris_channels_file_path)) + + if import_env.loris_bids_path is not None: + copy_bids_meg_ctf_ds(acquisition.ctf_path, import_env.data_dir_path / loris_file_path) + + for source_path, destination_path in files_to_copy: + copy_loris_bids_file(import_env, source_path, destination_path) + + env.db.commit() + + log(env, f"MEG file successfully imported with ID: {physio_file.id}.") + + if get_ephys_visualization_enabled_config(env): + log(env, "Creating visualization chunks...") + create_physio_channels_chunks(env, physio_file) + + env.db.commit() + + return BidsImportFileResult(BidsImportFileStatus.SUCCESS, loris_file_path) + + +def check_bids_meg_metadata_files(env: Env, acquisition: MegCtfAcquisition, bids_info: BidsAcquisitionInfo): + """ + Check for the presence of BIDS metadata files for the BIDS MEG acquisition and warn the user if + that is not the case. + """ + + if acquisition.channels_file is None: + log_warning(env, f"No channels file found for acquisition '{bids_info.name}'.") + + if acquisition.events_file is None: + log_warning(env, f"No events file found for acquisition '{bids_info.name}'.") + + if acquisition.events_file is not None and acquisition.events_file.dictionary is not None: + log_warning(env, f"No events dictionary file found for acquisition '{bids_info.name}'.") diff --git a/python/loris_bids_importer/src/loris_bids_importer/physio.py b/python/loris_bids_importer/src/loris_bids_importer/physio.py index c85c6f153..d80b13c95 100644 --- a/python/loris_bids_importer/src/loris_bids_importer/physio.py +++ b/python/loris_bids_importer/src/loris_bids_importer/physio.py @@ -5,7 +5,7 @@ from lib.db.queries.physio import try_get_physio_modality_with_name, try_get_physio_output_type_with_name from lib.db.queries.physio_file import try_get_physio_file_with_hash from lib.env import Env -from loris_utils.crypto import compute_file_blake2b_hash +from loris_utils.crypto import compute_directory_blake2b_hash, compute_file_blake2b_hash def get_check_bids_physio_modality(env: Env, data_type_name: str) -> DbPhysioModality: @@ -40,7 +40,10 @@ def get_check_bids_physio_file_hash(env: Env, file_path: Path) -> str: registered in the database. """ - file_hash = compute_file_blake2b_hash(file_path) + if file_path.is_dir(): + file_hash = compute_directory_blake2b_hash(file_path) + else: + file_hash = compute_file_blake2b_hash(file_path) file = try_get_physio_file_with_hash(env.db, file_hash) if file is not None: diff --git a/python/loris_bids_utils/pyproject.toml b/python/loris_bids_utils/pyproject.toml index a2528c8c0..754d93ca5 100644 --- a/python/loris_bids_utils/pyproject.toml +++ b/python/loris_bids_utils/pyproject.toml @@ -10,6 +10,10 @@ dependencies = [ "pybids", ] +[project.scripts] +copy-bids-ctf-ds = "loris_bids_utils.scripts.copy_bids_ctf_ds:main" +copy-ctf-ds = "loris_bids_utils.scripts.copy_ctf_ds:main" + [build-system] requires = ["hatchling"] build-backend = "hatchling.build" diff --git a/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/acquisition.py b/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/acquisition.py new file mode 100644 index 000000000..971ad40a9 --- /dev/null +++ b/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/acquisition.py @@ -0,0 +1,37 @@ + +import re +from pathlib import Path + +from loris_utils.path import remove_path_extension + +from loris_bids_utils.eeg.channels import BidsEegChannelsTsvFile +from loris_bids_utils.files.events import BidsEventsTsvFile +from loris_bids_utils.meg.ctf.head_shape import MegCtfHeadShapeFile +from loris_bids_utils.meg.sidecar import BidsMegSidecarJsonFile + + +class MegCtfAcquisition: + ctf_path: Path + sidecar_file: BidsMegSidecarJsonFile + channels_file: BidsEegChannelsTsvFile | None + events_file: BidsEventsTsvFile | None + head_shape_file: MegCtfHeadShapeFile | None + + def __init__(self, ctf_path: Path, head_shape_file: MegCtfHeadShapeFile | None): + self.ctf_path = ctf_path + + path = remove_path_extension(ctf_path) + + sidecar_path = path.with_suffix('.json') + if not sidecar_path.exists(): + raise Exception("No MEG JSON sidecar file.") + + self.sidecar_file = BidsMegSidecarJsonFile(sidecar_path) + + channels_path = path.parent / re.sub(r'_meg$', '_channels.tsv', path.name) + self.channels_file = BidsEegChannelsTsvFile(channels_path) if channels_path.exists() else None + + events_path = path.parent / re.sub(r'_meg$', '_events.tsv', path.name) + self.events_file = BidsEventsTsvFile(events_path) if events_path.exists() else None + + self.head_shape_file = head_shape_file diff --git a/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/copy_ds.py b/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/copy_ds.py new file mode 100644 index 000000000..1252f4f27 --- /dev/null +++ b/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/copy_ds.py @@ -0,0 +1,450 @@ +""" +Copy CTF MEG .ds directories while updating their internal BIDS-derived metadata and CTF references. + +This code is a Python port of the following Matlab code: +https://github.com/Moo-Marc/CtfMegBids/ +""" + +import re +import shutil +import warnings +from collections.abc import Iterable +from pathlib import Path + +from loris_bids_utils.path import build_bids_modality_path, parse_bids_entities + +KNOWN_DATASET_FILES = { + "processing.cfg", + "processing.cfg.bak", + "ClassFile.cls", + "ClassFile.cls.bak", + "MarkerFile.mrk", + "MarkerFile.mrk.bak", + "default.de", + "default.de.bak", + "bad.segments", + "bad.segments.bak", + "BadChannels", + "BadChannels.bak", + "ChannelGroupSet.cfg", + "ChannelGroupSet.cfg.bak", + "VirtualChannels", + "VirtualChannels.bak", + "DigTrigChannelInfo.txt", + "DigTrigChannelInfo.txt.bak", + "params.dsc", + "params.dsc.bak", + "hz.ds", + ".lock", +} + + +def copy_meg_ctf_ds(source_path: Path, destination_path: Path) -> list[Path]: + """Copy a CTF ``.ds`` dataset and update internal CTF references. + + This is a Python port of the non-anonymizing behavior in + ``Bids_ctf_rename_ds.m``, except it copies the source dataset instead of + moving it. It renames files whose names start with the old dataset stem in + the copy, removes empty ``.bak`` and ``.eeg`` files from the copy, updates + ``.hist`` references, and rewrites the dataset path line in ``.cls`` and + ``.mrk`` files. + + Args: + source_path: Path to the source CTF dataset directory. + destination_path: Exact path of the destination CTF dataset directory. + + Returns: + Unknown files found after the rename, using paths relative to the new + dataset directory. + """ + + return _copy_meg_ctf_ds( + source_path, + destination_path, + Path(destination_path.name), + ) + + +def _copy_meg_ctf_ds(source_path: Path, destination_path: Path, dataset_reference_path: Path) -> list[Path]: + _validate_dataset_paths(source_path, destination_path) + + orig_name = _strip_ds_suffix(source_path.name) + new_name = _strip_ds_suffix(destination_path.name) + + _copy_dataset(source_path, destination_path) + _rename_dataset_files(destination_path, orig_name, new_name) + unknown_files = _find_unknown_files(destination_path, new_name) + + if unknown_files: + warnings.warn(f"Unknown files in {destination_path}", stacklevel=2) + + _delete_empty_files(destination_path.glob("*.bak")) + _delete_empty_files(destination_path.glob("*.eeg")) + _replace_text_in_first_file( + destination_path, "*.hist", orig_name, new_name, "history" + ) + _replace_dataset_path_line( + destination_path, "*.cls", dataset_reference_path, "class" + ) + _replace_dataset_path_line( + destination_path, "*.mrk", dataset_reference_path, "marker" + ) + _warn_if_old_name_remains(destination_path, orig_name) + + return unknown_files + + +def copy_bids_meg_ctf_ds(source_path: Path, destination_path: Path) -> list[Path]: + """Copy a BIDS CTF ``.ds`` dataset and update BIDS-derived metadata.""" + + _validate_dataset_paths(source_path, destination_path) + + orig_name = _strip_ds_suffix(source_path.name) + new_name = _strip_ds_suffix(destination_path.name) + relative_dataset_path = _relative_dataset_path(destination_path.name) + subject_rename = _get_subject_rename(orig_name, new_name) + if subject_rename is not None: + _encode_res4_subject(subject_rename[1]) + + unknown_files = _copy_meg_ctf_ds( + source_path, + destination_path, + relative_dataset_path, + ) + if subject_rename is not None: + _replace_subject_metadata(destination_path, subject_rename[1]) + + return unknown_files + + +def _validate_dataset_paths(source_path: Path, destination_path: Path): + if not source_path.is_dir(): + raise FileNotFoundError(f"Dataset not found: {source_path}") + if destination_path.suffix != ".ds": + raise ValueError( + f"Destination must be an explicit .ds dataset path: {destination_path}" + ) + if source_path.resolve() == destination_path.resolve(): + raise ValueError("Destination must differ from the source dataset path.") + + +def _strip_ds_suffix(name: str) -> str: + name_path = Path(name) + if name_path.suffix == ".ds": + return name_path.stem + return name_path.name + + +def _relative_dataset_path(new_ds_name: str) -> Path: + if new_ds_name == 'hz.ds': + return Path(new_ds_name) + + if not new_ds_name.endswith('_meg.ds'): + raise ValueError( + f"New dataset name should follow BIDS specification: {new_ds_name}" + ) + entities = parse_bids_entities(new_ds_name) + if 'sub' not in entities or 'task' not in entities: + raise ValueError( + f"New dataset name should follow BIDS specification: {new_ds_name}" + ) + return build_bids_modality_path( + entities['sub'], entities.get('ses'), 'meg', new_ds_name + ) + + +def _copy_dataset(original_ds_path: Path, new_ds_path: Path): + if original_ds_path == new_ds_path: + raise ValueError("Destination must differ from the source dataset path.") + + if new_ds_path.exists(): + if not new_ds_path.is_dir(): + raise FileExistsError( + f"Destination exists and is not a directory: {new_ds_path}" + ) + for child_path in original_ds_path.iterdir(): + destination_child_path = new_ds_path / child_path.name + if child_path.is_dir(): + shutil.copytree( + child_path, destination_child_path, dirs_exist_ok=True + ) + else: + shutil.copy2(child_path, destination_child_path) + else: + new_ds_path.parent.mkdir(parents=True, exist_ok=True) + shutil.copytree(original_ds_path, new_ds_path) + + +def _rename_dataset_files(ds_path: Path, orig_name: str, new_name: str): + if orig_name == new_name: + return + + for file_path in sorted(ds_path.glob(f"{orig_name}*")): + renamed_path = file_path.with_name( + file_path.name.replace(orig_name, new_name) + ) + shutil.move(str(file_path), str(renamed_path)) + + +def _find_unknown_files(ds_path: Path, new_name: str) -> list[Path]: + expected_names = { + candidate_path.name for candidate_path in ds_path.glob(f"{new_name}*") + } + expected_names.update(KNOWN_DATASET_FILES) + + unknown = [ + candidate_path.relative_to(ds_path) + for candidate_path in ds_path.iterdir() + if candidate_path.name not in expected_names + ] + return sorted(unknown) + + +def _delete_empty_files(files: Iterable[Path]): + for file_path in files: + if file_path.is_file() and file_path.stat().st_size == 0: + file_path.unlink() + + +def _replace_text_in_first_file(ds_path: Path, pattern: str, old: str, new: str, file_kind: str): + matches = sorted(ds_path.glob(pattern)) + if len(matches) > 1: + warnings.warn(f"Multiple {file_kind} files found in {ds_path}", stacklevel=2) + matches = matches[:1] + if not matches: + return + + file_path = matches[0] + content = _read_text_lossless(file_path) + file_path.write_text(content.replace(old, new), encoding="latin-1") + + +def _replace_dataset_path_line(ds_path: Path, pattern: str, relative_dataset_path: Path, file_kind: str): + matches = sorted(ds_path.glob(pattern)) + if len(matches) > 1: + warnings.warn(f"Multiple {file_kind} files found in {ds_path}", stacklevel=2) + matches = matches[:1] + if not matches: + return + + file_path = matches[0] + content = _read_text_lossless(file_path) + line_breaks = [match.start() for match in re.finditer(r"\n", content)] + if len(line_breaks) < 2: + warnings.warn( + f"Could not update dataset path in {file_path}: expected at least 2 lines", + stacklevel=2, + ) + return + + start = line_breaks[0] + 1 + end = line_breaks[1] + old_path = content[start:end] + updated = content.replace(old_path, str(relative_dataset_path)) + file_path.write_text(updated, encoding="latin-1") + + +def _get_subject_rename(orig_name: str, new_name: str) -> tuple[str, str] | None: + orig_name_with_ext = f"{orig_name}.ds" + new_name_with_ext = f"{new_name}.ds" + if not orig_name.endswith("_meg") or not new_name.endswith("_meg"): + return None + try: + orig_info = parse_bids_entities(orig_name_with_ext) + new_info = parse_bids_entities(new_name_with_ext) + except ValueError: + return None + + old_subject_label = orig_info.get("sub") + new_subject_label = new_info.get("sub") + if ( + old_subject_label is None + or new_subject_label is None + or old_subject_label == new_subject_label + ): + return None + + return f"sub-{old_subject_label}", f"sub-{new_subject_label}" + + +def _replace_subject_metadata(ds_path: Path, new_subject: str): + _replace_res4_subject(ds_path, new_subject) + _replace_infods_patient_id(ds_path, new_subject) + + +def _encode_res4_subject(subject: str) -> bytes: + try: + encoded = subject.encode("latin-1") + except UnicodeEncodeError as exc: + raise ValueError(f"CTF res4 subject id must be Latin-1: {subject}") from exc + if len(encoded) > 32: + raise ValueError(f"CTF res4 subject id is limited to 32 bytes: {subject}") + return encoded + (b"\x00" * (32 - len(encoded))) + + +def _replace_res4_subject(ds_path: Path, subject: str): + padded = _encode_res4_subject(subject) + + for file_path in sorted(ds_path.glob("*.res4")): + with file_path.open("r+b") as fid: + fid.seek(1712) + fid.write(padded) + + +def _replace_infods_patient_id(ds_path: Path, new_subject: str): + replacements = {"_PATIENT_ID": new_subject} + for file_path in sorted(ds_path.glob("*.infods")): + updated = _replace_cpersist_strings(file_path.read_bytes(), replacements) + file_path.write_bytes(updated) + + +def _replace_cpersist_strings(content: bytes, replacements: dict[str, str]) -> bytes: + start_value = int.from_bytes(b"WS1_", byteorder="big", signed=True) + end_tag = "EndOfParameters" + out = bytearray() + pos = 0 + + while pos + 4 <= len(content): + tag_start = pos + tag_length = int.from_bytes(content[pos : pos + 4], "big", signed=True) + pos += 4 + + if tag_length == start_value: + out.extend(content[tag_start:pos]) + continue + if tag_length <= 0 or pos + tag_length > len(content): + out.extend(content[tag_start:]) + return bytes(out) + + tag_name = content[pos : pos + tag_length].decode("latin-1") + pos += tag_length + out.extend(content[tag_start:pos]) + + if tag_name == end_tag: + continue + if pos + 4 > len(content): + out.extend(content[pos:]) + return bytes(out) + + tag_type = int.from_bytes(content[pos : pos + 4], "big", signed=True) + out.extend(content[pos : pos + 4]) + pos += 4 + + if tag_type == 10: + if pos + 4 > len(content): + out.extend(content[pos:]) + return bytes(out) + string_length = int.from_bytes(content[pos : pos + 4], "big", signed=True) + value_start = pos + 4 + value_end = value_start + string_length + if string_length < 0 or value_end > len(content): + out.extend(content[pos:]) + return bytes(out) + + value = content[value_start:value_end].decode("latin-1") + if tag_name in replacements: + value = replacements[tag_name] + encoded = value.encode("latin-1") + out.extend(len(encoded).to_bytes(4, "big", signed=True)) + out.extend(encoded) + pos = value_end + continue + + value_end = _cpersist_value_end(content, pos, tag_type, tag_name) + if value_end is None: + out.extend(content[pos:]) + return bytes(out) + out.extend(content[pos:value_end]) + pos = value_end + + out.extend(content[pos:]) + return bytes(out) + + +def _cpersist_value_end(content: bytes, pos: int, tag_type: int, tag_name: str) -> int | None: + fixed_widths = { + 2: 0, + 4: 8, + 5: 4, + 6: 2, + 7: 2, + 8: 1, + 9: 32, + 14: 4, + 15: 4, + 16: 4, + 17: 4, + } + if tag_type == 1: + return pos + 4 if tag_name == "DatasetFiles" else pos + if tag_type in fixed_widths: + return pos + fixed_widths[tag_type] + if tag_type == 3: + if pos + 4 > len(content): + return None + byte_count = int.from_bytes(content[pos : pos + 4], "big", signed=True) + return pos + 4 + byte_count + if tag_type == 11: + if pos + 4 > len(content): + return None + count = int.from_bytes(content[pos : pos + 4], "big", signed=True) + value_pos = pos + 4 + for _ in range(count): + if value_pos + 4 > len(content): + return None + string_length = int.from_bytes( + content[value_pos : value_pos + 4], + "big", + signed=True, + ) + value_pos += 4 + string_length + return value_pos + if tag_type == 12: + if pos + 4 > len(content): + return None + count = int.from_bytes(content[pos : pos + 4], "big", signed=True) + return pos + 4 + (32 * count) + if tag_type == 13: + if pos + 4 > len(content): + return None + count = int.from_bytes(content[pos : pos + 4], "big", signed=True) + return pos + 4 + (4 * count) + return None + + +def _warn_if_old_name_remains(ds_path: Path, orig_name: str): + if orig_name == _strip_ds_suffix(ds_path.name): + return + + old_name = orig_name.encode("latin-1") + matching_paths = [ + candidate_path.relative_to(ds_path) + for candidate_path in ds_path.iterdir() + if candidate_path.is_file() and _file_contains(candidate_path, old_name) + ] + if matching_paths: + warnings.warn( + "Old dataset name still appears in copied files: " + + ", ".join(str(match_path) for match_path in sorted(matching_paths)), + stacklevel=2, + ) + + +def _file_contains(file_path: Path, needle: bytes) -> bool: + chunk_size = 1024 * 1024 + overlap = max(len(needle) - 1, 0) + previous = b"" + + with file_path.open("rb") as fid: + while True: + chunk = fid.read(chunk_size) + if not chunk: + return False + data = previous + chunk + if needle in data: + return True + previous = data[-overlap:] if overlap else b"" + + +def _read_text_lossless(file_path: Path) -> str: + return file_path.read_text(encoding="latin-1") diff --git a/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/head_shape.py b/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/head_shape.py new file mode 100644 index 000000000..9966b8104 --- /dev/null +++ b/python/loris_bids_utils/src/loris_bids_utils/meg/ctf/head_shape.py @@ -0,0 +1,111 @@ +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +from loris_utils.iter import get_first + +VecF = np.ndarray[tuple[int], np.dtype[np.float64]] + + +@dataclass +class MegCtfHeadShapePoint: + """ + A point in a MEG CTF `headshape.pos` file. + """ + + x: float + y: float + z: float + + def scale(self, factor: float) -> 'MegCtfHeadShapePoint': + """ + Scale the point coordinates by a factor. Notably useful to convert the head shape point from + one unit to another. + """ + + return MegCtfHeadShapePoint( + x = self.x * factor, + y = self.y * factor, + z = self.z * factor, + ) + + def to_numpy(self) -> VecF: + """ + Convert the point to a numpy array. + """ + + return np.array([self.x, self.y, self.z], dtype=np.float64) + + +@dataclass +class MegCtfHeadShapeFile: + """ + A MEG CTF `headshape.pos` file. + """ + + path: Path + """ + The path of this head shape file. + """ + + points: dict[str, MegCtfHeadShapePoint] + """ + The points of this head shape file. + """ + + @staticmethod + def read(path: Path) -> 'MegCtfHeadShapeFile': + """ + Read and parse a MEG CTF head shape file. + """ + + with path.open() as file: + lines = file.readlines() + + points: dict[str, MegCtfHeadShapePoint] = {} + # The first line simply gives the number of points. + for line in lines[1:]: + parts = line.split() + # The first column contains the sensor name or index. + # The second column may or may not be empty. + # The last three columns contain the point coordimates. + points[parts[0]] = MegCtfHeadShapePoint(float(parts[-3]), float(parts[-2]), float(parts[-1])) + + return MegCtfHeadShapeFile( + path = path, + points = points, + ) + + def scale(self, factor: float) -> 'MegCtfHeadShapeFile': + """ + Create a new head shape file with all points scaled by a factor. + """ + + return MegCtfHeadShapeFile( + path = self.path, + points = {name: point.scale(factor) for name, point in self.points.items()}, + ) + + @property + def nasion(self) -> MegCtfHeadShapePoint | None: + """ + Get the nasion fiducial point. + """ + + return get_first(self.points, ['NAS', 'Nasion', 'nasion']) + + @property + def lpa(self) -> MegCtfHeadShapePoint | None: + """ + Get the LPA fiducial point. + """ + + return get_first(self.points, ['LPA', 'left']) + + @property + def rpa(self) -> MegCtfHeadShapePoint | None: + """ + Get the RPA fiducial point. + """ + + return get_first(self.points, ['RPA', 'right']) diff --git a/python/loris_bids_utils/src/loris_bids_utils/meg/reader.py b/python/loris_bids_utils/src/loris_bids_utils/meg/reader.py new file mode 100644 index 000000000..2fd212d7c --- /dev/null +++ b/python/loris_bids_utils/src/loris_bids_utils/meg/reader.py @@ -0,0 +1,75 @@ +import re +from collections.abc import Iterator +from dataclasses import dataclass +from functools import cached_property +from pathlib import Path + +from loris_bids_utils.info import BidsAcquisitionInfo +from loris_bids_utils.meg.ctf.acquisition import MegCtfAcquisition +from loris_bids_utils.meg.ctf.head_shape import MegCtfHeadShapeFile +from loris_bids_utils.reader import BidsDataTypeReader +from loris_bids_utils.utils import get_pybids_file_path, try_get_pybids_value + + +@dataclass +class BidsMegDataTypeReader(BidsDataTypeReader): + path: Path + + @cached_property + def acquisitions(self) -> list[tuple[MegCtfAcquisition, BidsAcquisitionInfo]]: + """ + The MEG acquisitions found in the MEG data type. + """ + + acquisitions: list[tuple[MegCtfAcquisition, BidsAcquisitionInfo]] = [] + for ctf_name in find_dir_meg_acquisition_names(self.path): + scan_row = self.session.scans_file.get_row(self.path / ctf_name) \ + if self.session.scans_file is not None else None + + acquisition = MegCtfAcquisition(self.path / ctf_name, self.head_shape_file) + + info = BidsAcquisitionInfo( + subject = self.session.subject.label, + participant_row = self.session.subject.participant_row, + session = self.session.label, + scans_file = self.session.scans_file, + data_type = self.name, + scan_row = scan_row, + name = ctf_name, + suffix = 'meg', + ) + + acquisitions.append((acquisition, info)) + + return acquisitions + + @cached_property + def head_shape_file(self) -> MegCtfHeadShapeFile | None: + """ + The MEG CTF file of this acquisition if it exists. + """ + + head_shape_file = try_get_pybids_value( + self.session.subject.dataset.layout, + subject=self.session.subject.label, + session=self.session.label, + datatype=self.name, + suffix='headshape', + extension='.pos', + ) + + if head_shape_file is None: + return None + + return MegCtfHeadShapeFile.read(get_pybids_file_path(head_shape_file)) + + +def find_dir_meg_acquisition_names(dir_path: Path) -> Iterator[str]: + """ + Iterate over the Path objects of the NIfTI files found in a directory. + """ + + for item_path in dir_path.iterdir(): + name_match = re.search(r'.+_meg\.ds$', item_path.name) + if name_match is not None: + yield name_match.group(0) diff --git a/python/loris_bids_utils/src/loris_bids_utils/meg/sidecar.py b/python/loris_bids_utils/src/loris_bids_utils/meg/sidecar.py new file mode 100644 index 000000000..706054855 --- /dev/null +++ b/python/loris_bids_utils/src/loris_bids_utils/meg/sidecar.py @@ -0,0 +1,12 @@ +from loris_bids_utils.json import BidsJsonFile + + +class BidsMegSidecarJsonFile(BidsJsonFile): + """ + Class representing a BIDS EEG or iEEG sidecar JSON file. + + Documentation: + - https://bids-specification.readthedocs.io/en/stable/modality-specific-files/magnetoencephalography.html#sidecar-json-_megjson + """ + + pass diff --git a/python/loris_bids_utils/src/loris_bids_utils/reader.py b/python/loris_bids_utils/src/loris_bids_utils/reader.py index 3d010b3bb..14c6f05d6 100644 --- a/python/loris_bids_utils/src/loris_bids_utils/reader.py +++ b/python/loris_bids_utils/src/loris_bids_utils/reader.py @@ -16,6 +16,7 @@ # Circular imports if TYPE_CHECKING: + from loris_bids_utils.meg.reader import BidsMegDataTypeReader from loris_bids_utils.mri.reader import BidsMriDataTypeReader PYBIDS_IGNORE = ['.git', 'code/', 'log/', 'sourcedata/'] @@ -297,13 +298,38 @@ def eeg_data_types(self) -> list['BidsDataTypeReader']: ) ] + @cached_property + def meg_data_types(self) -> list['BidsMegDataTypeReader']: + """ + Get the MEG data type directory readers of this session. + """ + + from loris_bids_utils.meg.reader import BidsMegDataTypeReader + + return [ + BidsMegDataTypeReader( + session=self, + name=data_type, # type: ignore + path=( + self.subject.dataset.path + / f'sub-{self.subject.label}' + / (f'ses-{self.label}' if self.label is not None else '') + / data_type # type: ignore + ), + ) for data_type in self.subject.dataset.layout.get_datatypes( # type: ignore + subject=self.subject.label, + session=self.label, + datatype=['meg'], + ) + ] + @cached_property def data_types(self) -> Sequence['BidsDataTypeReader']: """ Get all the data type directory readers of this session. """ - return self.eeg_data_types + self.mri_data_types + return self.eeg_data_types + self.meg_data_types + self.mri_data_types @cached_property def info(self) -> BidsSessionInfo: diff --git a/python/loris_bids_utils/src/loris_bids_utils/scripts/copy_bids_ctf_ds.py b/python/loris_bids_utils/src/loris_bids_utils/scripts/copy_bids_ctf_ds.py new file mode 100755 index 000000000..7d9065ec5 --- /dev/null +++ b/python/loris_bids_utils/src/loris_bids_utils/scripts/copy_bids_ctf_ds.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python + +import argparse +from pathlib import Path + +from loris_bids_utils.meg.ctf.copy_ds import copy_bids_meg_ctf_ds + + +def main() -> None: + parser = argparse.ArgumentParser( + description=( + "Copy and rename a BIDS MEG CTF .ds dataset, updating its " + "BIDS-derived metadata and internal CTF references in the process." + ) + ) + + parser.add_argument( + 'source', + type=Path, + help="Source BIDS CTF .ds path", + ) + + parser.add_argument( + 'destination', + type=Path, + help="Destination .ds path or parent directory path", + ) + + args = parser.parse_args() + + source_path: Path = args.source + destination_path: Path = args.destination + + if destination_path.suffix != ".ds": + destination_path = destination_path / source_path.name + + unknown_paths = copy_bids_meg_ctf_ds(source_path, destination_path) + + print(f"Created: {destination_path}") + + if unknown_paths: + print("Unknown files:") + for unknown_path in unknown_paths: + print(f" {unknown_path}") + + +if __name__ == '__main__': + main() diff --git a/python/loris_bids_utils/src/loris_bids_utils/scripts/copy_ctf_ds.py b/python/loris_bids_utils/src/loris_bids_utils/scripts/copy_ctf_ds.py new file mode 100755 index 000000000..6299ac0a1 --- /dev/null +++ b/python/loris_bids_utils/src/loris_bids_utils/scripts/copy_ctf_ds.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python + +import argparse +from pathlib import Path + +from loris_bids_utils.meg.ctf.copy_ds import copy_meg_ctf_ds + + +def main() -> None: + parser = argparse.ArgumentParser( + description=( + "Copy and rename a MEG CTF .ds directory, updating its internal " + "CTF references in the process." + ) + ) + + parser.add_argument( + 'source', + type=Path, + help="Source CTF .ds path", + ) + + parser.add_argument( + 'destination', + type=Path, + help="Destination .ds path or parent directory path", + ) + + args = parser.parse_args() + + source_path: Path = args.source + destination_path: Path = args.destination + + if destination_path.suffix != ".ds": + destination_path = destination_path / source_path.name + + unknown_paths = copy_meg_ctf_ds(source_path, destination_path) + + print(f"Created: {destination_path}") + + if unknown_paths: + print("Unknown files:") + for unknown_path in unknown_paths: + print(f" {unknown_path}") + + +if __name__ == "__main__": + main() diff --git a/python/loris_bids_utils/src/loris_bids_utils/tsv.py b/python/loris_bids_utils/src/loris_bids_utils/tsv.py index a18618e0a..aaa4c7224 100644 --- a/python/loris_bids_utils/src/loris_bids_utils/tsv.py +++ b/python/loris_bids_utils/src/loris_bids_utils/tsv.py @@ -3,6 +3,9 @@ from typing import Generic, TypeVar from loris_utils.parse import nullify_empty_string +from loris_utils.path import replace_path_extension + +from loris_bids_utils.json import BidsJsonFile class BidsTsvRow: @@ -27,12 +30,19 @@ class BidsTsvFile(Generic[T]): """ path: Path + dictionary: BidsJsonFile | None rows: list[T] def __init__(self, model: type[T], path: Path): self.path = path self.rows = [] + dictionary_path = replace_path_extension(self.path, 'json') + if dictionary_path.exists(): + self.dictionary = BidsJsonFile(dictionary_path) + else: + self.dictionary = None + # The 'utf-8-sig' encoding is used to support some datasets where metadata files may contain # a byte-order mark (BOM). with open(self.path, encoding='utf-8-sig') as file: diff --git a/python/loris_bids_utils/src/loris_bids_utils/utils.py b/python/loris_bids_utils/src/loris_bids_utils/utils.py index 9d633f697..d3aaed347 100644 --- a/python/loris_bids_utils/src/loris_bids_utils/utils.py +++ b/python/loris_bids_utils/src/loris_bids_utils/utils.py @@ -13,7 +13,7 @@ def try_get_pybids_value(layout: BIDSLayout, **args: Any) -> Any | None: values are found. """ - match layout.get(args): # type: ignore + match layout.get(**args): # type: ignore case []: return None case [value]: # type: ignore diff --git a/python/loris_utils/src/loris_utils/crypto.py b/python/loris_utils/src/loris_utils/crypto.py index 84c25910c..64d1aff4d 100644 --- a/python/loris_utils/src/loris_utils/crypto.py +++ b/python/loris_utils/src/loris_utils/crypto.py @@ -1,19 +1,55 @@ import hashlib +from hashlib import blake2b from pathlib import Path -def compute_file_blake2b_hash(file_path: Path | str) -> str: +def compute_file_blake2b_hash(file_path: Path) -> str: """ Compute the BLAKE2b hash of a file. """ + hash = blake2b() + update_file_blake2b_hash(Path(file_path), hash) + return hash.hexdigest() + + +def compute_directory_blake2b_hash(dir_path: Path) -> str: + """ + Compute the BLAKE2b hash of a directory. + """ + + hash = blake2b() + update_directory_blake2b_hash(dir_path, hash) + return hash.hexdigest() + + +def update_file_blake2b_hash(file_path: Path, hash: blake2b): + """ + Update a BLAKE2b hash with the contents of a file. + """ + # Since the file given to this function may be large, we read it in chunks to avoid running # out of memory. - hash = hashlib.blake2b() with open(file_path, 'rb') as file: while chunk := file.read(1048576): hash.update(chunk) - return hash.hexdigest() + + +def update_directory_blake2b_hash(dir_path: Path, hash: blake2b): + """ + Update a BLAKE2b hash with the contents of a directory. + """ + + # The paths are sorted to ensure the hash is deterministic regardless of iteration order. + for path in sorted(dir_path.iterdir()): + # The file name is included in the hash to ensure the directory structure is reflected in + # the hash. + hash.update(path.name.encode()) + # Symlinks are currently not included in the hash. + if path.is_file(): + update_file_blake2b_hash(path, hash) + elif path.is_dir(): + update_directory_blake2b_hash(path, hash) def compute_file_md5_hash(file_path: Path | str) -> str: diff --git a/python/loris_utils/src/loris_utils/iter.py b/python/loris_utils/src/loris_utils/iter.py index d5d8535c8..4f596e35b 100644 --- a/python/loris_utils/src/loris_utils/iter.py +++ b/python/loris_utils/src/loris_utils/iter.py @@ -1,4 +1,4 @@ -from collections.abc import Callable, Iterable, Iterator, Sized +from collections.abc import Callable, Iterable, Iterator, Mapping, Sized from typing import TypeVar T = TypeVar('T') @@ -89,9 +89,27 @@ def map_non_none(value: T | None, function: Callable[[T], U]) -> U | None: def filter_non_none(iterable: Iterable[T | None]) -> Iterator[T]: """ - Filter the `None` elements out of an iterator. + Filter an iterator by removing its `None` elements. """ for element in iterable: if element is not None: yield element + + +K = TypeVar('K') +V = TypeVar('V') + + +def get_first(mapping: Mapping[K, V], keys: Iterable[K]) -> V | None: + """ + Get the first non-`None` value of a mapping using a list of keys, or return `None` if no + value is found. + """ + + for key in keys: + value = mapping.get(key) + if value is not None: + return value + + return None