# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import calendar
import datetime
import os.path as op
import numpy as np
from scipy.spatial.distance import cdist
from ..._fiff._digitization import DigPoint, _make_dig_points
from ..._fiff.constants import FIFF
from ..._fiff.meas_info import _empty_info
from ..._fiff.utils import _read_segments_file
from ...transforms import Transform, apply_trans, get_ras_to_neuromag_trans
from ...utils import _check_fname, logger, verbose, warn
from ..base import BaseRaw
from .utils import _load_mne_locs, _read_pos
@verbose
def read_raw_artemis123(
input_fname, preload=False, verbose=None, pos_fname=None, add_head_trans=True
) -> "RawArtemis123":
"""Read Artemis123 data as raw object.
Parameters
----------
input_fname : path-like
Path to the data file (extension ``.bin``). The header file with the
same file name stem and an extension ``.txt`` is expected to be found
in the same directory.
%(preload)s
%(verbose)s
pos_fname : path-like | None
If not None, load digitized head points from this file.
add_head_trans : bool (default True)
If True attempt to perform initial head localization. Compute initial
device to head coordinate transform using HPI coils. If no
HPI coils are in info['dig'] hpi coils are assumed to be in canonical
order of fiducial points (nas, rpa, lpa).
Returns
-------
raw : instance of Raw
A Raw object containing the data.
See Also
--------
mne.io.Raw : Documentation of attributes and methods.
"""
return RawArtemis123(
input_fname,
preload=preload,
verbose=verbose,
pos_fname=pos_fname,
add_head_trans=add_head_trans,
)
def _get_artemis123_info(fname, pos_fname=None):
"""Generate info struct from artemis123 header file."""
fname = op.splitext(fname)[0]
header = fname + ".txt"
logger.info("Reading header...")
# key names for artemis channel info...
chan_keys = [
"name",
"scaling",
"FLL_Gain",
"FLL_Mode",
"FLL_HighPass",
"FLL_AutoReset",
"FLL_ResetLock",
]
header_info = dict()
header_info["filter_hist"] = []
header_info["comments"] = ""
header_info["channels"] = []
with open(header) as fid:
# section flag
# 0 - None
# 1 - main header
# 2 - channel header
# 3 - comments
# 4 - length
# 5 - filtering History
sectionFlag = 0
for line in fid:
# skip emptylines or header line for channel info
if (not line.strip()) or (sectionFlag == 2 and line.startswith("DAQ Map")):
continue
# set sectionFlag
if line.startswith("<end"):
sectionFlag = 0
elif line.startswith("<start main header>"):
sectionFlag = 1
elif line.startswith("<start per channel header>"):
sectionFlag = 2
elif line.startswith("<start comments>"):
sectionFlag = 3
elif line.startswith("<start length>"):
sectionFlag = 4
elif line.startswith("<start filtering history>"):
sectionFlag = 5
else:
# parse header info lines
# part of main header - lines are name value pairs
if sectionFlag == 1:
values = line.strip().split("\t")
if len(values) == 1:
values.append("")
header_info[values[0]] = values[1]
# part of channel header - lines are Channel Info
elif sectionFlag == 2:
values = line.strip().split("\t")
if len(values) != 7:
raise OSError(
f"Error parsing line \n\t:{line}\nfrom file {header}"
)
tmp = dict()
for k, v in zip(chan_keys, values):
tmp[k] = v
header_info["channels"].append(tmp)
elif sectionFlag == 3:
header_info["comments"] = f"{header_info['comments']}{line.strip()}"
elif sectionFlag == 4:
header_info["num_samples"] = int(line.strip())
elif sectionFlag == 5:
header_info["filter_hist"].append(line.strip())
for k in [
"Temporal Filter Active?",
"Decimation Active?",
"Spatial Filter Active?",
]:
if header_info[k] != "FALSE":
warn(f"{k} - set to but is not supported")
if header_info["filter_hist"]:
warn("Non-Empty Filter history found, BUT is not supported")
# build mne info struct
info = _empty_info(float(header_info["DAQ Sample Rate"]))
# Attempt to get time/date from fname
# Artemis123 files saved from the scanner observe the following
# naming convention 'Artemis_Data_YYYY-MM-DD-HHh-MMm_[chosen by user].bin'
try:
date = datetime.datetime.strptime(
op.basename(fname).split("_")[2], "%Y-%m-%d-%Hh-%Mm"
)
meas_date = (calendar.timegm(date.utctimetuple()), 0)
except Exception:
meas_date = None
# build subject info must be an integer (as per FIFF)
try:
subject_info = {"id": int(header_info["Subject ID"])}
except ValueError:
subject_info = {"id": 0}
# build description
desc = ""
for k in ["Purpose", "Notes"]:
desc += f"{k} : {header_info[k]}\n"
desc += f"Comments : {header_info['comments']}"
info.update(
{
"meas_date": meas_date,
"description": desc,
"subject_info": subject_info,
"proj_name": header_info["Project Name"],
}
)
# Channel Names by type
ref_mag_names = ["REF_001", "REF_002", "REF_003", "REF_004", "REF_005", "REF_006"]
ref_grad_names = ["REF_007", "REF_008", "REF_009", "REF_010", "REF_011", "REF_012"]
# load mne loc dictionary
loc_dict = _load_mne_locs()
info["chs"] = []
bads = []
for i, chan in enumerate(header_info["channels"]):
# build chs struct
t = {
"cal": float(chan["scaling"]),
"ch_name": chan["name"],
"logno": i + 1,
"scanno": i + 1,
"range": 1.0,
"unit_mul": FIFF.FIFF_UNITM_NONE,
"coord_frame": FIFF.FIFFV_COORD_DEVICE,
}
# REF_018 has a zero cal which can cause problems. Let's set it to
# a value of another ref channel to make writers/readers happy.
if t["cal"] == 0:
t["cal"] = 4.716e-10
bads.append(t["ch_name"])
t["loc"] = loc_dict.get(chan["name"], np.zeros(12))
if chan["name"].startswith("MEG"):
t["coil_type"] = FIFF.FIFFV_COIL_ARTEMIS123_GRAD
t["kind"] = FIFF.FIFFV_MEG_CH
# While gradiometer units are T/m, the meg sensors referred to as
# gradiometers report the field difference between 2 pick-up coils.
# Therefore the units of the measurements should be T
# *AND* the baseline (difference between pickup coils)
# should not be used in leadfield / forwardfield computations.
t["unit"] = FIFF.FIFF_UNIT_T
t["unit_mul"] = FIFF.FIFF_UNITM_F
# 3 axis reference magnetometers
elif chan["name"] in ref_mag_names:
t["coil_type"] = FIFF.FIFFV_COIL_ARTEMIS123_REF_MAG
t["kind"] = FIFF.FIFFV_REF_MEG_CH
t["unit"] = FIFF.FIFF_UNIT_T
t["unit_mul"] = FIFF.FIFF_UNITM_F
# reference gradiometers
elif chan["name"] in ref_grad_names:
t["coil_type"] = FIFF.FIFFV_COIL_ARTEMIS123_REF_GRAD
t["kind"] = FIFF.FIFFV_REF_MEG_CH
# While gradiometer units are T/m, the meg sensors referred to as
# gradiometers report the field difference between 2 pick-up coils.
# Therefore the units of the measurements should be T
# *AND* the baseline (difference between pickup coils)
# should not be used in leadfield / forwardfield computations.
t["unit"] = FIFF.FIFF_UNIT_T
t["unit_mul"] = FIFF.FIFF_UNITM_F
# other reference channels are unplugged and should be ignored.
elif chan["name"].startswith("REF"):
t["coil_type"] = FIFF.FIFFV_COIL_NONE
t["kind"] = FIFF.FIFFV_MISC_CH
t["unit"] = FIFF.FIFF_UNIT_V
bads.append(t["ch_name"])
elif chan["name"].startswith(("AUX", "TRG", "MIO")):
t["coil_type"] = FIFF.FIFFV_COIL_NONE
t["unit"] = FIFF.FIFF_UNIT_V
if chan["name"].startswith("TRG"):
t["kind"] = FIFF.FIFFV_STIM_CH
else:
t["kind"] = FIFF.FIFFV_MISC_CH
else:
raise ValueError(
f'Channel does not match expected channel Types:"{chan["name"]}"'
)
# incorporate multiplier (unit_mul) into calibration
t["cal"] *= 10 ** t["unit_mul"]
t["unit_mul"] = FIFF.FIFF_UNITM_NONE
# append this channel to the info
info["chs"].append(t)
if chan["FLL_ResetLock"] == "TRUE":
bads.append(t["ch_name"])
# HPI information
# print header_info.keys()
hpi_sub = dict()
# Don't know what event_channel is don't think we have it HPIs are either
# always on or always off.
# hpi_sub['event_channel'] = ???
hpi_sub["hpi_coils"] = [dict(), dict(), dict(), dict()]
hpi_coils = [dict(), dict(), dict(), dict()]
drive_channels = ["MIO_001", "MIO_003", "MIO_009", "MIO_011"]
key_base = "Head Tracking %s %d"
# set default HPI frequencies
if info["sfreq"] == 1000:
default_freqs = [140, 150, 160, 40]
else:
default_freqs = [700, 750, 800, 40]
for i in range(4):
# build coil structure
hpi_coils[i]["number"] = i + 1
hpi_coils[i]["drive_chan"] = drive_channels[i]
this_freq = header_info.pop(key_base % ("Frequency", i + 1), default_freqs[i])
hpi_coils[i]["coil_freq"] = this_freq
# check if coil is on
if header_info[key_base % ("Channel", i + 1)] == "OFF":
hpi_sub["hpi_coils"][i]["event_bits"] = [0]
else:
hpi_sub["hpi_coils"][i]["event_bits"] = [256]
info["hpi_subsystem"] = hpi_sub
info["hpi_meas"] = [{"hpi_coils": hpi_coils}]
# read in digitized points if supplied
if pos_fname is not None:
info["dig"] = _read_pos(pos_fname)
else:
info["dig"] = []
info._unlocked = False
info._update_redundant()
# reduce info['bads'] to unique set
info["bads"] = list(set(bads))
del bads
return info, header_info
class RawArtemis123(BaseRaw):
"""Raw object from Artemis123 file.
Parameters
----------
input_fname : path-like
Path to the Artemis123 data file (ending in ``'.bin'``).
%(preload)s
%(verbose)s
See Also
--------
mne.io.Raw : Documentation of attributes and methods.
"""
@verbose
def __init__(
self,
input_fname,
preload=False,
verbose=None,
pos_fname=None,
add_head_trans=True,
):
from ...chpi import (
_fit_coil_order_dev_head_trans,
compute_chpi_amplitudes,
compute_chpi_locs,
)
input_fname = str(_check_fname(input_fname, "read", True, "input_fname"))
fname, ext = op.splitext(input_fname)
if ext == ".txt":
input_fname = fname + ".bin"
elif ext != ".bin":
raise RuntimeError(
'Valid artemis123 files must end in "txt"' + ' or ".bin".'
)
if not op.exists(input_fname):
raise RuntimeError(f"{input_fname} - Not Found")
info, header_info = _get_artemis123_info(input_fname, pos_fname=pos_fname)
last_samps = [header_info.get("num_samples", 1) - 1]
super().__init__(
info,
preload,
filenames=[input_fname],
raw_extras=[header_info],
last_samps=last_samps,
orig_format="single",
verbose=verbose,
)
if add_head_trans:
n_hpis = 0
for d in info["hpi_subsystem"]["hpi_coils"]:
if d["event_bits"] == [256]:
n_hpis += 1
if n_hpis < 3:
warn(
f"{n_hpis:d} HPIs active. At least 3 needed to perform"
"head localization\n *NO* head localization performed"
)
else:
# Localized HPIs using the 1st 250 milliseconds of data.
with info._unlock():
info["hpi_results"] = [
dict(
dig_points=[
dict(
r=np.zeros(3),
coord_frame=FIFF.FIFFV_COORD_DEVICE,
ident=ii + 1,
)
for ii in range(n_hpis)
],
coord_trans=Transform("meg", "head"),
)
]
coil_amplitudes = compute_chpi_amplitudes(
self, tmin=0, tmax=0.25, t_window=0.25, t_step_min=0.25
)
assert len(coil_amplitudes["times"]) == 1
coil_locs = compute_chpi_locs(self.info, coil_amplitudes)
with info._unlock():
info["hpi_results"] = None
hpi_g = coil_locs["gofs"][0]
hpi_dev = coil_locs["rrs"][0]
# only use HPI coils with localizaton goodness_of_fit > 0.98
bad_idx = []
for i, g in enumerate(hpi_g):
msg = f"HPI coil {i + 1} - location goodness of fit ({g:0.3f})"
if g < 0.98:
bad_idx.append(i)
msg += " *Removed from coregistration*"
logger.info(msg)
hpi_dev = np.delete(hpi_dev, bad_idx, axis=0)
hpi_g = np.delete(hpi_g, bad_idx, axis=0)
if pos_fname is not None:
# Digitized HPI points are needed.
hpi_head = np.array(
[
d["r"]
for d in self.info.get("dig", [])
if d["kind"] == FIFF.FIFFV_POINT_HPI
]
)
if len(hpi_head) != len(hpi_dev):
raise RuntimeError(
f"number of digitized ({len(hpi_head)}) and active "
f"({len(hpi_dev)}) HPI coils are not the same."
)
# compute initial head to dev transform and hpi ordering
head_to_dev_t, order, trans_g = _fit_coil_order_dev_head_trans(
hpi_dev, hpi_head
)
# set the device to head transform
self.info["dev_head_t"] = Transform(
FIFF.FIFFV_COORD_DEVICE, FIFF.FIFFV_COORD_HEAD, head_to_dev_t
)
# add hpi_meg_dev to dig...
for idx, point in enumerate(hpi_dev):
d = {
"r": point,
"ident": idx + 1,
"kind": FIFF.FIFFV_POINT_HPI,
"coord_frame": FIFF.FIFFV_COORD_DEVICE,
}
self.info["dig"].append(DigPoint(d))
dig_dists = cdist(hpi_head[order], hpi_head[order])
dev_dists = cdist(hpi_dev, hpi_dev)
tmp_dists = np.abs(dig_dists - dev_dists)
dist_limit = tmp_dists.max() * 1.1
logger.info(
"HPI-Dig corrregsitration\n"
f"\tGOF : {trans_g:0.3f}\n"
f"\tMax Coil Error : {100 * tmp_dists.max():0.3f} cm\n"
)
else:
logger.info("Assuming Cardinal HPIs")
nas = hpi_dev[0]
lpa = hpi_dev[2]
rpa = hpi_dev[1]
t = get_ras_to_neuromag_trans(nas, lpa, rpa)
with self.info._unlock():
self.info["dev_head_t"] = Transform(
FIFF.FIFFV_COORD_DEVICE, FIFF.FIFFV_COORD_HEAD, t
)
# transform fiducial points
nas = apply_trans(t, nas)
lpa = apply_trans(t, lpa)
rpa = apply_trans(t, rpa)
hpi = apply_trans(self.info["dev_head_t"], hpi_dev)
with self.info._unlock():
self.info["dig"] = _make_dig_points(
nasion=nas, lpa=lpa, rpa=rpa, hpi=hpi
)
order = np.array([0, 1, 2])
dist_limit = 0.005
# fill in hpi_results
hpi_result = dict()
# add HPI points in device coords...
dig = []
for idx, point in enumerate(hpi_dev):
dig.append(
{
"r": point,
"ident": idx + 1,
"kind": FIFF.FIFFV_POINT_HPI,
"coord_frame": FIFF.FIFFV_COORD_DEVICE,
}
)
hpi_result["dig_points"] = dig
# attach Transform
hpi_result["coord_trans"] = self.info["dev_head_t"]
# 1 based indexing
hpi_result["order"] = order + 1
hpi_result["used"] = np.arange(3) + 1
hpi_result["dist_limit"] = dist_limit
hpi_result["good_limit"] = 0.98
# Warn for large discrepancies between digitized and fit
# cHPI locations
if hpi_result["dist_limit"] > 0.005:
warn(
"Large difference between digitized geometry"
" and HPI geometry. Max coil to coil difference"
f" is {100.0 * tmp_dists.max():0.2f} cm\n"
"beware of *POOR* head localization"
)
# store it
with self.info._unlock():
self.info["hpi_results"] = [hpi_result]
def _read_segment_file(self, data, idx, fi, start, stop, cals, mult):
"""Read a chunk of raw data."""
_read_segments_file(self, data, idx, fi, start, stop, cals, mult, dtype=">f4")
Datasets
Models
