source

tit.source

EEG source-forward preparation and fsaverage field mapping.

Two pipelines that put EEG source reconstruction and TI stimulation fields on a common cortical grid:

• :func:prepare_forward -- build an MNE-compatible EEG forward solution (leadfield + source space + fsaverage morph) from a SimNIBS head model.
• :func:project_fields_to_fsaverage -- project existing simulation field outputs (TI_max, TI_normal, hf_peak, hf_rms) onto an fsaverage template.

Both run under the SimNIBS interpreter::

simnibs_python -m tit.source config.json

See Also

tit.source.config : ForwardConfig and FsavgMapConfig dataclasses.

ForwardConfig dataclass

ForwardConfig(eeg_net: str = 'GSN-HydroCel-185', fsaverage_spacing: int = 5, cpus: int = 1, overwrite: bool = False)

Parameters for rebuilding a SimNIBS/MNE EEG forward solution.

Attributes

eeg_net : str EEG cap name without the .csv suffix, as found in the subject's m2m_<id>/eeg_positions/ directory (e.g. "GSN-HydroCel-185"). fsaverage_spacing : int fsaverage subdivision factor (5, 6, or 7) for the morph target. cpus : int SimNIBS FEM workers used while computing the leadfield. overwrite : bool Recompute even when valid outputs already exist. The expensive FEM leadfield is still reused when present unless this is True.

FsavgMapConfig dataclass

FsavgMapConfig(fields: tuple[str, ...] = (lambda: VALID_FSAVG_FIELDS)(), fsaverage_spacing: int = 5, workers: int = 1, overwrite: bool = False)

Parameters for projecting simulation field outputs onto fsaverage.

Attributes

fields : tuple of str Which field quantities to project. Any of :data:VALID_FSAVG_FIELDS ("TI_max", "TI_normal", "hf_peak", "hf_rms"). fsaverage_spacing : int fsaverage subdivision factor (5, 6, or 7) to morph onto. workers : int Number of subjects projected in parallel (1 = serial). overwrite : bool Re-project even when the fsaverage .msh already exists.

project_fields_to_fsaverage

project_fields_to_fsaverage(subjects: list[tuple[str, str]], cfg: FsavgMapConfig) -> list[tuple[str, str, str]]

Project field outputs to fsaverage for many (subject, simulation) pairs.

Parameters

subjects : list of (str, str) (subject_id, simulation_name) pairs to project. cfg : FsavgMapConfig Field selection, fsaverage spacing, worker count, overwrite flag.

Returns

list of (str, str, str) Per-pair (subject_id, status, message) results.

Source code in tit/source/fsaverage.py

def project_fields_to_fsaverage(
    subjects: list[tuple[str, str]],
    cfg: FsavgMapConfig,
) -> list[tuple[str, str, str]]:
    """Project field outputs to fsaverage for many (subject, simulation) pairs.

    Parameters
    ----------
    subjects : list of (str, str)
        ``(subject_id, simulation_name)`` pairs to project.
    cfg : FsavgMapConfig
        Field selection, fsaverage spacing, worker count, overwrite flag.

    Returns
    -------
    list of (str, str, str)
        Per-pair ``(subject_id, status, message)`` results.
    """
    workers = max(1, min(cfg.workers, len(subjects))) if subjects else 1
    logger.info("Projecting %d simulation(s) with %d worker(s)", len(subjects), workers)

    results: list[tuple[str, str, str]] = []
    if workers == 1:
        for subject_id, sim in subjects:
            res = project_subject(subject_id, sim, cfg)
            _log_result(res)
            results.append(res)
    else:
        with ProcessPoolExecutor(max_workers=workers) as pool:
            futures = {
                pool.submit(project_subject, sid, sim, cfg): (sid, sim)
                for sid, sim in subjects
            }
            for future in as_completed(futures):
                res = future.result()
                _log_result(res)
                results.append(res)
    return results

project_subject

project_subject(subject_id: str, sim: str, cfg: FsavgMapConfig) -> tuple[str, str, str]

Project one (subject, simulation) and write it as an fsaverage .msh.

Returns (subject_id, status, message) where status is one of {"ok", "partial", "cached", "failed"} so batch runs can record and continue; "partial" means some requested fields could not be projected.

Source code in tit/source/fsaverage.py

def project_subject(
    subject_id: str,
    sim: str,
    cfg: FsavgMapConfig,
) -> tuple[str, str, str]:
    """Project one (subject, simulation) and write it as an fsaverage ``.msh``.

    Returns ``(subject_id, status, message)`` where status is one of
    ``{"ok", "partial", "cached", "failed"}`` so batch runs can record and
    continue; "partial" means some requested fields could not be projected.
    """
    from dataclasses import replace

    from tit.source.config import VALID_FSAVG_FIELDS
    from tit.tools.msh_fields import msh_field_names

    pm = get_path_manager()
    out_path = Path(pm.sim_fsaverage_fields(subject_id, sim, cfg.fsaverage_spacing))
    try:
        carrier_only = _recorded_carrier_only(pm, subject_id, sim)
        if carrier_only:
            # A carrier-only control has no TI envelope: TI_max / TI_normal cannot
            # exist. The default (all-fields) request means "what this run has";
            # an explicit selection naming a TI field is refused.
            ti = {"TI_max", "TI_normal"} & set(cfg.fields)
            if ti and tuple(cfg.fields) != VALID_FSAVG_FIELDS:
                raise ValueError(
                    f"{sorted(ti)} cannot be projected: {sim} is a carrier-only run "
                    "with no TI envelope; request hf_peak / hf_rms."
                )
            cfg = replace(cfg, fields=tuple(f for f in cfg.fields if f not in ti))
        if out_path.exists() and not cfg.overwrite:
            names = msh_field_names(out_path)  # None: could not tell, keep the file
            if names is None or set(cfg.fields) <= names:
                return subject_id, "cached", out_path.name
        maps = _compute_fields(pm, subject_id, sim, cfg, carrier_only)
    except Exception as exc:  # noqa: BLE001 - record per-subject and continue
        return subject_id, "failed", repr(exc)
    write_fsaverage_msh(
        out_path,
        maps,
        cfg.fsaverage_spacing,
        {"subject_id": subject_id, "simulation": sim, "carrier_only": carrier_only},
    )
    medians = ", ".join(f"{k} med={np.median(v):.3f}" for k, v in maps.items())
    missing = [f for f in cfg.fields if f not in maps]
    if missing:
        return subject_id, "partial", f"{medians} -> {out_path.name}; missing {missing}"
    return subject_id, "ok", f"{medians} -> {out_path.name}"

prepare_forward

prepare_forward(subject_id: str, cfg: ForwardConfig, *, output_dir: str | Path | None = None, info: Info | None = None, trans: Transform | None = None) -> tuple[Path, Path, Path]

Rebuild one subject's forward, source-space, and fsaverage morph files.

Parameters

subject_id : str Subject identifier without the sub- prefix (e.g. "101"). cfg : ForwardConfig EEG net, fsaverage spacing, FEM worker count, and overwrite flag. output_dir : str or pathlib.Path or None Destination directory. Defaults to pm.forward(subject_id) (derivatives/SimNIBS/sub-<id>/forward/). info : mne.Info or None A recording's Info whose dig carries the LPA, Nasion and RPA cardinal points. The head<->MRI transform is coregistered from them, so the forward is in the recording's head frame. Default: the net's own template fiducials. trans : mne.transforms.Transform or None An explicit head->MRI transform (MNE head -> MRI frames); used as is, never re-estimated. Pass info or trans, not both.

Returns

tuple of pathlib.Path (fwd_path, src_path, morph_path).

Source code in tit/source/forward.py

def prepare_forward(
    subject_id: str,
    cfg: ForwardConfig,
    *,
    output_dir: str | Path | None = None,
    info: mne.Info | None = None,
    trans: mne.transforms.Transform | None = None,
) -> tuple[Path, Path, Path]:
    """Rebuild one subject's forward, source-space, and fsaverage morph files.

    Parameters
    ----------
    subject_id : str
        Subject identifier without the ``sub-`` prefix (e.g. ``"101"``).
    cfg : ForwardConfig
        EEG net, fsaverage spacing, FEM worker count, and overwrite flag.
    output_dir : str or pathlib.Path or None
        Destination directory.  Defaults to ``pm.forward(subject_id)``
        (``derivatives/SimNIBS/sub-<id>/forward/``).
    info : mne.Info or None
        A recording's Info whose ``dig`` carries the LPA, Nasion and RPA
        cardinal points.  The head<->MRI transform is coregistered from them,
        so the forward is in the recording's head frame.  Default: the net's
        own template fiducials.
    trans : mne.transforms.Transform or None
        An explicit head->MRI transform (MNE head -> MRI frames); used as is,
        never re-estimated.  Pass ``info`` or ``trans``, not both.

    Returns
    -------
    tuple of pathlib.Path
        ``(fwd_path, src_path, morph_path)``.
    """
    _check_recording_geometry(info, trans)
    pm = get_path_manager()
    m2m_dir = Path(pm.m2m(subject_id))
    if not m2m_dir.is_dir():
        raise FileNotFoundError(f"m2m directory not found: {m2m_dir}")

    _check_forward_dependencies()

    forward_dir = (
        Path(output_dir) if output_dir is not None else Path(pm.forward(subject_id))
    )
    forward_dir.mkdir(parents=True, exist_ok=True)

    stem = f"sub-{subject_id}_net-{cfg.eeg_net}"
    expected = (
        forward_dir / f"{stem}-fwd.fif",
        forward_dir / f"{stem}-src.fif",
        forward_dir / f"{stem}-morph.h5",
    )
    electrodes, fiducials = _read_simnibs_montage(m2m_dir, cfg.eeg_net)

    # Cached files are reused only if built from these exact inputs: a net CSV
    # rewritten in place (the v3.1.1 cap re-apply) or a new head mesh rebuilds.
    eeg_csv = m2m_dir / "eeg_positions" / f"{cfg.eeg_net}.csv"
    provenance_path = forward_dir / f"{stem}-provenance.json"
    inputs = {
        "net_csv_sha256": sha256_file(eeg_csv),
        "head_mesh_sha256": sha256_file(m2m_dir / f"{subject_id}.msh"),
    }
    recorded = {} if cfg.overwrite else _read_provenance(provenance_path)
    same_inputs = all(recorded.get(k) == v for k, v in inputs.items())
    # Cached outputs may hold another recording's geometry, so a supplied
    # ``info``/``trans`` always rebuilds (the FEM leadfield is still reused).
    default_frame = info is None and trans is None
    if (
        same_inputs
        and default_frame
        and recorded.get("fsaverage_spacing") == cfg.fsaverage_spacing
        and _forward_outputs_valid(*expected)
    ):
        logger.info("Forward outputs already present for %s; skipping.", subject_id)
        return expected

    # SimNIBS names the leadfield <head mesh stem>_leadfield_<net>.hdf5.
    leadfield_hdf5 = forward_dir / f"{subject_id}_leadfield_{cfg.eeg_net}.hdf5"
    if same_inputs and _exists_real(leadfield_hdf5):
        logger.info("Reusing cached leadfield: %s", leadfield_hdf5.name)
    else:
        provenance_path.unlink(missing_ok=True)
        leadfield_hdf5.unlink(missing_ok=True)  # never append to a stale HDF5
        leadfield_hdf5 = _compute_leadfield(m2m_dir, forward_dir, eeg_csv, cfg.cpus)
        _write_provenance(provenance_path, inputs)

    # The Info is built *after* the leadfield, in the leadfield's channel order:
    # SimNIBS's assembly asserts `info["ch_names"] == forward["ch_names"]` element by
    # element, and the leadfield does not use the net CSV's order (below).
    electrodes = _in_leadfield_order(electrodes, leadfield_hdf5)
    montage_info, trans = _build_montage_info(
        electrodes, fiducials, info=info, trans=trans
    )

    info_path = forward_dir / f"{stem}-info.fif"
    trans_path = forward_dir / f"{stem}-trans.fif"
    mne.io.RawArray(
        np.zeros((len(montage_info.ch_names), 1)), montage_info, verbose=False
    ).save(str(info_path), overwrite=True, verbose=False)
    mne.write_trans(str(trans_path), trans, overwrite=True)

    # Not SimNIBS's own ``prepare_eeg_forward`` console script: that wrapper runs
    # ``python -E``, so no compatibility shim can reach it, and its MNE/cortech
    # bridge dies *after* the FEM leadfield above has been paid for. The worker
    # calls the same ``simnibs.eeg.forward.make_forward`` with the shims applied --
    # still a child process, so the ~0.9 GB gain matrix is reclaimed on exit.
    # See tit/source/_simnibs_compat.py.
    _run(
        [
            sys.executable,
            "-m",
            "tit.source._prepare_forward",
            str(m2m_dir),
            str(leadfield_hdf5),
            str(info_path),
            str(trans_path),
            "--fsaverage",
            str(cfg.fsaverage_spacing),
        ],
        f"MNE forward assembly for sub-{subject_id}",
        cwd=forward_dir,
    )
    outputs = _rename_generated_outputs(forward_dir, stem, leadfield_hdf5.stem)
    provenance = dict(inputs)
    if default_frame:  # a recording-frame forward is never reused as the default
        provenance["fsaverage_spacing"] = cfg.fsaverage_spacing
    _write_provenance(provenance_path, provenance)
    return outputs