Fix reading of large multichannel Nihon Kohden EEG-1200A files

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`.

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

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

mne/io/nihon/tests/test_nihon.py

@@ -2,13 +2,16 @@
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
+import numpy as np
 import pytest
 from numpy.testing import assert_allclose
 
 from mne.datasets import testing
 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("<u2")
+
+    fname = tmp_path / "SYNTH1200A.EEG"
+    with open(fname, "wb") as fid:
+        fid.write(bytes(buf))
+        fid.write(stored.tobytes(order="F"))
+
+    pnt = bytearray(0x40 + 14)
+    _write_str(pnt, 0x0000, version)
+    _write_str(pnt, 0x40, "20260101120000")
+    with open(fname.with_suffix(".PNT"), "wb") as fid:
+        fid.write(bytes(pnt))
+
+    log = bytearray(0x92)
+    _write_str(log, 0x0000, version)
+    _write_uint(log, 0x91, 0, np.uint8)  # n_logblocks
+    with open(fname.with_suffix(".LOG"), "wb") as fid:
+        fid.write(bytes(log))
+
+    with open(fname.with_suffix(".21E"), "w", encoding="ascii") as fid:
+        fid.write("[ELECTRODE]\n")
+        for idx in ch_indices:
+            fid.write(f"{idx}={_default_chan_labels[idx]}\n")
+
+    return fname
+
+
 @testing.requires_testing_data
 def test_nihon_eeg():
     """Test reading Nihon Kohden EEG files."""
@@ -77,7 +160,8 @@ def test_nihon_duplicate_channels(monkeypatch):
     fname = data_path / "NihonKohden" / "MB0400FU.EEG"
 
     def return_channel_duplicates(fname):
-        ch_names = nihon._default_chan_labels
+        # copy so we don't permanently mutate the shared module-level list
+        ch_names = list(nihon._default_chan_labels)
         ch_names[1] = ch_names[0]
         return ch_names
 
@@ -125,3 +209,52 @@ def test_nihon_calibration():
     assert raw.info["sfreq"] == raw_edf.info["sfreq"]
 
     assert_allclose(raw.get_data(), raw_edf.get_data())
+
+
+def test_nihon_1200a_extended_block(tmp_path):
+    """Test reading an EEG-1200A file via the extended channel block."""
+    rng = np.random.RandomState(0)
+    ch_names = ["FP1", "FP2", "F3", "F4"]
+    ch_indices = [_default_chan_labels.index(name) for name in ch_names]
+    sfreq = 500
+    n_samples = 200
+    # last row is the marker channel, which read_raw_nihon drops
+    samples = rng.randint(-1000, 1000, size=(len(ch_names) + 1, n_samples))
+    samples = samples.astype(np.int16)
+
+    fname = _write_synthetic_1200a_eeg(tmp_path, ch_indices, sfreq, samples)
+    raw = read_raw_nihon(fname, preload=True)
+
+    assert raw.ch_names == ch_names
+    assert raw.info["sfreq"] == sfreq
+    assert raw.n_times == n_samples
+
+    chan_labels_upper = [x.upper() for x in _default_chan_labels]
+    specs = [_map_ch_to_specs(name, chan_labels_upper) for name in ch_names]
+    expected = np.stack(
+        [
+            (samples[i].astype(np.float64) * spec["cal"] + spec["offset"])
+            * spec["unit"]
+            for i, spec in enumerate(specs)
+        ]
+    )
+    assert_allclose(raw.get_data(), expected)
+
+
+def test_nihon_1200a_multi_datablock_not_supported(tmp_path):
+    """EEG-1200A files with >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)