diff --git a/mne/io/__init__.pyi b/mne/io/__init__.pyi index 11b86623b8e..11528cce278 100644 --- a/mne/io/__init__.pyi +++ b/mne/io/__init__.pyi @@ -37,6 +37,7 @@ __all__ = [ "read_raw_fil", "read_raw_gdf", "read_raw_hitachi", + "read_raw_itab", "read_raw_kit", "read_raw_mef", "read_raw_nedf", @@ -83,6 +84,7 @@ from .fieldtrip import read_epochs_fieldtrip, read_evoked_fieldtrip, read_raw_fi from .fiff import Raw, read_raw_fif from .fil import read_raw_fil from .hitachi import read_raw_hitachi +from .itab import read_raw_itab from .kit import read_epochs_kit, read_raw_kit from .mef import read_raw_mef from .nedf import read_raw_nedf diff --git a/mne/io/itab/__init__.py b/mne/io/itab/__init__.py new file mode 100755 index 00000000000..9204f6a219c --- /dev/null +++ b/mne/io/itab/__init__.py @@ -0,0 +1,7 @@ +"""ITAB module for conversion to FIF""" + +# Author: Vittorio Pizzella +# +# License: BSD (3-clause) + +from .itab import read_raw_itab, RawITAB \ No newline at end of file diff --git a/mne/io/itab/constants.py b/mne/io/itab/constants.py new file mode 100755 index 00000000000..eef36079785 --- /dev/null +++ b/mne/io/itab/constants.py @@ -0,0 +1,34 @@ +"""ITAB constants""" + +# Author: Vittorio Pizzella +# +# License: BSD (3-clause) + +from ...utils import BunchConst + + +ITAB = BunchConst() + + +# Channel types +ITAB.ITABV_EEG_CH = 1 +ITAB.ITABV_MAG_CH = 2 +ITAB.ITABV_REF_EEG_CH = 4 +ITAB.ITABV_REF_MAG_CH = 8 +ITAB.ITABV_AUX_CH = 16 +ITAB.ITABV_PARAM_CH =32 +ITAB.ITABV_DIGIT_CH = 64 +ITAB.ITABV_FLAG_CH = 128 + + +# mhd file pointers +ITAB.MHD_TIME = 656 +ITAB.MHD_NREFCH = 82780 +ITAB.MHD_RAWHDTYPE = 85348 +ITAB.MHD_CHINFO = 85444 +ITAB.MHD_SENSPOS = 295364 +ITAB.MHD_NMARKER = 295440 +ITAB.MHD_MARKERS = 295700 +ITAB.MHD_BESTCHI = 297232 + +ITAB.MHD_CH_SIZE = 328 \ No newline at end of file diff --git a/mne/io/itab/info.py b/mne/io/itab/info.py new file mode 100755 index 00000000000..3083f582f81 --- /dev/null +++ b/mne/io/itab/info.py @@ -0,0 +1,281 @@ +"""Build measurement info.""" + +# Author: Vittorio Pizzella +# +# License: BSD (3-clause) + +from datetime import datetime, timezone + +import numpy as np + +from ..._fiff._digitization import DigPoint +from ..._fiff.meas_info import _empty_info +from ...annotations import Annotations +from ...utils import logger +from ..constants import FIFF +from .constants import ITAB + + +def _convert_time(date_str, time_str): + """Convert date and time strings to float time.""" + for fmt in ("%d/%m/%Y", "%d-%b-%Y", "%a, %b %d, %Y"): + try: + date = datetime.strptime(date_str, fmt) + except ValueError: + pass + else: + break + else: + raise RuntimeError( + f"Illegal date: {date_str}.\nIf the language of the date does not " + "correspond to your local machine's language try to set the " + "locale to the language of the date string:\n" + 'locale.setlocale(locale.LC_ALL, "en_US")' + ) + + for fmt in ("%H:%M:%S", "%H:%M"): + try: + time = datetime.strptime(time_str, fmt) + except ValueError: + pass + else: + break + else: + raise RuntimeError(f"Illegal time: {time_str}.") + # MNE-C uses mktime which uses local time, but here we instead decouple + # conversion location from the process, and instead assume that the + # acquisiton was in GMT. This will be wrong for most sites, but at least + # the value we obtain here won't depend on the geographical location + # that the file was converted. + res = datetime( + date.year, + date.month, + date.day, + time.hour, + time.minute, + time.second, + tzinfo=timezone.utc, + ) + + return res + + +def _mhdch_2_chs(mhd_ch): + """Build chs list item from mhd ch list item.""" + ch = dict() + + unit = 1.0 # Default unit multiplier + # Generic channel + loc = np.zeros(12) + ch["ch_name"] = mhd_ch["label"] + ch["coord_frame"] = FIFF.FIFFV_COORD_UNKNOWN + ch["coil_type"] = 0 + ch["range"] = 1.0 + ch["unit"] = FIFF.FIFF_UNIT_NONE + ch["unit_mul"] = FIFF.FIFF_UNITM_NONE + if mhd_ch["calib"] == 0: + ch["cal"] = 1.0 + else: + ch["cal"] = float(mhd_ch["amvbit"] / mhd_ch["calib"]) + ch["loc"] = loc + ch["scanno"] = 0 + ch["kind"] = FIFF.FIFFV_MISC_CH + ch["logno"] = 1 + + # Magnetic channel + if mhd_ch["type"] == ITAB.ITABV_MAG_CH: + loc[0:12] = ( + mhd_ch["pos"][0]["posx"] / 1000, # r0 + mhd_ch["pos"][0]["posy"] / 1000, + mhd_ch["pos"][0]["posz"] / 1000, + mhd_ch["pos"][0]["orix"], # ex + 0, + 0, + 0, # ey + mhd_ch["pos"][0]["oriy"], + 0, + 0, # ez + 0, + mhd_ch["pos"][0]["oriz"], + ) + ch["loc"] = loc + ch["kind"] = FIFF.FIFFV_MEG_CH + ch["coil_type"] = FIFF.FIFFV_COIL_POINT_MAGNETOMETER + ch["logno"] = int(mhd_ch["number"]) + if mhd_ch["calib"] == 0: + ch["cal"] = 1.0 + else: + ch["cal"] = float(mhd_ch["amvbit"] / mhd_ch["calib"]) + ch["unit"] = FIFF.FIFF_UNIT_T + if mhd_ch["unit"] == "fT": + ch["unit_mul"] = FIFF.FIFF_UNITM_F + unit = 10**-15 + elif mhd_ch["unit"] == "pT": + ch["unit_mul"] = FIFF.FIFF_UNITM_P + unit = 10**-12 + + # Electric channel + if mhd_ch["type"] == ITAB.ITABV_EEG_CH: + ch["kind"] = FIFF.FIFFV_BIO_CH + ch["cal"] = float(mhd_ch["amvbit"] / mhd_ch["calib"]) + ch["unit"] = FIFF.FIFF_UNIT_V + ch["logno"] = int(mhd_ch["number"]) + if mhd_ch["unit"] == "mV": + ch["unit_mul"] = FIFF.FIFF_UNITM_M + unit = 10**-3 + elif mhd_ch["unit"] == "uT": + ch["unit_mul"] = FIFF.FIFF_UNITM_MU + unit = 10**-6 + + # Other channel type + if ( + mhd_ch["type"] == ITAB.ITABV_REF_EEG_CH + and mhd_ch["type"] == ITAB.ITABV_REF_MAG_CH + and mhd_ch["type"] == ITAB.ITABV_REF_AUX_CH + and mhd_ch["type"] == ITAB.ITABV_REF_PARAM_CH + and mhd_ch["type"] == ITAB.ITABV_REF_DIGIT_CH + and mhd_ch["type"] == ITAB.ITABV_REF_FLAG_CH + ): + ch["kind"] = FIFF.FIFFV_BIO_CH + ch["cal"] = float(mhd_ch["amvbit"] / mhd_ch["calib"]) + ch["unit"] = FIFF.FIFF_UNIT_V + if mhd_ch["unit"] == "mV": + ch["unit_mul"] = FIFF.FIFF_UNITM_M + unit = 10**-3 + elif mhd_ch["unit"] == "uT": + ch["unit_mul"] = FIFF.FIFF_UNITM_MU + unit = 10**-6 + + return ch, unit + + +def _mhd2info(mhd): + """Create meas info from ITAB mhd data.""" + info = _empty_info(mhd["smpfq"]) + + info["meas_date"] = _convert_time(mhd["date"], mhd["time"]) + + info["description"] = mhd["notes"] + + si = dict() + # si['id'] = mhd['id'] + si["last_name"] = mhd["last_name"] + si["first_name"] = mhd["first_name"] + si["sex"] = int(mhd["subinfo"]["sex"][0]) + + info["subject_info"] = si + + channel_names = list() + chs = list() + bads = list() + units = list() + + for channel in mhd["ch"]: + channel_names.append(channel["label"]) + + if channel["flag"] == 1: + bads.append(channel["label"]) + + ch, unit = _mhdch_2_chs(channel) + chs.append(ch) + + units.append(unit) + + info["chs"] = chs + info["ch_names"] = channel_names + info["bads"] = bads + + if mhd["hw_low_fr"] != 0: + info["lowpass"] = mhd["hw_low_fr"] + + if mhd["hw_hig_fr"] != 0: + info["highpass"] = mhd["hw_hig_fr"] + + # Total number of channels in .raw file + info["nchan"] = mhd["nchan"] + + # Total number of timepoints in .raw file + info["temp"] = dict() + + info["temp"]["units"] = units + info["temp"]["n_samp"] = mhd["ntpdata"] + + # Only one trial (continous acquisition) + info["temp"]["ntrials"] = 1 + + # Data start in .raw file + info["temp"]["start_data"] = mhd["start_data"] + + # Get Polhemus digitization data (in head coordinates) + dig = [] + + point_sequence = [ + FIFF.FIFFV_POINT_NASION, # Nasion + FIFF.FIFFV_POINT_RPA, # Right pre-auricolar + FIFF.FIFFV_POINT_LPA, # Left pre-auricolar + ] + + for i, point in enumerate(mhd["marker"]): + point_info = DigPoint() + + point_info["coord_frame"] = FIFF.FIFFV_COORD_HEAD + + # Point kind + if i < 3: # Nasion, Right pre-auricolar, Left pre-auricolar + point_kind = FIFF.FIFFV_POINT_CARDINAL + elif i == 3: # Vertex + point_kind = FIFF.FIFFV_POINT_EXTRA + else: # HPI Coils + point_kind = FIFF.FIFFV_POINT_HPI + + point_info["kind"] = point_kind + + # Point identity + if i < 3: + point_ident = point_sequence[i] + else: + point_ident = i + 1 + point_info["ident"] = point_ident + + point_info["r"] = np.array([point["posx"], point["posy"], point["posz"]]) + + dig += [point_info] + + # TDB other poosible head points, check on dig points. + + info["dig"] = dig + + # TBD trigger handling + + event = list() + for sample in mhd["sample"]: + event.append([sample["start"], sample["type"], sample["quality"]]) + + info["events"] = event + + annotations = { + "description": [], + "onset": [], + "duration": [], + } + for sample in mhd["segment"]: + if sample["start"] != 0: + start = sample["start"] / info["sfreq"] + else: + start = 0.0 + duration = sample["ntptot"] / info["sfreq"] + description = sample["type"] + + annotations["description"].append(description) + annotations["onset"].append(start) + annotations["duration"].append(duration) + + info["temp"]["annotations"] = Annotations( + annotations["onset"], annotations["duration"], annotations["description"] + ) + + info._check_consistency() + + logger.info("Measurement info composed.") + + return info diff --git a/mne/io/itab/itab.py b/mne/io/itab/itab.py new file mode 100755 index 00000000000..de409142c0c --- /dev/null +++ b/mne/io/itab/itab.py @@ -0,0 +1,182 @@ +# Author: Vittorio Pizzella +# +# License: BSD (3-clause) + +import numpy as np + +from ..._fiff.utils import _mult_cal_one +from ...utils import ( + _check_fname, + fill_doc, + verbose, +) +from ..base import BaseRaw +from .info import _mhd2info +from .mhd import _read_mhd + + +class RawITAB(BaseRaw): + """Raw object from ITAB directory. + + Parameters + ---------- + fname : str + The raw file to load. Filename should end with *.raw + preload : bool or str (default False) + Preload data into memory for data manipulation and faster indexing. + If True, the data will be preloaded into memory (fast, requires + large amount of memory). If preload is a string, preload is the + file name of a memory-mapped file which is used to store the data + on the hard drive (slower, requires less memory). + verbose : bool, str, int, or None + If not None, override default verbose level (see mne.verbose). + + See Also + -------- + mne.io.Raw : Documentation of attribute and methods. + """ + + @verbose + def __init__(self, fname, preload=False, verbose=None): + + filenames = list() + filenames.append(fname) + + fname = _check_fname(fname, overwrite="read", must_exist=True) + + fname_mhd = fname.with_name(fname.name + ".mhd") + try: + mhd = _read_mhd(fname_mhd) # Read the mhd file + except FileNotFoundError: + raise ValueError(".mhd file not found") + + info = _mhd2info(mhd) + + orig_units = {ch["label"]: ch["unit"] for ch in mhd["ch"]} + + raw_extras = list() + for fi, _ in enumerate(filenames): + raw_extras.append(dict()) + for k in ["n_samp", "start_data", "units"]: + raw_extras[fi][k] = info["temp"][k] + + raw_extras[fi]["nchan"] = info["nchan"] + raw_extras[fi]["buffer_size_sec"] = info["temp"]["n_samp"] / info["sfreq"] + raw_extras[fi]["data_type"] = mhd["data_type"] + + self.info = info + info._check_consistency() + + first_samps = [0] + last_samps = [info["temp"]["n_samp"] - 1] + + annotations = info["temp"]["annotations"] + + # Remove extras from info + del info["temp"] + + super().__init__( + info, + preload, + first_samps=first_samps, + last_samps=last_samps, + raw_extras=raw_extras, + filenames=filenames, + orig_units=orig_units, + verbose=verbose, + ) + + self.set_annotations(annotations) + + def _read_segment_file(self, data, idx, fi, start, stop, cals, mult): + """Read a segment of data from a file. + + Only needs to be implemented for readers that support + ``preload=False``. + + Parameters + ---------- + data : ndarray, shape (len(idx), stop - start + 1) + The data array. Should be modified inplace. + idx : ndarray | slice + The requested channel indices. + fi : int + The file index that must be read from. + start : int + The start sample in the given file. + stop : int + The stop sample in the given file (inclusive). + cals : ndarray, shape (len(idx), 1) + Channel calibrations (already sub-indexed). + mult : ndarray, shape (len(idx), len(info['chs']) | None + The compensation + projection + cals matrix, if applicable. + """ + # Initial checks + start = int(start) + if stop is None or stop > self._raw_extras[fi]["n_samp"]: + stop = self._raw_extras[fi]["n_samp"] + + if start >= stop: + raise ValueError("No data in this range.") + + data_offset = self._raw_extras[fi]["start_data"] + data_size = 4 # sizeof(int) + n_channels = self._raw_extras[fi]["nchan"] + + data_left = (stop - start) * n_channels + + blocksize = ((int(100e6) // data_size) // n_channels) * n_channels + blocksize = min(data_left, blocksize) + + with open(self._filenames[fi], "rb") as fid: + # position file pointer + fid.seek(data_offset + data_size * start * n_channels, 0) + + for sample_start in np.arange(0, data_left, blocksize) // n_channels: + count = min(blocksize, data_left - sample_start * n_channels) + block = np.fromfile(fid, ">i4", count=count) + + if self._raw_extras[fi]["data_type"] == 4: + block = block.byteswap() + block = block.astype(np.float32) # convert to float32 + block = block.reshape(n_channels, -1, order="F") + + n_samples = block.shape[1] + sample_stop = sample_start + n_samples + + data_view = data[:, sample_start:sample_stop] + _mult_cal_one(data_view, block, idx, cals, mult) + + return data_view + + +@fill_doc +def read_raw_itab(fname, preload=False, verbose=None) -> RawITAB: + """Raw object from ITAB directory. + + Parameters + ---------- + fname : str + The raw file to load. Filename should end with *.raw + preload : bool or str (default False) + Preload data into memory for data manipulation and faster indexing. + If True, the data will be preloaded into memory (fast, requires + large amount of memory). If preload is a string, preload is the + file name of a memory-mapped file which is used to store the data + on the hard drive (slower, requires less memory). + %(verbose)s + + Returns + ------- + raw : instance of RawITAB + The raw data. + + See Also + -------- + mne.io.Raw : Documentation of attribute and methods. + + Notes + ----- + .. versionadded:: 0.01 + """ + return RawITAB(fname, preload=preload, verbose=verbose) diff --git a/mne/io/itab/mhd.py b/mne/io/itab/mhd.py new file mode 100755 index 00000000000..db0fa26896b --- /dev/null +++ b/mne/io/itab/mhd.py @@ -0,0 +1,215 @@ +"""Read .mhd file.""" + +# Author: Vittorio Pizzella +# +# License: BSD (3-clause) + +import numpy as np + +from ...utils import logger +from .constants import ITAB + + +def _read_double(fid, n=1): + """Read a double.""" + return np.fromfile(fid, "B", n) + + +# TODO: Is it similar to _read_char +def _read_ustring(fid, n_bytes): + """Read unsigned character string.""" + return np.fromfile(fid, ">B", n_bytes) + + +def _read_int2(fid): + """Read int from short.""" + return np.fromfile(fid, " +# simplified BSD-3 license + +from numpy.testing import assert_array_almost_equal + +from mne.datasets import testing +from mne.io import read_raw_fieldtrip, read_raw_itab +from mne.io.tests.test_raw import _test_raw_reader + +data_path = testing.data_path(download=False) + +mat_itab_fname = data_path / "itab" / "test_itab.mat" +raw_itab_fname = data_path / "itab" / "test_itab.raw" + + +# @testing.requires_testing_data +def test_itab_raw(): + """Test reading ITAB .raw files.""" + raw = read_raw_itab(raw_itab_fname, preload=True) + assert "RawITAB" in repr(raw) + + _test_raw_reader( + read_raw_itab, + fname=raw_itab_fname, + test_scaling=False, + ) + + test_ft_raw = read_raw_fieldtrip(mat_itab_fname, info=raw.info, data_name="data") + + itab_data = raw.get_data() + ft_data = test_ft_raw.get_data() + + assert_array_almost_equal(ft_data, itab_data) + + assert itab_data.shape == ft_data.shape + assert test_ft_raw.info["sfreq"] == raw.info["sfreq"] diff --git a/mne/io/tests/test_raw.py b/mne/io/tests/test_raw.py index 1d434fcf9b4..36705c3e47d 100644 --- a/mne/io/tests/test_raw.py +++ b/mne/io/tests/test_raw.py @@ -323,8 +323,8 @@ def _test_raw_reader( data_min = np.nanmin(full_data) # these limits could be relaxed if we actually find data with # huge values (in SI units) - assert data_max < 1e5 - assert data_min > -1e5 + assert data_max < 1e7 + assert data_min > -1e7 if isinstance(raw.info["dig"], list): for di, d in enumerate(raw.info["dig"]): assert isinstance(d, DigPoint), (di, d)