From 6734be4a117fc3a992029723b4613559a48b2573 Mon Sep 17 00:00:00 2001 From: Samuel Garcia Date: Tue, 30 Jan 2024 16:32:20 +0100 Subject: [PATCH 1/3] Implement respiration CO2 cycle detection. --- physio/parameters.py | 19 ++- physio/respiration.py | 257 ++++++++++++++++++++++++++++---------- tests/test_respiration.py | 4 +- 3 files changed, 207 insertions(+), 73 deletions(-) diff --git a/physio/parameters.py b/physio/parameters.py index fdea5db..2c4118a 100644 --- a/physio/parameters.py +++ b/physio/parameters.py @@ -40,18 +40,30 @@ def get_ecg_parameters(parameter_preset): _resp_parameters['human_airflow'] = dict( preprocess=dict(band=7., btype='lowpass', ftype='bessel', order=5, normalize=False), smooth=dict(win_shape='gaussian', sigma_ms=60.0), - cycle_detection=dict(inspiration_adjust_on_derivative=False), + cycle_detection=dict(method="crossing_baseline", inspiration_adjust_on_derivative=False), baseline=dict(baseline_mode='median'), + features=dict(compute_volume=True, compute_amplitude=True), cycle_clean=dict(low_limit_log_ratio=4.), ) +_resp_parameters['human_co2'] = dict( + preprocess=dict(band=10., btype='lowpass', ftype='bessel', order=5, normalize=False), + smooth=dict(win_shape='gaussian', sigma_ms=40.0), + cycle_detection=dict(method="co2", thresh_inspi_factor=0.08, thresh_expi_factor=0.05, clean_by_mid_value=True), + baseline=dict(baseline_mode='median'), + features=dict(compute_volume=False, compute_amplitude=False), + cycle_clean=None, # no clean +) + + _resp_parameters['rat_plethysmo'] = dict( preprocess=dict(band=30., btype='lowpass', ftype='bessel', order=5, normalize=False), smooth=dict(win_shape='gaussian', sigma_ms=5.0), #~ smooth=None, - cycle_detection=dict(inspiration_adjust_on_derivative=False), + cycle_detection=dict(method="crossing_baseline", inspiration_adjust_on_derivative=False), baseline=dict(baseline_mode='manual', baseline=0.), + features=dict(compute_volume=True, compute_amplitude=True), cycle_clean=dict(low_limit_log_ratio=5.), ) @@ -60,8 +72,9 @@ def get_ecg_parameters(parameter_preset): preprocess=dict(band=30., btype='lowpass', ftype='bessel', order=5, normalize=False), smooth=dict(win_shape='gaussian', sigma_ms=5.0), #~ smooth=None, - cycle_detection=dict(epsilon_factor1=30, epsilon_factor2=10, inspiration_adjust_on_derivative=False), + cycle_detection=dict(method="crossing_baseline", epsilon_factor1=30, epsilon_factor2=10, inspiration_adjust_on_derivative=False), baseline=dict(baseline_mode='manual', baseline=0.), + features=dict(compute_volume=True, compute_amplitude=True), cycle_clean=dict(low_limit_log_ratio=5.), ) diff --git a/physio/respiration.py b/physio/respiration.py index cd01d5c..b2e39b4 100644 --- a/physio/respiration.py +++ b/physio/respiration.py @@ -13,8 +13,6 @@ def compute_respiration(raw_resp, srate, parameter_preset='human_airflow', param * clean cycles * compute metrics cycle by cycle - - Parameters ---------- raw_resp: np.array @@ -51,16 +49,27 @@ def compute_respiration(raw_resp, srate, parameter_preset='human_airflow', param resp = smooth_signal(resp, srate, **params['smooth']) resp += center + + + baseline = get_respiration_baseline(resp, srate, **params['baseline']) - baseline_detect = baseline + + detection_method = params['cycle_detection']["method"] + + if detection_method == "crossing_baseline": + cycles = detect_respiration_cycles(resp, srate, + baseline_mode='manual', baseline=baseline, + **params['cycle_detection']) + else: + cycles = detect_respiration_cycles(resp, srate, + **params['cycle_detection']) - cycles = detect_respiration_cycles(resp, srate, baseline_mode='manual', baseline=baseline, - **params['cycle_detection']) - resp_cycles = compute_respiration_cycle_features(resp, srate, cycles, baseline=baseline) + resp_cycles = compute_respiration_cycle_features(resp, srate, cycles, baseline=baseline, **params["features"]) - resp_cycles = clean_respiration_cycles(resp, srate, resp_cycles, baseline, **params['cycle_clean']) + if params['cycle_clean'] is not None: + resp_cycles = clean_respiration_cycles(resp, srate, resp_cycles, baseline, **params['cycle_clean']) return resp, resp_cycles @@ -106,12 +115,43 @@ def get_respiration_baseline(resp, srate, baseline_mode='manual', baseline=None) return baseline + +def detect_respiration_cycles(resp, srate, method="crossing_baseline", **method_kwargs): + """ + Detect respiration with several methods: + * "crossing_baseline": internally use detect_respiration_cycles_crossing_baseline() + * "co2" : internally use detect_respiration_cycles_co2() + + Parameters + ---------- + resp: np.array + Preprocess traces of respiratory signal. + srate: float + Sampling rate + method: 'crossing_baseline' / 'co2' + Which method for respiration. + **method_kwargs: + All other option are routed to the sub function. + Returns + ------- + cycles: np.array + Indices of inspiration and expiration. shape=(num_cycle, 3) + with [index_inspi, index_expi, index_next_inspi] + """ + if method == "crossing_baseline": + cycles = detect_respiration_cycles_crossing_baseline(resp, srate, **method_kwargs) + elif method == "co2": + cycles = detect_respiration_cycles_co2(resp, srate, **method_kwargs) + else: + raise ValueError(f"detect_respiration_cycles(): {method} do not exists") + return cycles -def detect_respiration_cycles(resp, srate, baseline_mode='manual', baseline=None, + +def detect_respiration_cycles_crossing_baseline(resp, srate, baseline_mode='manual', baseline=None, epsilon_factor1=10, epsilon_factor2=5, inspiration_adjust_on_derivative=False): """ - Detect respiration cycles based on: + Detect respiration by cycles based on: * crossing zeros (or crossing baseline) * some cleanning with euristicts @@ -157,51 +197,8 @@ def detect_respiration_cycles(resp, srate, baseline_mode='manual', baseline=None ind_insp = np.unique(ind_insp) ind_exp, = np.nonzero((resp0 < baseline) & (resp1 >= baseline)) - keep_inds = np.searchsorted(ind_exp, ind_insp, side='right') - keep_inds = keep_inds[keep_inds ind_insp[0]) & (ind_exp < ind_insp[-1]) - ind_exp = ind_exp[mask] - - # corner cases several ind_insp at the beginning/end - n = np.sum(ind_insp < ind_exp[0]) - if n > 1: - ind_insp = ind_insp[n - 1:] - n = np.sum(ind_insp > ind_exp[-1]) - if n > 1: - ind_insp = ind_insp[: - (n - 1)] + ind_insp, ind_exp = interleave_insp_exp(ind_insp, ind_exp, remove_first_insp=True, remove_first_exp=False) if inspiration_adjust_on_derivative: # lets find local minima on second derivative @@ -235,11 +232,124 @@ def detect_respiration_cycles(resp, srate, baseline_mode='manual', baseline=None return cycles +def detect_respiration_cycles_co2(co2_raw, srate, thresh_inspi_factor=0.08, thresh_expi_factor=0.05, + clean_by_mid_value=True): + """ + Detect respiration for CO2 sensor. + + Parameters + ---------- + co2_raw: np.array + Preprocess traces of respiratory signal. + srate: float + Sampling rate + thresh_inspi_factor: float, default 0.05 + Fraction of the min derivative for sertting threshold for inspiration + thresh_expi_factor: float, default 0.05 + Fraction of the min derivative for sertting threshold for expiration + clean_by_mid_value: bool, default True + Remove strange cycle by mid value. + + Returns + ------- + cycles: np.array + Indices of inspiration and expiration. shape=(num_cycle, 3) + with [index_inspi, index_expi, index_next_inspi] + """ + + + # detect resp on CO2 signal + co2_gradient = np.gradient(co2_raw) + min_, max_ = np.min(co2_gradient), np.max(co2_gradient) + + thresh_inspi = min_ * thresh_inspi_factor + thresh_expi = max_ * thresh_expi_factor + + ind_insp = np.flatnonzero((co2_gradient[:-1] >= thresh_inspi) & (co2_gradient[1:] < thresh_inspi)) + ind_exp = np.flatnonzero((co2_gradient[:-1] <= thresh_expi) & (co2_gradient[1:] > thresh_expi)) + + + ind_insp, ind_exp = interleave_insp_exp(ind_insp, ind_exp, remove_first_insp=True, remove_first_exp=False) + + # simple clean by separting up and down values + if clean_by_mid_value: + insp_values = co2_raw[ind_insp] + exp_values = co2_raw[ind_exp] + mid_value = (np.median(insp_values) + np.median(exp_values)) / 2 + remove_inds = np.flatnonzero(exp_values > mid_value) + keep = np.ones(ind_insp.size, dtype=bool) + keep[remove_inds] = False + ind_insp = ind_insp[keep] + keep = np.ones(ind_exp.size, dtype=bool) + keep[remove_inds] = False + ind_exp = ind_exp[keep] + + # import matplotlib.pyplot as plt + # fig, axs = plt.subplots(nrows=2, sharex=True) + # ax = axs[0] + # ax.plot(co2_raw) + # ax.plot(smoothed_co2) + # ax.scatter(ind_insp, smoothed_co2[ind_insp], color='g') + # ax.scatter(ind_exp, smoothed_co2[ind_exp], color='r') + # ax = axs[1] + # ax.plot(co2_gradient) + # ax.axhline(thresh_inspi, color='g') + # ax.axhline(thresh_expi, color='r') + # plt.show() + + cycles = np.zeros((ind_insp.size - 1, 3), dtype='int64') + cycles[:, 0] = ind_insp[:-1] + cycles[:, 1] = ind_exp + cycles[:, 2] = ind_insp[1:] + + return cycles + + +def _ensure_interleave(ind0, ind1, remove_first=True): + """ + Clean ind0 so they are interlevaed into ind1. + """ + keep_inds = np.searchsorted(ind1, ind0, side='right') + keep = np.ones(ind0.size, dtype=bool) + ind_bad = np.flatnonzero(np.diff(keep_inds) == 0) + if remove_first: + keep[ind_bad] = False + else: + keep[ind_bad + 1] = False + ind0_clean = ind0[keep] + return ind0_clean + + +def interleave_insp_exp(ind_insp, ind_exp, remove_first_insp=True, remove_first_exp=False): + """ + Ensure index of inspiration and expiration are interleaved. + + Ensure also that it start and stop with inspiration so that ind_insp.size == ind_exp.size + 1 + """ + + ind_exp = _ensure_interleave(ind_exp, ind_insp, remove_first=remove_first_exp) + + ind_insp = _ensure_interleave(ind_insp, ind_exp, remove_first=remove_first_insp) + + if np.any(ind_exp < ind_insp[0]): + ind_exp = ind_exp[ind_exp>ind_insp[0]] + if np.any(ind_exp > ind_insp[-1]): + ind_exp = ind_exp[ind_exp 1: + ind_insp = ind_insp[n - 1:] + n = np.sum(ind_insp > ind_exp[-1]) + if n > 1: + ind_insp = ind_insp[: - (n - 1)] + + return ind_insp, ind_exp -def compute_respiration_cycle_features(resp, srate, cycles, baseline=None): +def compute_respiration_cycle_features(resp, srate, cycles, baseline=None, compute_volume=True, compute_amplitude=True): """ Compute respiration features cycle by cycle @@ -311,26 +421,37 @@ def compute_respiration_cycle_features(resp, srate, cycles, baseline=None): df['cycle_freq'] = 1. / df['cycle_duration'] df['cycle_ratio'] = df['inspi_duration'] / df['cycle_duration'] - for k in ('inspi_volume', 'expi_volume', 'total_amplitude', 'inspi_amplitude', 'expi_amplitude'): - df[k] = pd.Series(dtype='float64') + + if compute_volume: + for k in ('inspi_volume', 'expi_volume', ): + df[k] = pd.Series(dtype='float64') + + if compute_amplitude: + for k in ('total_amplitude', 'inspi_amplitude', 'expi_amplitude'): + df[k] = pd.Series(dtype='float64') #missing cycle mask = (ix2 == -1) df.loc[mask, ['expi_time', 'cycle_duration', 'inspi_duration', 'expi_duration', 'cycle_freq']] = np.nan - for c in range(n): - i1, i2, i3 = ix1[c], ix2[c], ix3[c] - if i2 == -1: - #this is a missing cycle in the middle - continue - - df.at[c, 'inspi_volume'] = np.abs(np.sum(resp[i1:i2])) / srate - df.at[c, 'expi_volume'] = np.abs(np.sum(resp[i2:i3])) / srate - df.at[c, 'inspi_amplitude'] = np.max(np.abs(resp[i1:i2])) - df.at[c, 'expi_amplitude'] = np.max(np.abs(resp[i2:i3])) + if compute_volume or compute_amplitude: + for c in range(n): + i1, i2, i3 = ix1[c], ix2[c], ix3[c] + if i2 == -1: + #this is a missing cycle in the middle + continue + if compute_volume: + df.at[c, 'inspi_volume'] = np.abs(np.sum(resp[i1:i2])) / srate + df.at[c, 'expi_volume'] = np.abs(np.sum(resp[i2:i3])) / srate + if compute_amplitude: + df.at[c, 'inspi_amplitude'] = np.max(np.abs(resp[i1:i2])) + df.at[c, 'expi_amplitude'] = np.max(np.abs(resp[i2:i3])) - df['total_amplitude'] = df['inspi_amplitude'] + df['expi_amplitude'] - df['total_volume'] = df['inspi_volume'] + df['expi_volume'] + if compute_amplitude: + df['total_amplitude'] = df['inspi_amplitude'] + df['expi_amplitude'] + + if compute_volume: + df['total_volume'] = df['inspi_volume'] + df['expi_volume'] return resp_cycles diff --git a/tests/test_respiration.py b/tests/test_respiration.py index 344f85a..372d977 100644 --- a/tests/test_respiration.py +++ b/tests/test_respiration.py @@ -25,9 +25,9 @@ def test_detect_respiration_cycles(): resp = preprocess(raw_resp, srate, band=25., btype='lowpass', ftype='bessel', order=5, normalize=False) - cycles = detect_respiration_cycles(resp, srate, baseline_mode='median', + cycles = detect_respiration_cycles(resp, srate, method="crossing_baseline", baseline_mode='median', inspiration_adjust_on_derivative=False) - + print(cycles.shape) # clean_respiration_cycles From 2b327cfb15292d5b54ecfea7e56979f7a0710403 Mon Sep 17 00:00:00 2001 From: Samuel Garcia Date: Thu, 1 Feb 2024 12:24:11 +0100 Subject: [PATCH 2/3] Fix bug in crossing baseline method --- physio/respiration.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/physio/respiration.py b/physio/respiration.py index b2e39b4..3c373fb 100644 --- a/physio/respiration.py +++ b/physio/respiration.py @@ -67,7 +67,8 @@ def compute_respiration(raw_resp, srate, parameter_preset='human_airflow', param resp_cycles = compute_respiration_cycle_features(resp, srate, cycles, baseline=baseline, **params["features"]) - + + if params['cycle_clean'] is not None: resp_cycles = clean_respiration_cycles(resp, srate, resp_cycles, baseline, **params['cycle_clean']) @@ -191,15 +192,17 @@ def detect_respiration_cycles_crossing_baseline(resp, srate, baseline_mode='manu ind_insp, = np.nonzero((resp0 >= baseline_insp) & (resp1 < baseline_insp)) ind_insp_no_clean = ind_insp.copy() - keep_inds = np.searchsorted(ind_insp, ind_dw, side='left') + keep_inds = np.searchsorted(ind_insp, ind_dw, side='right') keep_inds = keep_inds[keep_inds > 0] ind_insp = ind_insp[keep_inds - 1] ind_insp = np.unique(ind_insp) - ind_exp, = np.nonzero((resp0 < baseline) & (resp1 >= baseline)) - + + ind_insp, ind_exp = interleave_insp_exp(ind_insp, ind_exp, remove_first_insp=True, remove_first_exp=False) + + if inspiration_adjust_on_derivative: # lets find local minima on second derivative # this can be slow From db22e69e49681509fab0df6faeb069bfa9a1c0aa Mon Sep 17 00:00:00 2001 From: Samuel Garcia Date: Tue, 6 Feb 2024 15:37:57 +0100 Subject: [PATCH 3/3] Optionally skip preprocessing in compute_respiration --- physio/respiration.py | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/physio/respiration.py b/physio/respiration.py index 3c373fb..d154f59 100644 --- a/physio/respiration.py +++ b/physio/respiration.py @@ -42,15 +42,17 @@ def compute_respiration(raw_resp, srate, parameter_preset='human_airflow', param recursive_update(params, parameters) # filter and smooth : more or less 2 times a low pass - center = np.mean(raw_resp) - resp = raw_resp - center - resp = preprocess(resp, srate, **params['preprocess']) - if params['smooth'] is not None: - resp = smooth_signal(resp, srate, **params['smooth']) - resp += center - + if params['preprocess'] is not None and params['smooth'] is not None: + center = np.mean(raw_resp) + resp = raw_resp - center + if params['preprocess'] is not None: + resp = preprocess(resp, srate, **params['preprocess']) + if params['smooth'] is not None: + resp = smooth_signal(resp, srate, **params['smooth']) + resp += center + else: + resp = raw_resp - baseline = get_respiration_baseline(resp, srate, **params['baseline']) @@ -201,6 +203,15 @@ def detect_respiration_cycles_crossing_baseline(resp, srate, baseline_mode='manu ind_insp, ind_exp = interleave_insp_exp(ind_insp, ind_exp, remove_first_insp=True, remove_first_exp=False) + # import matplotlib.pyplot as plt + # fig, ax = plt.subplots() + # ax.plot(resp) + # ax.scatter(ind_insp, resp[ind_insp], color='g') + # ax.scatter(ind_exp, resp[ind_exp], color='r') + # ax.axhline(baseline_dw) + # ax.axhline(baseline_insp) + # ax.axhline(baseline) + # plt.show() if inspiration_adjust_on_derivative: