From 198408176f5a308764cb52d472066094cffbfefc Mon Sep 17 00:00:00 2001 From: Selma Lugtmeijer Date: Thu, 11 Jun 2026 16:12:53 +0200 Subject: [PATCH 01/11] Update io.py to include Siemens physlogs Function added to include physiological measures from a Siemens scanner. This takes into account different sampling rates for the different measures and that physiological might have been started manually so have different start/end times. By default is uses the dicom file to get the start trigger time of the aquisition. Co-Authored-By: Stefano Moia Co-Authored-By: Brandon Ingram --- phys2bids/io.py | 165 ++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 158 insertions(+), 7 deletions(-) diff --git a/phys2bids/io.py b/phys2bids/io.py index d4a5e3903..8196bce96 100644 --- a/phys2bids/io.py +++ b/phys2bids/io.py @@ -1,11 +1,13 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """phys2bids interfaces for loading extension files.""" - import logging +import os import warnings from copy import deepcopy +from glob import glob from itertools import groupby +from pathlib import Path import numpy as np @@ -485,10 +487,6 @@ def load_gep(filename): -------- physio_obj.BlueprintInput """ - import os - from glob import glob - from pathlib import Path - # Initiate lists of column names and units with time and trigger names = ["time", "trigger"] units = ["s", "mV"] # Assuming recording units are mV... @@ -510,7 +508,7 @@ def load_gep(filename): # Calculate time in seconds for first input (starts from -30s) interval = 1 / freq[0] duration = data[0].shape[0] * interval - t_ch = np.ogrid[-30 : duration - 30 : interval] + time_ch = np.ogrid[-30 : duration - 30 : interval] # Find and add additional data files filename = Path(filename) @@ -535,10 +533,163 @@ def load_gep(filename): trigger = np.hstack((np.zeros(int(30 / interval)), np.ones(int((duration - 30) / interval)))) # Create final list of timeseries - timeseries = [t_ch, trigger] + timeseries = [time_ch, trigger] timeseries.extend(data) return BlueprintInput(timeseries, freq, names, units, 1) +# SELMA +def acqtime_to_seconds(acq_time_str): + """ + """ + hh = int(acq_time_str[0:2]) + mm = int(acq_time_str[2:4]) + ss = int(acq_time_str[4:6]) + frac = float("0." + acq_time_str.split(".")[1]) if "." in acq_time_str else 0 + return hh*3600 + mm*60 + ss + frac + + +def load_siemens(filename, dicomfolder=None): + """ + Populate object phys_input from SIEMENS physiological files. + + Uses the filename that the user provides to find any matching inputs + from other recording types (.puls, .resp, or .ecg). + + **Note that the filename must not be altered from how it is output from + the scanner.** + + Populates physio_obj with all identified recording types (note that one + or more of these may not be true recordings as the scanner outputs all + possible types in all cases). The modality corresponding to the filename + entered by the user is put first (after time and trigger). + + Parameters + ---------- + filename: str + path to the SIEMENS scanner physiological file + + Returns + ------- + BlueprintInput + + Note + ---- + + SIEMENS physiological files do not record a trigger so a column is added at + position 1. This has a value of zero up to the scan start time and then + a value of one for the duration of the scan. + + See Also + -------- + physio_obj.BlueprintInput + """ + import fmri_physio_log as fpl + import pydicom + + # initialise lists of metadata and data + names = [] + units = [] + freq = [] + timeseries = [] + starttime = [] + stoptime = [] + + # Find and add additional data files + filename = Path(filename) + fnames = glob(os.path.join(filename.parent, f"*{filename.name}.*")) + if not len(fnames) == 0: + for fname in fnames: + if fname.endswith(".puls"): + names.append("PPG") + units.append("") + data = fpl.PhysioLog.from_filename(filename) + timeseries.append(data.ts) + freq.append(data.rate) + starttime.append(data.mdh.start_time) + stoptime.append(data.mdh.stop_time) + elif fname.endswith(".resp"): + names.append("respiratory") + units.append("") + data = fpl.PhysioLog.from_filename(filename) + timeseries.append(data.ts) + freq.append(data.rate) + starttime.append(data.mdh.start_time) + stoptime.append(data.mdh.stop_time) + elif fname.endswith(".ecg"): + names.append("ECG") + units.append("") + data = fpl.PhysioLog.from_filename(filename) + timeseries.append(data.ts) + freq.append(data.rate) + starttime.append(data.mdh.start_time) + stoptime.append(data.mdh.stop_time) + elif fname.endswith(".ext"): + names.append("EXT") + units.append("") + data = fpl.PhysioLog.from_filename(filename) + timeseries.append(data.ts) + freq.append(data.rate) + starttime.append(data.mdh.start_time) + stoptime.append(data.mdh.stop_time) + elif fname.endswith(".ext2"): + names.append("EXT2") + units.append("") + data = fpl.PhysioLog.from_filename(filename) + timeseries.append(data.ts) + freq.append(data.rate) + starttime.append(data.mdh.start_time) + stoptime.append(data.mdh.stop_time) + + checkstartlen = np.unique(starttime) + checkstoplen = np.unique(stoptime) + checkfreqlen = np.unique(freq) + + time_ch = np.ogrid[checkstartlen[0]: checkstoplen[-1]: 1/checkfreqlen[-1]*1000] + if checkstartlen.size > 1 or checkstoplen.size > 1: + for n, t in enumerate(timeseries): + time = np.ogrid[checkstartlen[0]: checkstoplen[-1]: 1/freq[n]*1000] + timeseries[n] = np.zeros(time) + phys_start = int(np.argmax(np.isclose(time, starttime[n]))) + phys_stop = phys_start + len(t) + timeseries[n][phys_start:phys_stop] = t + + # Calculate time in seconds + time_ch = time_ch / 1000 + + # Read the dicom to find start time, if it's given, otherwise assume it's + trigger = np.zeros_like(time_ch) + + try: + trigger_samples = np.asarray(data.trigger_times, dtype=int) + valid = trigger_samples[(trigger_samples >= 0) & (trigger_samples < trigger.size)] + trigger[valid] = 1.0 + LGR.info(f"Found {len(valid)} trigger events in log file.") + except AttributeError: + LGR.warning("Log file trigger events not found — attempting opening DICOM instead.") + + if dicomfolder is not None: + files = sorted([os.path.join(dicomfolder, f) for f in os.listdir(dicomfolder) if os.path.isfile(os.path.join(dicomfolder, f))]) + first_dcm = pydicom.dcmread(files[0], stop_before_pixels=True) + last_dcm = pydicom.dcmread(files[-1], stop_before_pixels=True) + + epi_start = acqtime_to_seconds(first_dcm.AcquisitionTime) - data.mdh.start_time / 1000 + epi_stop = acqtime_to_seconds(last_dcm.AcquisitionTime) - data.mdh.start_time / 1000 + + take_start = int(np.argmax(np.isclose(time_ch, epi_start))) + take_end = int(np.argmax(np.isclose(time_ch, epi_stop))) + + trigger[take_start:take_end] = 1 + else: + LGR.warning("DICOM folder not provided. Setting start_time to 0 (although " + "it might be wrong).") + + # Initiate lists of column names and units with time and trigger + names = ["time", "trigger"] + names + units = ["s", ""] + units + freq = [checkfreqlen[-1], checkfreqlen[-1]] + freq + + return BlueprintInput(timeseries, freq, names, units, 1) + def load_smr(filename, chtrig=0): """Load Spike2 smr file and populate object phys_input. From 3d1f3e4f6218fc2b007a8e9a2613079b6c692147 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 11 Jun 2026 14:41:02 +0000 Subject: [PATCH 02/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- phys2bids/io.py | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/phys2bids/io.py b/phys2bids/io.py index 8196bce96..561c331d5 100644 --- a/phys2bids/io.py +++ b/phys2bids/io.py @@ -537,15 +537,15 @@ def load_gep(filename): timeseries.extend(data) return BlueprintInput(timeseries, freq, names, units, 1) + # SELMA def acqtime_to_seconds(acq_time_str): - """ - """ + """ """ hh = int(acq_time_str[0:2]) mm = int(acq_time_str[2:4]) ss = int(acq_time_str[4:6]) frac = float("0." + acq_time_str.split(".")[1]) if "." in acq_time_str else 0 - return hh*3600 + mm*60 + ss + frac + return hh * 3600 + mm * 60 + ss + frac def load_siemens(filename, dicomfolder=None): @@ -599,7 +599,7 @@ def load_siemens(filename, dicomfolder=None): fnames = glob(os.path.join(filename.parent, f"*{filename.name}.*")) if not len(fnames) == 0: for fname in fnames: - if fname.endswith(".puls"): + if fname.endswith(".puls"): names.append("PPG") units.append("") data = fpl.PhysioLog.from_filename(filename) @@ -644,10 +644,10 @@ def load_siemens(filename, dicomfolder=None): checkstoplen = np.unique(stoptime) checkfreqlen = np.unique(freq) - time_ch = np.ogrid[checkstartlen[0]: checkstoplen[-1]: 1/checkfreqlen[-1]*1000] + time_ch = np.ogrid[checkstartlen[0] : checkstoplen[-1] : 1 / checkfreqlen[-1] * 1000] if checkstartlen.size > 1 or checkstoplen.size > 1: for n, t in enumerate(timeseries): - time = np.ogrid[checkstartlen[0]: checkstoplen[-1]: 1/freq[n]*1000] + time = np.ogrid[checkstartlen[0] : checkstoplen[-1] : 1 / freq[n] * 1000] timeseries[n] = np.zeros(time) phys_start = int(np.argmax(np.isclose(time, starttime[n]))) phys_stop = phys_start + len(t) @@ -656,7 +656,7 @@ def load_siemens(filename, dicomfolder=None): # Calculate time in seconds time_ch = time_ch / 1000 - # Read the dicom to find start time, if it's given, otherwise assume it's + # Read the dicom to find start time, if it's given, otherwise assume it's trigger = np.zeros_like(time_ch) try: @@ -668,7 +668,13 @@ def load_siemens(filename, dicomfolder=None): LGR.warning("Log file trigger events not found — attempting opening DICOM instead.") if dicomfolder is not None: - files = sorted([os.path.join(dicomfolder, f) for f in os.listdir(dicomfolder) if os.path.isfile(os.path.join(dicomfolder, f))]) + files = sorted( + [ + os.path.join(dicomfolder, f) + for f in os.listdir(dicomfolder) + if os.path.isfile(os.path.join(dicomfolder, f)) + ] + ) first_dcm = pydicom.dcmread(files[0], stop_before_pixels=True) last_dcm = pydicom.dcmread(files[-1], stop_before_pixels=True) @@ -680,8 +686,10 @@ def load_siemens(filename, dicomfolder=None): trigger[take_start:take_end] = 1 else: - LGR.warning("DICOM folder not provided. Setting start_time to 0 (although " - "it might be wrong).") + LGR.warning( + "DICOM folder not provided. Setting start_time to 0 (although " + "it might be wrong)." + ) # Initiate lists of column names and units with time and trigger names = ["time", "trigger"] + names From 213ea19a0dd14eb169862c6d3624cb0dc94d7b37 Mon Sep 17 00:00:00 2001 From: Selma Lugtmeijer Date: Thu, 11 Jun 2026 17:19:39 +0200 Subject: [PATCH 03/11] Added docstrings specific for Siemens --- phys2bids/io.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/phys2bids/io.py b/phys2bids/io.py index 561c331d5..388af3cc9 100644 --- a/phys2bids/io.py +++ b/phys2bids/io.py @@ -555,18 +555,17 @@ def load_siemens(filename, dicomfolder=None): Uses the filename that the user provides to find any matching inputs from other recording types (.puls, .resp, or .ecg). - **Note that the filename must not be altered from how it is output from - the scanner.** - - Populates physio_obj with all identified recording types (note that one - or more of these may not be true recordings as the scanner outputs all - possible types in all cases). The modality corresponding to the filename - entered by the user is put first (after time and trigger). + Populates physio_obj with all identified recording types. + Takes into account possible differences in start and stop times across channels. + Uses DICOM metadata to find scan start and stop times and create trigger channel. + The modalities are in the order of 'names' entered by the user (after time and trigger). Parameters ---------- filename: str path to the SIEMENS scanner physiological file + dicomfolder: str or None + path to the folder containing the DICOM files Returns ------- From 779c16332982c9971b1ebcfaf2b9026f58be6857 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 11 Jun 2026 15:19:54 +0000 Subject: [PATCH 04/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- phys2bids/io.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/phys2bids/io.py b/phys2bids/io.py index 388af3cc9..c3cde1479 100644 --- a/phys2bids/io.py +++ b/phys2bids/io.py @@ -555,9 +555,9 @@ def load_siemens(filename, dicomfolder=None): Uses the filename that the user provides to find any matching inputs from other recording types (.puls, .resp, or .ecg). - Populates physio_obj with all identified recording types. + Populates physio_obj with all identified recording types. Takes into account possible differences in start and stop times across channels. - Uses DICOM metadata to find scan start and stop times and create trigger channel. + Uses DICOM metadata to find scan start and stop times and create trigger channel. The modalities are in the order of 'names' entered by the user (after time and trigger). Parameters From ff80c21c4d73c323d0179d215ca896aae53bbff8 Mon Sep 17 00:00:00 2001 From: Selma Lugtmeijer Date: Thu, 11 Jun 2026 17:35:55 +0200 Subject: [PATCH 05/11] Deleted my name --- phys2bids/io.py | 1 - setup.cfg | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/phys2bids/io.py b/phys2bids/io.py index c3cde1479..dfa09f941 100644 --- a/phys2bids/io.py +++ b/phys2bids/io.py @@ -538,7 +538,6 @@ def load_gep(filename): return BlueprintInput(timeseries, freq, names, units, 1) -# SELMA def acqtime_to_seconds(acq_time_str): """ """ hh = int(acq_time_str[0:2]) diff --git a/setup.cfg b/setup.cfg index 20204cd93..2216a142a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -126,6 +126,7 @@ skip = venvs,.venv,versioneer.py,.git,build,./docs/_build write-changes = count = quiet-level = 3 +ignore-words-list = puls [tool:pytest] doctest_optionflags = NORMALIZE_WHITESPACE From b891ad676697a68cff3b0fbe8851865ffd80ec22 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 11 Jun 2026 17:44:17 +0200 Subject: [PATCH 06/11] Add dicomdir input to main workflow and CLI --- phys2bids/cli/run.py | 8 ++++++++ phys2bids/phys2bids.py | 7 +++++++ setup.cfg | 1 + 3 files changed, 16 insertions(+) diff --git a/phys2bids/cli/run.py b/phys2bids/cli/run.py index b1a3dec14..55f007f1d 100644 --- a/phys2bids/cli/run.py +++ b/phys2bids/cli/run.py @@ -164,6 +164,14 @@ def _get_parser(): help="Column header (for json file output).", default=[], ) + optional.add_argument( + "-ddir", + "--dicomdir", + dest="dicomdir", + type=str, + help="full path to MRI dicom folder corresponding to physiological data. Use this in conjunction with SIEMENS log files ONLY", + default="", + ) optional.add_argument( "-yml", "--participant-yml", diff --git a/phys2bids/phys2bids.py b/phys2bids/phys2bids.py index 5161792bf..a8632767e 100644 --- a/phys2bids/phys2bids.py +++ b/phys2bids/phys2bids.py @@ -144,6 +144,7 @@ def phys2bids( thr=None, pad=9, ch_name=[], + dicomdir=None, yml="", debug=False, quiet=False, @@ -260,6 +261,12 @@ def phys2bids( from phys2bids.io import load_gep phys_in = load_gep(infile) + elif ftype in [ + "puls", + ]: + from phys2bids.io import load_siemens + + phys_in = load_siemens(infile, dicomdir) LGR.info("Checking that units of measure are BIDS compatible") for index, unit in enumerate(phys_in.units): diff --git a/setup.cfg b/setup.cfg index 20204cd93..2216a142a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -126,6 +126,7 @@ skip = venvs,.venv,versioneer.py,.git,build,./docs/_build write-changes = count = quiet-level = 3 +ignore-words-list = puls [tool:pytest] doctest_optionflags = NORMALIZE_WHITESPACE From f88287788e01493e3a6c9716a6d73f4011067b39 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 11 Jun 2026 17:44:37 +0200 Subject: [PATCH 07/11] Add siemens python modules --- setup.cfg | 3 +++ 1 file changed, 3 insertions(+) diff --git a/setup.cfg b/setup.cfg index 2216a142a..ccf1e0a02 100644 --- a/setup.cfg +++ b/setup.cfg @@ -49,6 +49,9 @@ doc = sphinx >=2.0 sphinx-argparse sphinx_rtd_theme +siemens = + fmri_physio_log + pydicom style = flake8 >=4.0 black From 42f66002405a8c303f9f44e9edbcf645935797cf Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 11 Jun 2026 17:50:36 +0200 Subject: [PATCH 08/11] Add filetypes supported --- phys2bids/phys2bids.py | 4 +--- phys2bids/utils.py | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/phys2bids/phys2bids.py b/phys2bids/phys2bids.py index a8632767e..492fa7337 100644 --- a/phys2bids/phys2bids.py +++ b/phys2bids/phys2bids.py @@ -261,9 +261,7 @@ def phys2bids( from phys2bids.io import load_gep phys_in = load_gep(infile) - elif ftype in [ - "puls", - ]: + elif ftype in ["puls", "resp", "ecg", "ext", "ext2"]: from phys2bids.io import load_siemens phys_in = load_siemens(infile, dicomdir) diff --git a/phys2bids/utils.py b/phys2bids/utils.py index 05949afa9..0b38fa85d 100644 --- a/phys2bids/utils.py +++ b/phys2bids/utils.py @@ -9,7 +9,7 @@ LGR = logging.getLogger(__name__) -SUPPORTED_FTYPES = ("acq", "txt", "mat", "gep") +SUPPORTED_FTYPES = ("acq", "txt", "mat", "gep", "puls", "resp", "ecg", "ext", "ext2") def check_input_ext(filename, ext): From d33c0260fc15331637553c4c42a5179dd8e8d97f Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Fri, 12 Jun 2026 11:56:55 +0200 Subject: [PATCH 09/11] Start adding tests for SIEMENS Co-Authored-By: Selma Lugtmeijer <134110542+slugtmeijer@users.noreply.github.com> Co-Authored-By: Carsten Schmidt-Samoa <51449795+CarstenSchmidt-Samoa@users.noreply.github.com> --- phys2bids/tests/conftest.py | 19 ++++++++++++++----- phys2bids/tests/test_io.py | 19 +++++++++++++++++++ 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/phys2bids/tests/conftest.py b/phys2bids/tests/conftest.py index 0c43e6dde..92c0ff8c5 100644 --- a/phys2bids/tests/conftest.py +++ b/phys2bids/tests/conftest.py @@ -116,13 +116,13 @@ def ge_one_gep_file(testpath): @pytest.fixture def ge_two_gep_files_ppg(testpath): - tmp = fetch_file("qawjv", testpath, "RESPData_epiRT_0000000000_00_00_000.gep") + _ = fetch_file("qawjv", testpath, "RESPData_epiRT_0000000000_00_00_000.gep") return fetch_file("wb84d", testpath, "PPGData_epiRT_0000000000_00_00_000.gep") @pytest.fixture def ge_two_gep_files_resp(testpath): - tmp = fetch_file("wb84d", testpath, "PPGData_epiRT_0000000000_00_00_000.gep") + _ = fetch_file("wb84d", testpath, "PPGData_epiRT_0000000000_00_00_000.gep") return fetch_file("qawjv", testpath, "RESPData_epiRT_0000000000_00_00_000.gep") @@ -133,17 +133,26 @@ def ge_one_raw_file(testpath): @pytest.fixture def ge_two_raw_files(testpath): - tmp = fetch_file("49xpw", testpath, "RESPData_epiRT_0000000000_00_00_000") + _ = fetch_file("49xpw", testpath, "RESPData_epiRT_0000000000_00_00_000") return fetch_file("u9wsr", testpath, "PPGData_epiRT_0000000000_00_00_000") @pytest.fixture def ge_badfiles(testpath): - tmp = fetch_file("tdmyn", testpath, "PPGData_epiRT_columnscsv_00_00_000") - tmp = fetch_file("b6skq", testpath, "PPGData_epiRT_columnstsv_00_00_000") + _ = fetch_file("tdmyn", testpath, "PPGData_epiRT_columnscsv_00_00_000") + _ = fetch_file("b6skq", testpath, "PPGData_epiRT_columnstsv_00_00_000") return fetch_file("8235b", testpath, "PPGData_epiRT_string0000_00_00_000") +@pytest.fixture +def SIEMENS_files(testpath): + _ = fetch_file(" ", testpath, "275_pulse_part_1.puls") + _ = fetch_file(" ", testpath, "275_pulse_part_1.ecg") + _ = fetch_file(" ", testpath, "275_pulse_part_1.ext") + _ = fetch_file(" ", testpath, "275_pulse_part_1.ext2") + return fetch_file(" ", testpath, "275_resp_part_1.resp") + + @pytest.fixture def spike2_smrx_file(testpath): return fetch_file("7x5qw", testpath, "Test_ppg_pulse_spike2.smrx") diff --git a/phys2bids/tests/test_io.py b/phys2bids/tests/test_io.py index ffed64546..1390f584c 100644 --- a/phys2bids/tests/test_io.py +++ b/phys2bids/tests/test_io.py @@ -2,6 +2,7 @@ import os import sys +import fmri_physio_log as fpl import numpy as np import pytest from pytest import raises @@ -217,6 +218,24 @@ def test_load_gep_two_files_resp(ge_two_gep_files_resp, testpath): assert np.array_equal(gep_data2, phys_obj.timeseries[3]) +def test_load_SIEMENS_two_files_ppg(SIEMENS_files, testpath): + # Load data + phys_obj = io.load_siemens(SIEMENS_files) + + assert phys_obj.ch_name == ["time", "trigger", "ECG", "PPG", "respiratory"] + + # Check the channel data is as expected + siemens_ecg = fpl.PhysioLog.from_filename(os.path.join(testpath, "275_pulse_part_1.ecg")) + siemens_resp = fpl.PhysioLog.from_filename(SIEMENS_files) + siemens_pulse = fpl.PhysioLog.from_filename(os.path.join(testpath, "275_pulse_part_1.puls")) + assert np.array_equal(siemens_ecg, phys_obj.timeseries[2]) + assert np.array_equal(siemens_pulse, phys_obj.timeseries[3]) + assert np.array_equal(siemens_resp, phys_obj.timeseries[4]) + + +# check if data in phys_obj timeseries are padded, padding needs to be undone for the test to succeed + + @pytest.mark.skipif( sys.version_info < (3, 7) or sys.version_info > (3, 9), reason="Requires python between 3.7 and 3.9", From 28d93bd648b716c8f5f69be70fe9f9ab54087b35 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Fri, 12 Jun 2026 12:09:24 +0200 Subject: [PATCH 10/11] Change saving of file to object to deal with NaNs --- phys2bids/phys2bids.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/phys2bids/phys2bids.py b/phys2bids/phys2bids.py index 492fa7337..ed9ae99ae 100644 --- a/phys2bids/phys2bids.py +++ b/phys2bids/phys2bids.py @@ -520,6 +520,10 @@ def phys2bids( phys_out[key].filename = f"{phys_out[key].filename}_{uniq_freq:.0f}Hz" LGR.info(f"Exporting files for take {take} freq {uniq_freq}") + # Check if timeseries has NaNs and save them as n/a + phys_out[key].timeseries = phys_out[key].timeseries.astype(object) + phys_out[key].timeseries[np.isnan(phys_out[key].timeseries)] = "n/a" + np.savetxt( phys_out[key].filename + ".tsv.gz", phys_out[key].timeseries, From 3e2724f7c17202de01419cd95b52117c0fc8ba80 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Sat, 13 Jun 2026 16:54:57 +0200 Subject: [PATCH 11/11] Finishing all the issues and adding a little testing Co-Authored-By: Selma Lugtmeijer <134110542+slugtmeijer@users.noreply.github.com> Co-Authored-By: Carsten Schmidt-Samoa <51449795+CarstenSchmidt-Samoa@users.noreply.github.com> --- phys2bids/io.py | 148 ++++++++++++++++++++++++++----------- phys2bids/tests/test_io.py | 24 ++++-- 2 files changed, 125 insertions(+), 47 deletions(-) diff --git a/phys2bids/io.py b/phys2bids/io.py index dfa09f941..730886c4e 100644 --- a/phys2bids/io.py +++ b/phys2bids/io.py @@ -547,6 +547,81 @@ def acqtime_to_seconds(acq_time_str): return hh * 3600 + mm * 60 + ss + frac +def physiologtime_to_seconds(pulse_time): + """ """ + hh = pulse_time.hour + mm = pulse_time.minute + ss = pulse_time.second + frac = pulse_time.microsecond / 1000000 + return hh * 3600 + mm * 60 + ss + frac + + +def read_siemens_channel( + filename, channel_type, names, units, timeseries, freq, starttime, stoptime +): + """ """ + import fmri_physio_log as fpl + + data = fpl.PhysioLog.from_filename(filename) + + if np.all(data.ts != 0): + + if channel_type == "ECG": + freq_mod = round( + len(data.ts[::4]) + / ( + ( + physiologtime_to_seconds(data.mdh.stop_time) + - physiologtime_to_seconds(data.mdh.start_time) + ) + * data.rate + ) + ) + if freq_mod != 1: + LGR.warning( + f"Incoherency found between logged time and frequency by a factor of {freq_mod}. Please check frequency time." + ) + + frequency = data.freq * freq_mod + start = physiologtime_to_seconds(data.mdh.start_time) + stop = start + len(data.ts - 1) * frequency + for i in range(4): + names.append(f"{channel_type}{i+1}") + units.append("") + timeseries.append(data.ts[i::4]) + freq.append(data.rate * freq_mod) + starttime.append(physiologtime_to_seconds(data.mdh.start_time)) + stoptime.append(stop) + + else: + freq_mod = round( + len(data.ts[::4]) + / ( + ( + physiologtime_to_seconds(data.mdh.stop_time) + - physiologtime_to_seconds(data.mdh.start_time) + ) + * data.rate + ) + ) + if freq_mod != 1: + LGR.warning( + f"Incoherency found between logged time and frequency by a factor of {freq_mod}. Please check frequency time." + ) + + frequency = data.freq * freq_mod + start = physiologtime_to_seconds(data.mdh.start_time) + stop = start + len(data.ts - 1) * frequency + names.append(channel_type) + units.append("") + timeseries.append(data.ts) + freq.append(data.rate * freq_mod) + starttime.append(physiologtime_to_seconds(data.mdh.start_time)) + stoptime.append(stop) + + return names, units, timeseries, freq, starttime, stoptime, data + + def load_siemens(filename, dicomfolder=None): """ Populate object phys_input from SIEMENS physiological files. @@ -581,7 +656,6 @@ def load_siemens(filename, dicomfolder=None): -------- physio_obj.BlueprintInput """ - import fmri_physio_log as fpl import pydicom # initialise lists of metadata and data @@ -594,49 +668,33 @@ def load_siemens(filename, dicomfolder=None): # Find and add additional data files filename = Path(filename) - fnames = glob(os.path.join(filename.parent, f"*{filename.name}.*")) + fnames = sorted(glob(os.path.join(filename.parent, f"*{filename.name}.*"))) if not len(fnames) == 0: for fname in fnames: - if fname.endswith(".puls"): - names.append("PPG") - units.append("") - data = fpl.PhysioLog.from_filename(filename) - timeseries.append(data.ts) - freq.append(data.rate) - starttime.append(data.mdh.start_time) - stoptime.append(data.mdh.stop_time) - elif fname.endswith(".resp"): - names.append("respiratory") - units.append("") - data = fpl.PhysioLog.from_filename(filename) - timeseries.append(data.ts) - freq.append(data.rate) - starttime.append(data.mdh.start_time) - stoptime.append(data.mdh.stop_time) - elif fname.endswith(".ecg"): - names.append("ECG") - units.append("") - data = fpl.PhysioLog.from_filename(filename) - timeseries.append(data.ts) - freq.append(data.rate) - starttime.append(data.mdh.start_time) - stoptime.append(data.mdh.stop_time) + if fname.endswith(".ecg"): + names, units, timeseries, freq, starttime, stoptime, data = read_siemens_channel( + fname, "ECG", names, units, timeseries, freq, starttime, stoptime + ) + elif fname.endswith(".ext"): - names.append("EXT") - units.append("") - data = fpl.PhysioLog.from_filename(filename) - timeseries.append(data.ts) - freq.append(data.rate) - starttime.append(data.mdh.start_time) - stoptime.append(data.mdh.stop_time) + names, units, timeseries, freq, starttime, stoptime, data = read_siemens_channel( + fname, "EXT", names, units, timeseries, freq, starttime, stoptime + ) + elif fname.endswith(".ext2"): - names.append("EXT2") - units.append("") - data = fpl.PhysioLog.from_filename(filename) - timeseries.append(data.ts) - freq.append(data.rate) - starttime.append(data.mdh.start_time) - stoptime.append(data.mdh.stop_time) + names, units, timeseries, freq, starttime, stoptime, data = read_siemens_channel( + fname, "EXT2", names, units, timeseries, freq, starttime, stoptime + ) + + elif fname.endswith(".puls"): + names, units, timeseries, freq, starttime, stoptime, data = read_siemens_channel( + fname, "PPG", names, units, timeseries, freq, starttime, stoptime + ) + + elif fname.endswith(".resp"): + names, units, timeseries, freq, starttime, stoptime, data = read_siemens_channel( + fname, "respiratory", names, units, timeseries, freq, starttime, stoptime + ) checkstartlen = np.unique(starttime) checkstoplen = np.unique(stoptime) @@ -646,7 +704,7 @@ def load_siemens(filename, dicomfolder=None): if checkstartlen.size > 1 or checkstoplen.size > 1: for n, t in enumerate(timeseries): time = np.ogrid[checkstartlen[0] : checkstoplen[-1] : 1 / freq[n] * 1000] - timeseries[n] = np.zeros(time) + timeseries[n] = np.ones_like(time) * np.nan phys_start = int(np.argmax(np.isclose(time, starttime[n]))) phys_stop = phys_start + len(t) timeseries[n][phys_start:phys_stop] = t @@ -677,7 +735,13 @@ def load_siemens(filename, dicomfolder=None): last_dcm = pydicom.dcmread(files[-1], stop_before_pixels=True) epi_start = acqtime_to_seconds(first_dcm.AcquisitionTime) - data.mdh.start_time / 1000 - epi_stop = acqtime_to_seconds(last_dcm.AcquisitionTime) - data.mdh.start_time / 1000 + epi_stop = ( + acqtime_to_seconds(last_dcm.AcquisitionTime) + - data.mdh.start_time / 1000 + + float(last_dcm.get("RepetitionTime", 0)) / 1000 + ) + + # epi_stop needs to be incremented by TR in seconds take_start = int(np.argmax(np.isclose(time_ch, epi_start))) take_end = int(np.argmax(np.isclose(time_ch, epi_stop))) diff --git a/phys2bids/tests/test_io.py b/phys2bids/tests/test_io.py index 1390f584c..a3299bf9e 100644 --- a/phys2bids/tests/test_io.py +++ b/phys2bids/tests/test_io.py @@ -218,19 +218,33 @@ def test_load_gep_two_files_resp(ge_two_gep_files_resp, testpath): assert np.array_equal(gep_data2, phys_obj.timeseries[3]) -def test_load_SIEMENS_two_files_ppg(SIEMENS_files, testpath): +def test_SIEMENS(SIEMENS_files, testpath): # Load data phys_obj = io.load_siemens(SIEMENS_files) - assert phys_obj.ch_name == ["time", "trigger", "ECG", "PPG", "respiratory"] + assert phys_obj.ch_name == [ + "time", + "trigger", + "ECG1", + "ECG2", + "ECG3", + "ECG4", + "PPG", + "respiratory", + ] # Check the channel data is as expected siemens_ecg = fpl.PhysioLog.from_filename(os.path.join(testpath, "275_pulse_part_1.ecg")) + data_ecg1 = np.array(siemens_ecg.ts[::4]) + data_ecg3 = np.array(siemens_ecg.ts[2::4]) siemens_resp = fpl.PhysioLog.from_filename(SIEMENS_files) + data_resp = np.array(siemens_resp.ts) siemens_pulse = fpl.PhysioLog.from_filename(os.path.join(testpath, "275_pulse_part_1.puls")) - assert np.array_equal(siemens_ecg, phys_obj.timeseries[2]) - assert np.array_equal(siemens_pulse, phys_obj.timeseries[3]) - assert np.array_equal(siemens_resp, phys_obj.timeseries[4]) + data_pulse = np.array(siemens_pulse.ts) + assert np.array_equal(data_ecg1, phys_obj.timeseries[2][~np.isnan(phys_obj.timeseries[2])]) + assert np.array_equal(data_ecg3, phys_obj.timeseries[4][~np.isnan(phys_obj.timeseries[4])]) + assert np.array_equal(data_pulse, phys_obj.timeseries[6][~np.isnan(phys_obj.timeseries[6])]) + assert np.array_equal(data_resp, phys_obj.timeseries[7][~np.isnan(phys_obj.timeseries[7])]) # check if data in phys_obj timeseries are padded, padding needs to be undone for the test to succeed