mne-tools / mne-nirs

Process Near-Infrared Spectroscopy Data in MNE
https://mne.tools/mne-nirs/
BSD 3-Clause "New" or "Revised" License
79 stars 35 forks source link

read_snirf_aux_data() / fails when /nirs/aux/dataTimeSeries is 2D array #558

Open kdarti opened 2 months ago

kdarti commented 2 months ago

Describe the bug

The example aux dataset used in the aux tutorial has /nirs/aux/dataTimeSeries in 1D format.

This means that read_snirf_aux_data() runs fine on this example dataset.

However, the snirf spec states that /nirs/aux/dataTimeSeries should be a 2D array.

image

Trying to read a /nirs/aux/dataTimeSeries 2D array with read_snirf_aux_data() fails on line 56 due to d containing 2D arrays

https://github.com/mne-tools/mne-nirs/blob/2fdddfe9962ad000ebe60d71d0f9794d1800acbd/mne_nirs/io/snirf/_aux.py#L47-L57

image image

Steps to reproduce

Can't use the offical snirf samples for demonstration, because /nirs/aux/time in those files are 2D arrays, so this leads to another fail.

Use this file with the code below to reproduce .

import numpy as np
import logging
import h5py

from scipy import interpolate
from pandas import DataFrame
from mne.io import Raw

fname = r"04_25090.snirf"

dat = h5py.File(fname, 'r')

        if 'nirs' in dat:
            basename = "nirs"
        elif 'nirs1' in dat:
            basename = "nirs1"
        else:
            raise RuntimeError("Data does not contain nirs field")

        all_keys = list(dat.get(basename).keys())
        aux_keys = [i for i in all_keys if i.startswith('aux')]
        print(aux_keys)
        aux_names = [_decode_name(dat.get(f'{basename}/{k}/name'))
                     for k in aux_keys]
        logging.debug(f"Found auxiliary channels {aux_names}")

        d = {'times': raw.times}
        for idx, aux in enumerate(aux_keys):
            print("idx:",idx)
            print("aux:",aux)
            aux_data = np.array(dat.get(f'{basename}/{aux}/dataTimeSeries'))
            aux_time = np.array(dat.get(f'{basename}/{aux}/time'))
            aux_data_interp = interpolate.interp1d(aux_time, aux_data,
                                                   axis=0, bounds_error=False,
                                                   fill_value='extrapolate')
            aux_data_matched_to_raw = aux_data_interp(raw.times)
            d[aux_names[idx]] = aux_data_matched_to_raw

        df = DataFrame(data=d)
        df = df.set_index('times')

Expected results

Since /nirs/aux/dataTimeSeries should be 2D array it should work fine to read data stored in this way

Actual results

Fail due to code expecting a 1D array

Additional information

As a quick fix for my data I used .flatten() on the aux_data (line 48), but that's not really a fix in case the data is actually 2D.