diff --git a/doc/changes/dev/14060.bugfix.rst b/doc/changes/dev/14060.bugfix.rst new file mode 100644 index 00000000000..31baba279d0 --- /dev/null +++ b/doc/changes/dev/14060.bugfix.rst @@ -0,0 +1 @@ +Fix :func:`mne.io.read_raw_nihon` misreading large multichannel Nihon Kohden ``EEG-1200A`` recordings as a single ~1 second channel, by :newcontrib:`mgj05otf`. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 9c337b3d5c1..2c1030cd7c7 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -235,6 +235,7 @@ .. _Matti Toivonen: https://github.com/mattitoi .. _Mauricio Cespedes Tenorio: https://github.com/mcespedes99 .. _Melih Yayli: https://github.com/yaylim +.. _mgj05otf: https://github.com/mgj05otf .. _Michael Straube: https://github.com/mistraube .. _Michal Žák: https://github.com/michalrzak .. _Michiru Kaneda: https://github.com/rcmdnk diff --git a/mne/io/nihon/nihon.py b/mne/io/nihon/nihon.py index 2ea10c1d123..8df74385f5d 100644 --- a/mne/io/nihon/nihon.py +++ b/mne/io/nihon/nihon.py @@ -175,12 +175,45 @@ def _read_21e_file(fname): return _chan_labels +_nihon_1200a_version = "EEG-1200A V01.00" + + +def _read_nihon_1200a_ext_block(fid, _chan_labels): + # EEG-1200A (and newer) systems no longer store the real channel count, + # channel order, or start-of-data address in the data block itself (that + # header field is left over from older systems and is not meaningful for + # these files, which is why large multichannel EEG-1200A recordings were + # previously misread as a single ~1 second channel). Instead, that info + # lives in a chain of "extended blocks" reached via a separate pointer. + # Layout reverse-engineered by the Brainstorm project, see + # https://github.com/brainstorm-tools/brainstorm3/blob/master/toolbox/io/in_fopen_nk.m + fid.seek(0x03EE) + ext_address = np.fromfile(fid, np.uint32, 1)[0].astype(np.int64) + fid.seek(ext_address + 18) + extblock2_address = np.fromfile(fid, np.uint32, 1)[0].astype(np.int64) + fid.seek(extblock2_address + 20) + extblock3_address = np.fromfile(fid, np.uint32, 1)[0].astype(np.int64) + + fid.seek(extblock3_address + 68) + n_channels = np.fromfile(fid, np.uint16, 1)[0].astype(np.int64) + + channels = [] + for i_ch in range(n_channels): + fid.seek(extblock3_address + 72 + (i_ch * 10)) + t_idx = np.fromfile(fid, np.uint16, 1)[0] + channels.append(_chan_labels[t_idx]) + + data_start = extblock3_address + 72 + (n_channels * 10) + return dict(n_channels=n_channels, channels=channels, data_start=data_start) + + def _read_nihon_header(fname): # Read the Nihon Kohden EEG file header fname = _ensure_path(fname) _chan_labels = _read_21e_file(fname) header: dict[str, Any] = {} logger.info(f"Reading header from {fname}") + file_size = fname.stat().st_size with open(fname) as fid: version = np.fromfile(fid, "|S16", 1).astype("U16")[0] if version not in _valid_headers: @@ -199,6 +232,10 @@ def _read_nihon_header(fname): raise ValueError("Not a valid Nihon Kohden EEG file (waveform block)") header["version"] = version + ext_block = None + if version == _nihon_1200a_version: + ext_block = _read_nihon_1200a_ext_block(fid, _chan_labels) + fid.seek(0x0091) n_ctlblocks = np.fromfile(fid, np.uint8, 1)[0] header["n_ctlblocks"] = n_ctlblocks @@ -218,27 +255,51 @@ def _read_nihon_header(fname): t_data_address = np.fromfile(fid, np.uint32, 1)[0] t_datablock["address"] = t_data_address - fid.seek(t_data_address + 0x26) - t_n_channels = np.fromfile(fid, np.uint8, 1)[0].astype(np.int64) - t_datablock["n_channels"] = t_n_channels - - t_channels = [] - for i_ch in range(t_n_channels): - fid.seek(t_data_address + 0x27 + (i_ch * 10)) - t_idx = np.fromfile(fid, np.uint8, 1)[0] - t_channels.append(_chan_labels[t_idx]) - - t_datablock["channels"] = t_channels - - fid.seek(t_data_address + 0x1C) - t_record_duration = np.fromfile(fid, np.uint32, 1)[0].astype(np.int64) - t_datablock["duration"] = t_record_duration - fid.seek(t_data_address + 0x1A) sfreq = np.fromfile(fid, np.uint16, 1)[0] & 0x3FFF t_datablock["sfreq"] = sfreq.astype(np.int64) - t_datablock["n_samples"] = np.int64(t_record_duration * sfreq // 10) + if ext_block is not None: + if n_ctlblocks != 1 or n_datablocks != 1: + raise NotImplementedError( + "Reading EEG-1200A files with more than one " + "control block or more than one data block is " + "not supported." + ) + t_datablock["n_channels"] = ext_block["n_channels"] + t_datablock["channels"] = ext_block["channels"] + t_datablock["data_start"] = ext_block["data_start"] + # EEG-1200A data blocks don't store a usable record + # count, so the number of samples is inferred from how + # much data is left in the file after the header. + t_datablock["n_samples"] = np.int64( + (file_size - ext_block["data_start"]) + // (ext_block["n_channels"] + 1) + // 2 + ) + else: + fid.seek(t_data_address + 0x26) + t_n_channels = np.fromfile(fid, np.uint8, 1)[0].astype(np.int64) + t_datablock["n_channels"] = t_n_channels + + t_channels = [] + for i_ch in range(t_n_channels): + fid.seek(t_data_address + 0x27 + (i_ch * 10)) + t_idx = np.fromfile(fid, np.uint8, 1)[0] + t_channels.append(_chan_labels[t_idx]) + + t_datablock["channels"] = t_channels + t_datablock["data_start"] = ( + t_data_address + 0x27 + (t_n_channels * 10) + ) + + fid.seek(t_data_address + 0x1C) + t_record_duration = np.fromfile(fid, np.uint32, 1)[0].astype( + np.int64 + ) + t_datablock["duration"] = t_record_duration + + t_datablock["n_samples"] = np.int64(t_record_duration * sfreq // 10) t_controlblock["datablocks"].append(t_datablock) controlblocks.append(t_controlblock) header["controlblocks"] = controlblocks @@ -554,7 +615,7 @@ def _read_segment_file(self, data, idx, fi, start, stop, cals, mult): datablock = datablocks[start_block] n_channels = datablock["n_channels"] + 1 - datastart = datablock["address"] + 0x27 + (datablock["n_channels"] * 10) + datastart = datablock["data_start"] # Compute start offset based on the beginning of the block rel_start = start diff --git a/mne/io/nihon/tests/test_nihon.py b/mne/io/nihon/tests/test_nihon.py index 28996cfb49a..d8f0355d649 100644 --- a/mne/io/nihon/tests/test_nihon.py +++ b/mne/io/nihon/tests/test_nihon.py @@ -2,6 +2,7 @@ # License: BSD-3-Clause # Copyright the MNE-Python contributors. +import numpy as np import pytest from numpy.testing import assert_allclose @@ -9,6 +10,8 @@ from mne.io import read_raw_edf, read_raw_nihon from mne.io.nihon import nihon from mne.io.nihon.nihon import ( + _default_chan_labels, + _map_ch_to_specs, _read_nihon_annotations, _read_nihon_header, _read_nihon_metadata, @@ -18,6 +21,86 @@ data_path = testing.data_path(download=False) +def _write_uint(buf, offset, value, dtype): + buf[offset : offset + np.dtype(dtype).itemsize] = np.array( + value, dtype=dtype + ).tobytes() + + +def _write_str(buf, offset, s): + encoded = s.encode("ascii") + buf[offset : offset + len(encoded)] = encoded + + +def _write_synthetic_1200a_eeg(tmp_path, ch_indices, sfreq, samples): + """Write a minimal synthetic EEG-1200A file set (.EEG/.PNT/.21E/.LOG). + + EEG-1200A systems store the real channel count/order and the start of + the raw samples in a chain of "extended blocks" rather than in the + data block itself (see ``_read_nihon_1200a_ext_block``). This writes + just enough of that structure, with made-up (non-patient) sample data, + to exercise that code path without requiring a real recording. + """ + version = "EEG-1200A V01.00" + n_channels = len(ch_indices) + n_samples = samples.shape[1] + assert samples.shape == (n_channels + 1, n_samples) # +1 for marker channel + + ctl_address = 0x400 + # 0x17FE is a fixed, unconditional offset checked in every valid Nihon + # Kohden file (the "waveform block" signature byte), so the header must + # extend at least that far regardless of where the other blocks live. + data_address = 0x17FE + ext_address = 0x1900 + extblock2_address = 0x1A00 + extblock3_address = 0x1B00 + data_start = extblock3_address + 72 + n_channels * 10 + + buf = bytearray(data_start) + _write_str(buf, 0x0000, version) + _write_str(buf, 0x0081, version) + _write_uint(buf, 0x17FE, 1, np.uint8) # waveform sign + _write_uint(buf, 0x03EE, ext_address, np.uint32) + _write_uint(buf, 0x0091, 1, np.uint8) # n_ctlblocks + _write_uint(buf, 0x0092, ctl_address, np.uint32) + _write_uint(buf, ctl_address + 17, 1, np.uint8) # n_datablocks + _write_uint(buf, ctl_address + 18, data_address, np.uint32) + _write_uint(buf, data_address + 0x1A, sfreq, np.uint16) + + _write_uint(buf, ext_address + 18, extblock2_address, np.uint32) + _write_uint(buf, extblock2_address + 20, extblock3_address, np.uint32) + _write_uint(buf, extblock3_address + 68, n_channels, np.uint16) + for i_ch, idx in enumerate(ch_indices): + _write_uint(buf, extblock3_address + 72 + i_ch * 10, idx, np.uint16) + + # Samples are stored as the "NK int16" encoding: uint16 offset by 0x8000. + stored = (samples.view(np.uint16) + np.uint16(0x8000)).astype("1 data block should raise, not misread.""" + rng = np.random.RandomState(0) + ch_names = ["FP1", "FP2"] + ch_indices = [_default_chan_labels.index(name) for name in ch_names] + sfreq = 500 + samples = rng.randint(-1000, 1000, size=(len(ch_names) + 1, 10)).astype(np.int16) + + fname = _write_synthetic_1200a_eeg(tmp_path, ch_indices, sfreq, samples) + # Claim a second data block in the (only) control block. + with open(fname, "r+b") as fid: + ctl_address = 0x400 + fid.seek(ctl_address + 17) + fid.write(np.array(2, dtype=np.uint8).tobytes()) + + with pytest.raises(NotImplementedError, match="more than one"): + read_raw_nihon(fname)