Two different CTF examples.

examples/ctf/README.md → examples/ctf_with_headshape_pos/README.md

@@ -1,6 +1,8 @@
 # CTF Example
 
-In this example we:
+In this example we use headshape points to coregister the MEG and sMRI data.
+
+Full pipeline:
 
 - `1_preprocess.py`: Preprocess the sensor-level data. This includes standard signal processing such as downsampling and filtering as well as artefact removal.
 - `2_coregister.py`: Coregister the MEG and sMRI data and calculate the forward model.

examples/ctf_with_smri_fid/1_preprocess.py

@@ -0,0 +1,29 @@
+"""Preprocessing.
+
+"""
+
+# Authors: Chetan Gohil <chetan.gohil@psych.ox.ac.uk>
+
+from osl import preprocessing
+
+# Settings
+config = """
+    preproc:
+    - filter: {l_freq: 1, h_freq: 125, method: iir, iir_params: {order: 5, ftype: butter}}
+    - notch_filter: {freqs: 50 100}
+    - resample: {sfreq: 250}
+"""
+
+# Create a list of paths to files to preprocess
+inputs = ["data/raw/mg04938_BrainampDBS_20170504_01_raw.fif"]
+
+# Directory to save output to
+outdir = "data/preproc"
+
+# Do preprocessing
+preprocessing.run_proc_batch(
+    config,
+    inputs,
+    outdir=outdir,
+    overwrite=True,
+)

examples/ctf_with_smri_fid/2_coregister.py

@@ -0,0 +1,42 @@
+"""Coregistration.
+
+"""
+
+# Authors: Chetan Gohil <chetan.gohil@psych.ox.ac.uk>
+
+import numpy as np
+
+from osl import source_recon, utils
+
+
+config = """
+    source_recon:   
+    - extract_polhemus_from_info: {}
+    - save_mni_fiducials:
+        filepath: data/fiducials/{subject}_smri_fid.txt
+    - compute_surfaces:
+        include_nose: False    
+    - coregister:
+        use_nose: False
+        use_headshape: False
+"""
+
+# List of subject IDs
+subjects = ["LN_VTA2"]
+
+# Lists for input files
+preproc_files = ["data/preproc/mg04938_BrainampDBS_20170504_01_preproc_raw.fif"]
+smri_files = ["data/smri/LN_VTA2.nii"]
+
+# Output directory
+coreg_dir = "data/coreg"
+
+# Do coregistration
+source_recon.run_src_batch(
+    config,
+    src_dir=coreg_dir,
+    subjects=subjects,
+    preproc_files=preproc_files,
+    smri_files=smri_files,
+    extra_funcs=[save_mni_fids],
+)

examples/ctf_with_smri_fid/3_source_reconstruct.py

@@ -0,0 +1,49 @@
+"""Source reconstruction using a LCMV beamformer and parcellation.
+
+"""
+
+# Authors: Chetan Gohil <chetan.gohil@psych.ox.ac.uk>
+
+import os
+
+from osl import source_recon
+
+# Directories
+coreg_dir = "data/coreg"
+src_dir = "data/src"
+
+# First copy the coregistration directory
+if os.path.exists(src_dir):
+    print(f"Please first delete: {src_dir}")
+    exit()
+cmd = f"cp -r {coreg_dir} {src_dir}"
+print(cmd)
+os.system(cmd)
+
+# Settings
+config = """
+    source_recon:
+    - forward_model:
+        model: Single Layer
+    - beamform_and_parcellate:
+        freq_range: [1, 80]
+        chantypes: mag
+        rank: {mag: 120}
+        parcellation_file: Glasser52_binary_space-MNI152NLin6_res-8x8x8.nii.gz
+        method: spatial_basis
+        orthogonalisation: symmetric
+"""
+
+# List of subject IDs
+subjects = ["LN_VTA2"]
+
+# Fif files containing the sensor-level preprocessed data for each subject
+preproc_files = ["data/preproc/mg04938_BrainampDBS_20170504_01_preproc_raw.fif"]
+
+# Source reconstruction
+source_recon.run_src_batch(
+    config,
+    src_dir=src_dir,
+    subjects=subjects,
+    preproc_files=preproc_files,
+)

examples/ctf_with_smri_fid/4_sign_flip.py

@@ -0,0 +1,63 @@
+"""Align the sign of each parcel time course across subjects
+and save the data as a vanilla numpy file.
+
+"""
+
+# Authors: Chetan Gohil <chetan.gohil@psych.ox.ac.uk>
+
+import os
+import mne
+import numpy as np
+from glob import glob
+
+from osl.source_recon.sign_flipping import (
+    load_covariances,
+    find_template_subject,
+    find_flips,
+    apply_flips,
+)
+
+def load(filename):
+    """Load data without bad segments."""
+    raw = mne.io.read_raw_fif(filename, verbose=False)
+    raw = raw.pick("misc")
+    data = raw.get_data(reject_by_annotation="omit", verbose=False)
+    return data.T
+
+
+# Files to sign flip
+files = sorted(glob("data/src/*/parc/parc-raw.fif"))
+
+# Get covariance matrices
+covs = load_covariances(
+    files,
+    n_embeddings=15,
+    standardize=True,
+    loader=load,
+)
+
+# Load template covariance
+template_cov = np.load("../camcan_norm_model/template_cov.npy")
+
+# Output directory
+os.makedirs("data/npy", exist_ok=True)
+
+# Do sign flipping
+for i in range(len(files)):
+    print("Sign flipping", files[i])
+
+    # Find the channels to flip
+    flips, metrics = find_flips(
+        covs[i],
+        template_cov,
+        n_embeddings=15,
+        n_init=3,
+        n_iter=2500,
+        max_flips=20,
+    )
+
+    # Apply flips to the parcellated data and save
+    parc_data = load(files[i])
+    parc_data *= flips[np.newaxis, ...]
+    subject = files[i].split("/")[-3]
+    np.save(f"data/npy/{subject}.npy", parc_data)

examples/ctf_with_smri_fid/README.md

@@ -0,0 +1,16 @@
+# CTF Example
+
+In this example we coregister using sMRI fiducials.
+
+## Pipeline
+
+In this example we:
+
+- `1_preprocess.py`: Preprocess the sensor-level data.
+- `2_coregister.py`: Coregister the MEG and sMRI data.
+- `3_source_reconstruct.py`: Beamform the sensor-level data and parcellate to give us the source-level data.
+- `4_sign_flip.py`: Align the sign of each parcel time course to a template subject from the normative model.
+
+## Parallelisation
+
+See [here](https://github.com/OHBA-analysis/osl/tree/main/examples/parallelisation) for how to parallelise these scripts.