Merge pull request #260 from gkBCCN/spikegadgets_npix

NEW: probe reader for Neuropixels 1.0 in SpikeGadgets .rec file.

Author: alejoe91

src/probeinterface/__init__.py

@@ -16,6 +16,7 @@
     write_imro,
     read_BIDS_probe,
     write_BIDS_probe,
+    read_spikegadgets,
     read_spikeglx,
     parse_spikeglx_meta,
     get_saved_channel_indices_from_spikeglx_meta,

src/probeinterface/io.py

@@ -19,6 +19,7 @@
 from collections import OrderedDict
 from packaging.version import Version, parse
 import numpy as np
+from xml.etree import ElementTree
 
 from . import __version__
 from .probe import Probe
@@ -1263,12 +1264,194 @@ def write_imro(file: str | Path, probe: Probe):
         f.write("".join(ret))
 
 
+def read_spikegadgets(file: str | Path, raise_error: bool = True) -> ProbeGroup:
+    """
+    Find active channels of the given Neuropixels probe from a SpikeGadgets .rec file.
+    SpikeGadgets headstages support up to three Neuropixels 1.0 probes (as of March 28, 2024),
+    and information for all probes will be returned in a ProbeGroup object.
+
+
+    Parameters
+    ----------
+    file : Path or str
+        The .rec file path
+
+    Returns
+    -------
+    probe_group : ProbeGroup object
+
+    """
+    # ------------------------- #
+    #     Npix 1.0 constants    #
+    # ------------------------- #
+    TOTAL_NPIX_ELECTRODES = 960
+    MAX_ACTIVE_CHANNELS = 384
+    CONTACT_WIDTH = 16  # um
+    CONTACT_HEIGHT = 20  # um
+    # ------------------------- #
+
+    # Read the header and get Configuration elements
+    header_txt = parse_spikegadgets_header(file)
+    root = ElementTree.fromstring(header_txt)
+    hconf = root.find("HardwareConfiguration")
+    sconf = root.find("SpikeConfiguration")
+
+    # Get number of probes present (each has its own Device element)
+    probe_configs = [device for device in hconf if device.attrib["name"] == "NeuroPixels1"]
+    n_probes = len(probe_configs)
+
+    if n_probes == 0:
+        if raise_error:
+            raise Exception("No Neuropixels 1.0 probes found")
+        return None
+
+    # Container to store Probe objects
+    probe_group = ProbeGroup()
+
+    for curr_probe in range(1, n_probes + 1):
+        probe_config = probe_configs[curr_probe - 1]
+
+        # Get number of active channels from probe Device element
+        active_channel_str = [option for option in probe_config if option.attrib["name"] == "channelsOn"][0].attrib[
+            "data"
+        ]
+        active_channels = [int(ch) for ch in active_channel_str.split(" ") if ch]
+        n_active_channels = sum(active_channels)
+        assert len(active_channels) == TOTAL_NPIX_ELECTRODES
+        assert n_active_channels <= MAX_ACTIVE_CHANNELS
+
+        """
+        Within the SpikeConfiguration header element (sconf), there is a SpikeNTrode element
+        for each electrophysiology channel that contains information relevant to scaling and
+        otherwise displaying the information from that channel, as well as the id of the electrode
+        from which it is recording ('id').
+
+        Nested within each SpikeNTrode element is a SpikeChannel element with information about
+        the electrode dynamically connected to that channel.  This contains information relevant
+        for spike sorting, i.e., its spatial location along the probe shank and the hardware channel
+        to which it is connected.
+
+        Excerpt of a sample SpikeConfiguration element:
+
+        <SpikeConfiguration chanPerChip="1889715760" device="neuropixels1" categories="">
+            <SpikeNTrode viewLFPBand="0"
+                viewStimBand="0"
+                id="1384"  # @USE: The first digit is the probe number; the last three digits are the electrode number
+                lfpScalingToUv="0.018311105685598315"
+                LFPChan="1"
+                notchFreq="60"
+                rawRefOn="0"
+                refChan="1"
+                viewSpikeBand="1"
+                rawScalingToUv="0.018311105685598315"  # For Neuropixels 1.0, raw and spike scaling are identical
+                spikeScalingToUv="0.018311105685598315"  # Extracted when reading the raw data
+                refNTrodeID="1"
+                notchBW="10"
+                color="#c83200"
+                refGroup="2"
+                filterOn="1"
+                LFPHighFilter="200"
+                moduleDataOn="0"
+                groupRefOn="0"
+                lowFilter="600"
+                refOn="0"
+                notchFilterOn="0"
+                lfpRefOn="0"
+                lfpFilterOn="0"
+                highFilter="6000"
+            >
+                <SpikeChannel thresh="60"
+                    coord_dv="-480"  # @USE: dorsal-ventral coordinate in um (in pairs for staggered probe)
+                    spikeSortingGroup="1782505664"
+                    triggerOn="1"
+                    stimCapable="0"
+                    coord_ml="3192"  # @USE: medial-lateral coordinate in um
+                    coord_ap="3700"  # doesn't vary, assuming the shank's flat face is along the ML axis
+                    maxDisp="400"
+                    hwChan="735"  # @USE: unique device channel that is reading from electrode
+                />
+            </SpikeNTrode>
+            ...
+        </SpikeConfiguration>
+        """
+        # Find all channels/electrodes that belong to the current probe
+        contact_ids = []
+        device_channels = []
+        positions = np.zeros((n_active_channels, 2))
+
+        nt_i = 0  # Both probes are in sconf, so need an independent counter of probe electrodes while iterating through
+        for ntrode in sconf:
+            electrode_id = ntrode.attrib["id"]
+            if int(electrode_id[0]) == curr_probe:  # first digit of electrode id is probe number
+                contact_ids.append(electrode_id)
+                positions[nt_i, :] = (ntrode[0].attrib["coord_ml"], ntrode[0].attrib["coord_dv"])
+                device_channels.append(ntrode[0].attrib["hwChan"])
+                nt_i += 1
+        assert len(contact_ids) == n_active_channels
+
+        # Construct Probe object
+        probe = Probe(ndim=2, si_units="um", model_name="Neuropixels 1.0", manufacturer="IMEC")
+        probe.set_contacts(
+            contact_ids=contact_ids,
+            positions=positions,
+            shapes="square",
+            shank_ids=None,
+            shape_params={"width": CONTACT_WIDTH, "height": CONTACT_HEIGHT},
+        )
+
+        # Wire it (i.e., point contact/electrode ids to corresponding hardware/channel ids)
+        probe.set_device_channel_indices(device_channels)
+
+        # Create a nice polygon background when plotting the probes
+        x_min = positions[:, 0].min()
+        x_max = positions[:, 0].max()
+        x_mid = 0.5 * (x_max + x_min)
+        y_min = positions[:, 1].min()
+        y_max = positions[:, 1].max()
+        polygon_default = [
+            (x_min - 20, y_min - CONTACT_HEIGHT / 2),
+            (x_mid, y_min - 100),
+            (x_max + 20, y_min - CONTACT_HEIGHT / 2),
+            (x_max + 20, y_max + 20),
+            (x_min - 20, y_max + 20),
+        ]
+        probe.set_planar_contour(polygon_default)
+
+        # If there are multiple probes, they must be shifted such that they don't occupy the same coordinates.
+        probe.move([250 * (curr_probe - 1), 0])
+
+        # Add the probe to the probe container
+        probe_group.add_probe(probe)
+
+    return probe_group
+
+
+def parse_spikegadgets_header(file: str | Path) -> str:
+    """
+    Parse file (SpikeGadgets .rec format) into a string until "</Configuration>",
+    which is the last tag of the header, after which the binary data begins.
+    """
+    header_size = None
+    with open(file, mode="rb") as f:
+        while True:
+            line = f.readline()
+            if b"</Configuration>" in line:
+                header_size = f.tell()
+                break
+
+        if header_size is None:
+            ValueError("SpikeGadgets: the xml header does not contain '</Configuration>'")
+
+        f.seek(0)
+        return f.read(header_size).decode("utf8")
+
+
 def read_spikeglx(file: str | Path) -> Probe:
     """
     Read probe position for the meta file generated by SpikeGLX
 
     See http://billkarsh.github.io/SpikeGLX/#metadata-guides for implementation.
-    The x_pitch/y_pitch/width are set automatically depending the NP version.
+    The x_pitch/y_pitch/width are set automatically depending on the NP version.
 
     The shape is auto generated as a shank.
 
@@ -1333,7 +1516,7 @@ def read_spikeglx(file: str | Path) -> Probe:
 def parse_spikeglx_meta(meta_file: str | Path) -> dict:
     """
     Parse the "meta" file from spikeglx into a dict.
-    All fiields are kept in txt format and must also parsed themself.
+    All fields are kept in txt format and must also be parsed themselves.
     """
     meta_file = Path(meta_file)
     with meta_file.open(mode="r") as f: