From 3e0d83ca9b7ab9277b26c961ea3e9a4c9ee059c0 Mon Sep 17 00:00:00 2001 From: Tyler Date: Thu, 30 Jul 2026 16:15:33 -0700 Subject: [PATCH] implement the new preferences in the backend --- flares.py | 409 +++++++++++++++++++++++---------- main.py | 44 ++-- src/shared/flaresbasewidget.py | 63 +++-- 3 files changed, 351 insertions(+), 165 deletions(-) diff --git a/flares.py b/flares.py index 17a2421..405a3e7 100644 --- a/flares.py +++ b/flares.py @@ -135,27 +135,40 @@ OPTODE_PLACEMENT: bool SHOW_OPTODE_NAMES: bool SHORT_CHANNELS: bool +LONG_CHANNELS: bool SHORT_CHANNELS_THRESHOLD: float LONG_CHANNELS_THRESHOLD: float HEART_RATE: bool SECONDS_TO_STRIP_HR: int +HR_LOW_FREQ: int +HR_HIGH_FREQ: int +HR_SEARCH_MIN: int +HR_SEARCH_MAX: int MAX_LOW_HR: int MAX_HIGH_HR: int SMOOTHING_WINDOW_HR: int HEART_RATE_WINDOW: int SCI: bool +SCI_USE_HEART_RATE_BAND: bool +SCI_LOW_FREQ: float +SCI_HIGH_FREQ: float SCI_TIME_WINDOW: int SCI_THRESHOLD: float SNR: bool -# SNR_TIME_WINDOW : int #TODO: is this needed? SNR_THRESHOLD: float +SNR_SIGNAL_LOW_FREQ: float +SNR_SIGNAL_HIGH_FREQ: float +SNR_NOISE_LOW_FREQ: float +SNR_NOISE_HIGH_FREQ: float PSP: bool PSP_TIME_WINDOW: int PSP_THRESHOLD: float +PSP_LOW_FREQ: float +PSP_HIGH_FREQ: float COEFF_VAR: bool COEFF_VAR_THRESHOLD: int @@ -166,6 +179,8 @@ MAD_THRESHOLD: int PSD_NOISE: bool TARGET_FREQ_DIV: int DB_LIMIT: int +PSD_MIN_FREQ: float +PSD_TARGET_BANDWIDTH: float SENSOR_DROPOUT: bool SENSOR_DROPOUT_VARIANCE_THRESHOLD: float @@ -189,27 +204,48 @@ PPF_UPPER_WAVELENGTH: float ENHANCE_NEGATIVE_CORRELATION: bool FILTER: bool +FILTER_ALGORITHM: list L_FREQ: float H_FREQ: float L_TRANS_BANDWIDTH: float H_TRANS_BANDWIDTH: float +IIR_TYPE: list +IIR_ORDER: int +FILTER_LENGTH: str +FILTER_PHASE: list +FIR_WINDOW: list +FIR_DESIGN: list +IIR_OUTPUT: list +PASSBAND_RIPPLE: float +STOPBAND_ATTENUATION: float +FILTER_PAD: list +SKIP_BY_ANNOTATION: list +FILTER_N_JOBS: int +EVENTS: bool +EVENT_ID: str +EVENT_REGEX: str +EVENT_CHUNK_DURATION: float + +EPOCHS: bool EPOCH_HANDLING: str MAX_SHIFT: int T_MIN: int T_MAX: int +BASELINE: list +REJECT_EPOCHS: bool +REJECT_HBO_THRESHOLD: float RESAMPLE: bool RESAMPLE_FREQ: int -STIM_DUR: float HRF_MODEL: str +STIM_DUR: float +FIR_DELAYS: range DRIFT_MODEL: str HIGH_PASS: float DRIFT_ORDER: int -FIR_DELAYS: range MIN_ONSET: int OVERSAMPLING: int -REMOVE_EVENTS: list SHORT_CHANNEL_REGRESSION: bool NOISE_MODEL: str @@ -828,7 +864,7 @@ def scalp_coupling_index_windowed_raw(data, time_window: float = 3.0, l_freq: fl return data, scores, times -def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5): +def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5, time_window: int = 3, threshold: float = 0.6): """ Calculate the scalp coupling index (SCI) and identify bad channels based on a threshold. @@ -852,21 +888,21 @@ def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5): print("Calculating scalp coupling index...") # Compute the SCI - _, scores, times = scalp_coupling_index_windowed_raw(data, time_window=SCI_TIME_WINDOW, l_freq=l_freq, h_freq=h_freq) + _, scores, times = scalp_coupling_index_windowed_raw(data, time_window=time_window, l_freq=l_freq, h_freq=h_freq) # Identify channels that don't meet the provided threshold print("Identifying channels that do not meet the threshold...") sci = scores.mean(axis=1) - data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), sci < SCI_THRESHOLD)) + data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), sci < threshold)) # Determine the colors based on the threshold, and create the figures print("Creating the figures...") - color_stops = ([0.0, SCI_THRESHOLD, SCI_THRESHOLD+0.1, 0.8, 1.0], [0.0, SCI_THRESHOLD, SCI_THRESHOLD, 1.0]) - fig1, fig2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, SCI_THRESHOLD, "Scalp Coupling Index") + color_stops = ([0.0, threshold, threshold+0.1, 0.8, 1.0], [0.0, threshold, threshold, 1.0]) + fig1, fig2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, threshold, "Scalp Coupling Index") print("Successfully calculated scalp coupling index.") - return list(compress(cast(list[str], getattr(data, "ch_names")), sci < SCI_THRESHOLD)), fig1, fig2 + return list(compress(cast(list[str], getattr(data, "ch_names")), sci < threshold)), fig1, fig2 @@ -933,7 +969,7 @@ def get_hbo_hbr_picks(raw): return hbo_picks, hbr_picks, hbo_wl, hbr_wl -def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2): +def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2, short_channels_threshold=0.015): """ Interpolate bad fNIRS channels using a distance-weighted average of nearby good channels. @@ -984,7 +1020,7 @@ def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2) hbo_ch = raw.ch_names[hbo_picks[i]] hbr_ch = raw.ch_names[hbr_picks[i]] is_bad = (hbo_ch in raw.info['bads']) or (hbr_ch in raw.info['bads']) - is_short = pair_distances[i] < SHORT_CHANNELS_THRESHOLD + is_short = pair_distances[i] < short_channels_threshold if is_bad: bad_pairs.append(i) @@ -1168,7 +1204,7 @@ def calculate_signal_noise_ratio(data): -def calculate_peak_power(data: BaseRaw, l_freq: float = 0.7, h_freq: float = 1.5) -> tuple[list[str], Figure, Figure]: +def calculate_peak_power(data: BaseRaw, time_window: int = 3, threshold: float = 0.1, l_freq: float = 0.7, h_freq: float = 1.5) -> tuple[list[str], Figure, Figure]: """ Calculate peak spectral power (PSP) for fNIRS channels and identify bad channels. @@ -1191,18 +1227,18 @@ def calculate_peak_power(data: BaseRaw, l_freq: float = 0.7, h_freq: float = 1.5 # Compute the PSP - _, scores, times = cast(tuple[NDArray[float64], NDArray[float64], list[tuple[float]]], peak_power(data, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=l_freq, h_freq=h_freq)) + _, scores, times = cast(tuple[NDArray[float64], NDArray[float64], list[tuple[float]]], peak_power(data, time_window=time_window, threshold=threshold, l_freq=l_freq, h_freq=h_freq)) # Identify channels that don't meet the provided threshold psp = scores.mean(axis=1) - data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), psp < PSP_THRESHOLD)) + data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), psp < threshold)) # Determine the colors based on the threshold, and create the figures - color_stops = ([0.0, PSP_THRESHOLD, PSP_THRESHOLD+0.1, PSP_THRESHOLD+0.2, 1.0], [0.0, PSP_THRESHOLD, PSP_THRESHOLD, 1.0]) - psp1, psp2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, PSP_THRESHOLD, "Peak Spectral Power") + color_stops = ([0.0, threshold, threshold+0.1, threshold+0.2, 1.0], [0.0, threshold, threshold, 1.0]) + psp1, psp2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, threshold, "Peak Spectral Power") - return list(compress(cast(list[str], getattr(data, "ch_names")), psp < PSP_THRESHOLD)), psp1, psp2 + return list(compress(cast(list[str], getattr(data, "ch_names")), psp < threshold)), psp1, psp2 def mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_range, bad_noise, bad_disp): @@ -1279,18 +1315,49 @@ def mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_range, bad_noi -def filter_the_data(raw_haemo): +def filter_the_data( + raw_haemo, + filter_algorithm, + l_freq, + h_freq, + l_trans_bandwidth, + h_trans_bandwidth, + iir_type, + iir_order, + filter_length, + filter_phase, + fir_window, + fir_design, + iir_output, + passband_ripple, + stopband_attenuation, + filter_pad, + skip_by_annotation, + filter_n_jobs, + verbosity + ): # --- STEP 5: Filtering (0.01-0.2 Hz bandpass) --- fig_filter = raw_haemo.compute_psd(fmax=3).plot( average=True, color="r", show=False, amplitude=True ) - if L_FREQ == 0 and H_FREQ != 0: - raw_haemo = raw_haemo.filter(l_freq=None, h_freq=H_FREQ, h_trans_bandwidth=H_TRANS_BANDWIDTH) - elif L_FREQ != 0 and H_FREQ == 0: - raw_haemo = raw_haemo.filter(l_freq=L_FREQ, h_freq=None, l_trans_bandwidth=L_TRANS_BANDWIDTH) - elif L_FREQ != 0 and H_FREQ != 0: - raw_haemo = raw_haemo.filter(l_freq=L_FREQ, h_freq=H_FREQ, l_trans_bandwidth=L_TRANS_BANDWIDTH, h_trans_bandwidth=H_TRANS_BANDWIDTH) + if l_freq == 0 and h_freq != 0: + if filter_algorithm == "fir": + raw_haemo = raw_haemo.filter(l_freq=None, h_freq=h_freq, filter_length=filter_length, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, phase=filter_phase, fir_window=fir_window, fir_design=fir_design, pad=filter_pad, verbose=verbosity) + else: + raw_haemo = raw_haemo.filter(l_freq=None, h_freq=h_freq, filter_length=filter_length, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, pad=filter_pad, verbose=verbosity) + + elif l_freq != 0 and h_freq == 0: + if filter_algorithm == "fir": + raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=None, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, phase=filter_phase, fir_window=fir_window, fir_design=fir_design, pad=filter_pad, verbose=verbosity) + else: + raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=None, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, pad=filter_pad, verbose=verbosity) + + elif l_freq != 0 and h_freq != 0: + if filter_algorithm == "fir": + raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=h_freq, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, phase=filter_phase, fir_window=fir_window, fir_design=fir_design, pad=filter_pad, verbose=verbosity) + else: + raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=h_freq, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, pad=filter_pad, verbose=verbosity) else: print("No filter") #raw_haemo = raw_haemo.filter(l_freq=None, h_freq=0.4, h_trans_bandwidth=0.2) @@ -1307,23 +1374,30 @@ def filter_the_data(raw_haemo): -def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline): +def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline, max_shift, reject_epochs, reject_hbo_threshold): """ Attempts to create epochs, shifting event times slightly if sample collisions are detected. """ shift_increment = 1.0 / raw.info['sfreq'] # The duration of exactly one sample - #TODO: User expose this - reject_criteria = dict(hbo=80e-7) - for attempt in range(MAX_SHIFT): # Limit attempts to avoid infinite loops + + for attempt in range(max_shift): # Limit attempts to avoid infinite loops try: - epochs = Epochs( - raw, events, event_id=event_dict, - tmin=tmin, tmax=tmax, baseline=baseline, - reject=reject_criteria, - preload=True, verbose=False - ) + if reject_epochs: + epochs = Epochs( + raw, events, event_id=event_dict, + tmin=tmin, tmax=tmax, baseline=baseline, + reject=reject_hbo_threshold, + preload=True, verbose=False + ) + else: + epochs = Epochs( + raw, events, event_id=event_dict, + tmin=tmin, tmax=tmax, baseline=baseline, + preload=True, verbose=False + ) return epochs + except RuntimeError as e: if "Event time samples were not unique" in str(e): # Find duplicates in the events array (column 0 is the sample index) @@ -1344,13 +1418,16 @@ def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline): -def epochs_calculations(raw_haemo, events, event_dict): +def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift, t_min, t_max, baseline, reject_epochs, reject_hbo_threshold): fig_epochs = [] # List to store figures - if EPOCH_HANDLING == 'shift': - epochs = safe_create_epochs(raw=raw_haemo, events=events, event_dict=event_dict, tmin=T_MIN, tmax=T_MAX, baseline=(None, 0)) + if epoch_handling == 'shift': + epochs = safe_create_epochs(raw=raw_haemo, events=events, event_dict=event_dict, tmin=t_min, tmax=t_max, baseline=baseline, max_shift=max_shift, reject_epochs=reject_epochs, reject_hbo_threshold=reject_hbo_threshold) else: - epochs = Epochs(raw_haemo, events, event_id=event_dict, tmin=T_MIN, tmax=T_MAX, baseline=(None, 0)) + if reject_epochs: + epochs = Epochs(raw_haemo, events, event_id=event_dict, tmin=t_min, tmax=t_max, baseline=baseline, reject=reject_hbo_threshold) + else: + epochs = Epochs(raw_haemo, events, event_id=event_dict, tmin=t_min, tmax=t_max, baseline=baseline) # Make a copy of the epochs and drop bad ones epochs2 = epochs.copy() @@ -1502,7 +1579,25 @@ def epochs_calculations(raw_haemo, events, event_dict): -def make_design_matrix(raw_haemo): +def make_design_matrix( + raw_haemo, + resample, + resample_freq, + stim_dur, + hrf_model, + drift_model, + high_pass, + drift_order, + fir_delays, + min_onset, + oversampling, + short_channel_regression, + short_channels, + long_channels, + short_channels_threshold, + long_channels_threshold, + folding_bypass + ): # events_to_remove = REMOVE_EVENTS events_to_remove = "" @@ -1515,37 +1610,39 @@ def make_design_matrix(raw_haemo): description=[ann['description'] for ann in filtered_annotations] ) - if SHORT_CHANNELS: - short_chans = get_short_channels(raw_haemo, max_dist=SHORT_CHANNELS_THRESHOLD) - raw_haemo = get_long_channels(raw_haemo, min_dist=SHORT_CHANNELS_THRESHOLD, max_dist=LONG_CHANNELS_THRESHOLD) + if short_channels: + short_chans = get_short_channels(raw_haemo, max_dist=short_channels_threshold) + if long_channels: + raw_haemo = get_long_channels(raw_haemo, min_dist=long_channels_threshold, max_dist=long_channels_threshold) + else: short_chans = None # Set the new annotations raw_haemo.set_annotations(new_annot) - if RESAMPLE: - raw_haemo.resample(RESAMPLE_FREQ, npad="auto") + if resample: + raw_haemo.resample(resample_freq, npad="auto") raw_haemo._data = raw_haemo._data * 1e6 try: - short_chans.resample(RESAMPLE_FREQ) + short_chans.resample(resample_freq) except: pass design_matrix = make_first_level_design_matrix( raw=raw_haemo, - stim_dur=STIM_DUR, - hrf_model=HRF_MODEL, - drift_model=DRIFT_MODEL, - high_pass=HIGH_PASS, - drift_order=DRIFT_ORDER, - fir_delays=FIR_DELAYS, - min_onset=MIN_ONSET, - oversampling=OVERSAMPLING + stim_dur=stim_dur, + hrf_model=hrf_model, + drift_model=drift_model, + high_pass=high_pass, + drift_order=drift_order, + fir_delays=fir_delays, + min_onset=min_onset, + oversampling=oversampling ) # 3) Average and Append Short Channels - if SHORT_CHANNEL_REGRESSION and not FOLDING_BYP: + if short_channel_regression and not folding_bypass: if short_chans is not None and len(short_chans.ch_names) > 0: ch_types = short_chans.get_channel_types() @@ -2724,7 +2821,7 @@ def plot_fir_model_results( -def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]: +def load_snirf(file_path: str, downsample_frequency: int, verbosity: bool) -> tuple[BaseRaw, Figure]: """ Loads a snirf file, optionally drops channels, downsamples, and creates a figure showing the results. @@ -2745,7 +2842,7 @@ def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]: """ # Read the snirf file - raw = read_raw_snirf(file_path, preload=True, verbose=VERBOSITY) # type: ignore + raw = read_raw_snirf(file_path, preload=True, verbose=verbosity) # type: ignore #raw.load_data(verbose=VERBOSITY) # type: ignore redundant since preload is set to true # TODO: Why was this commented again? @@ -2772,7 +2869,7 @@ def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]: if DOWNSAMPLE: logger.info("Downsample was specified.") sfreq_old = getattr(raw, "info")["sfreq"] - raw.resample(DOWNSAMPLE_FREQUENCY, verbose=VERBOSITY) # type: ignore + raw.resample(downsample_frequency, verbose=verbosity) # type: ignore sfreq_new = getattr(raw, "info")["sfreq"] logger.info(f"Finished downsampling. Old frequency: {sfreq_old}. New frequency: {sfreq_new}.") @@ -4437,7 +4534,7 @@ def iqr_threshold(coeffs: NDArray[float64], k: float = 1.5) -> floating[Any]: -def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: int = 3) -> NDArray[float64]: +def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: int = 3, iqr: float = 1.5) -> NDArray[float64]: """ Denoises a signal using wavelet decomposition and IQR-based thresholding on detail coefficients. @@ -4463,7 +4560,7 @@ def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: i # Threshold detail coefficients to reduce noise for cD in coeffs[1:]: - threshold = iqr_threshold(cD, IQR) + threshold = iqr_threshold(cD, iqr) cD_thresh = np.sign(cD) * np.maximum(np.abs(cD) - threshold, 0.0) # np.where((cD < lower) | (cD > upper), 0, cD) cD_thresh = cD_thresh.astype(float64) denoised_coeffs.append(cD_thresh) @@ -4474,7 +4571,7 @@ def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: i -def calculate_and_apply_wavelet(data: BaseRaw) -> tuple[BaseRaw, Figure]: +def calculate_and_apply_wavelet(data: BaseRaw, wavelet_type: str, wavelet_level: int, iqr: float, verbosity: bool) -> tuple[BaseRaw, Figure]: """ Applies a wavelet IQR denoising filter to the data and generates a plot. @@ -4496,17 +4593,17 @@ def calculate_and_apply_wavelet(data: BaseRaw) -> tuple[BaseRaw, Figure]: # Denoise the data logger.info("Denoising the data...") - loaded_data: NDArray[float64] = data.get_data(verbose=VERBOSITY) # type: ignore + loaded_data: NDArray[float64] = data.get_data(verbose=verbosity) # type: ignore denoised_data = np.zeros_like(loaded_data) logger.info("Calculating the IQR, decomposing the signal, and thresholding the coefficients...") for ch in range(loaded_data.shape[0]): - denoised_data[ch, :] = wavelet_iqr_denoise(loaded_data[ch, :], wavelet=WAVELET_TYPE, level=WAVELET_LEVEL) + denoised_data[ch, :] = wavelet_iqr_denoise(loaded_data[ch, :], wavelet=wavelet_type, level=wavelet_level, iqr=iqr) # Reconstruct the data with the annotations logger.info("Reconstructing the data with annotations...") - raw_with_tddr_and_wavelet = RawArray(denoised_data, cast(Info, data.info), verbose=VERBOSITY) - raw_with_tddr_and_wavelet.set_annotations(data.annotations.copy(), verbose=VERBOSITY) # type: ignore + raw_with_tddr_and_wavelet = RawArray(denoised_data, cast(Info, data.info), verbose=verbosity) + raw_with_tddr_and_wavelet.set_annotations(data.annotations.copy(), verbose=verbosity) # type: ignore # Create a figure for the results logger.info("Creating the figure...") @@ -4521,7 +4618,7 @@ def calculate_and_apply_wavelet(data: BaseRaw) -> tuple[BaseRaw, Figure]: -def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None) -> tuple[float, NDArray[float64], NDArray[float64]]: +def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None, seconds_to_strip_hr: int, verbosity: bool) -> tuple[float, NDArray[float64], NDArray[float64]]: """ Extract and trim short-channel fNIRS signal for heart rate analysis. @@ -4546,16 +4643,16 @@ def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None) # If a short channel exists, use it for our signal. Otherwise just take the first channel in the data # TODO: Find a better way around this if short_chans is not None: - signal = cast(NDArray[float64], short_chans.get_data(picks=[0], verbose=VERBOSITY))[0] # type: ignore + signal = cast(NDArray[float64], short_chans.get_data(picks=[0], verbose=verbosity))[0] # type: ignore else: - signal = cast(NDArray[float64], data.get_data(picks=[0], verbose=VERBOSITY))[0] # type: ignore + signal = cast(NDArray[float64], data.get_data(picks=[0], verbose=verbosity))[0] # type: ignore # Calculate the sampling frequency sfreq = cast(int, data.info['sfreq']) # Trim start and end of the signal to remove edge artifacts - logger.info(f"Removing {SECONDS_TO_STRIP_HR} seconds from the beginning and end of the file...") - strip_samples = int(sfreq * SECONDS_TO_STRIP_HR) + logger.info(f"Removing {seconds_to_strip_hr} seconds from the beginning and end of the file...") + strip_samples = int(sfreq * seconds_to_strip_hr) signal_trimmed = signal[strip_samples:-strip_samples] times_trimmed = cast(NDArray[float64], getattr(data, "times"))[strip_samples:-strip_samples] @@ -4563,7 +4660,7 @@ def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None) -def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64]) -> tuple[NDArray[float64], float]: +def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64], hr_low_freq, hr_high_freq, max_low_hr, max_high_hr, smoothing_window_hr) -> tuple[NDArray[float64], float]: """ Calculate and smooth heart rate from a trimmed signal using NeuroKit. @@ -4585,19 +4682,19 @@ def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64] # Filter signal to isolate heart rate frequencies and detect peaks logger.info("Filtering the signal and detecting peaks...") - signal_filtered = cast(NDArray[float64], nk.signal_filter(signal_trimmed, sampling_rate=sfreq, lowcut=0.8, highcut=2.5)) # type: ignore + signal_filtered = cast(NDArray[float64], nk.signal_filter(signal_trimmed, sampling_rate=sfreq, lowcut=hr_low_freq, highcut=hr_high_freq)) # type: ignore peaks_dict = cast(dict[str, Any], nk.signal_findpeaks(signal_filtered)) # type: ignore peaks = peaks_dict['Peaks'] hr = cast(NDArray[float64], nk.signal_rate(peaks, sampling_rate=sfreq, desired_length=len(signal_trimmed))) # type: ignore - hr_clean = np.clip(hr, MAX_LOW_HR, MAX_HIGH_HR) + hr_clean = np.clip(hr, max_low_hr, max_high_hr) # Smooth heart rate time series by replacing spikes with local rolling mean and calculate the mean logger.info("Smoothing the signal and calculating the mean...") hr_series = pd.Series(hr_clean) - local_median = hr_series.rolling(window=SMOOTHING_WINDOW_HR, center=True, min_periods=1).median() + local_median = hr_series.rolling(window=smoothing_window_hr, center=True, min_periods=1).median() spikes = hr_series > (local_median + 10) smoothed_values = hr_series.copy() - smoothed_spikes = hr_series.rolling(window=SMOOTHING_WINDOW_HR, center=True, min_periods=1).mean() + smoothed_spikes = hr_series.rolling(window=smoothing_window_hr, center=True, min_periods=1).mean() smoothed_values[spikes] = smoothed_spikes[spikes] hr_smooth_nk = cast(NDArray[float64], smoothed_values.to_numpy()) # type: ignore mean_hr_nk = hr_smooth_nk.mean() @@ -4612,7 +4709,7 @@ def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64] -def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64]) -> tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float]: +def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64], search_min, search_max) -> tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float]: """ Estimate heart rate using spectral analysis on a high-pass filtered signal. @@ -4628,7 +4725,7 @@ def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64]) - tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float] - NDArray[floating[Any]]: Frequencies converted to beats per minute (BPM). - NDArray[float64]: Power spectral density (PSD) of the signal. - - np.ndarray[Any, np.dtype[np.bool_]]: Boolean mask indicating frequencies within heart rate range (30-300 BPM). + - np.ndarray[Any, np.dtype[np.bool_]]: Boolean mask indicating frequencies within heart rate range. - float: Estimated mean heart rate in BPM corresponding to the PSD peak within the range. """ @@ -4644,10 +4741,10 @@ def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64]) - nperseg = min(len(signal_hp), 4096) frequencies_scipy, psd_scipy = cast(tuple[NDArray[float64], NDArray[float64]], welch(signal_hp, fs=sfreq, nperseg=nperseg, noverlap=nperseg//2)) - # Convert frequency values to beats per minute (BPM) and set a heart rate range (30-300 BPM) + # Convert frequency values to beats per minute (BPM) and set a heart rate range logger.info("Converting to BPM...") freq_bpm_scipy = frequencies_scipy * 60 - freq_range_scipy = (freq_bpm_scipy > 30) & (freq_bpm_scipy < 300) + freq_range_scipy = (freq_bpm_scipy > search_min) & (freq_bpm_scipy < search_max) # Identify the peak frequency within the heart rate range and estimate the mean heart rate in BPM logger.info("Finding the mean...") @@ -4667,7 +4764,8 @@ def plot_heart_rate( hr_smooth_nk: NDArray[floating[Any]], mean_hr_nk: float, times_trimmed: NDArray[floating[Any]], - overruled: bool + overruled: bool, + hr_window: int ) -> tuple[Figure, Figure]: """ Generate plots comparing heart rate estimates from SciPy PSD and NeuroKit2. @@ -4704,7 +4802,7 @@ def plot_heart_rate( ax1.set_xlim(30, 300) ax1.plot(freq_bpm_scipy[freq_range_scipy], psd_scipy[freq_range_scipy]) # type: ignore ax1.axvline(x=mean_hr_scipy, color='red', linestyle='--', label=f'Mean HR: {mean_hr_scipy:.1f} BPM') # type: ignore - ax1.axvspan(min(mean_hr_nk - HEART_RATE_WINDOW, mean_hr_scipy - HEART_RATE_WINDOW), max(mean_hr_nk + HEART_RATE_WINDOW, mean_hr_scipy + HEART_RATE_WINDOW), color='yellow', alpha=0.3, label=f'HR Range ±{HEART_RATE_WINDOW} BPM') # type: ignore + ax1.axvspan(min(mean_hr_nk - hr_window, mean_hr_scipy - hr_window), max(mean_hr_nk + hr_window, mean_hr_scipy + hr_window), color='yellow', alpha=0.3, label=f'HR Range ±{hr_window} BPM') # type: ignore ax1.set_xlabel('Heart Rate (BPM)') # type: ignore ax1.set_ylabel('Power Spectral Density') # type: ignore ax1.set_title('PSD of fNIRS signal - Peak indicates Heart Rate') # type: ignore @@ -4962,7 +5060,7 @@ def detect_sensor_dropout(raw, threshold_ratio=0.05): -def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4): +def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4, min_freq=0.1, target_bandwith=0.2): """ Identifies channels with excessive power at high frequencies (sfreq/4), usually indicating electronic interference. @@ -4972,11 +5070,11 @@ def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4): target_freq = sfreq / freq_div # Compute PSD - spectrum = raw.compute_psd(fmin=0.1, fmax=sfreq/2) + spectrum = raw.compute_psd(fmin=min_freq, fmax=sfreq/2) psd_data, freqs = spectrum.get_data(return_freqs=True) # Find power near the target frequency - f_idx = np.where((freqs >= target_freq - 0.2) & (freqs <= target_freq + 0.2))[0] + f_idx = np.where((freqs >= target_freq - target_bandwith) & (freqs <= target_freq + target_bandwith))[0] power_at_target = np.mean(psd_data[:, f_idx], axis=1) abs_threshold = 10 ** (db_limit / 10) @@ -5003,7 +5101,7 @@ def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4): -def find_bad_channels_by_amplitude_range(raw, threshold=4.0): +def find_bad_channels_by_amplitude_range(raw, threshold=4): """Median absolute deviation""" picks = [ch for ch in raw.ch_names] data = raw.get_data(picks=picks) @@ -5094,14 +5192,14 @@ def find_bad_channels_coeff_var(raw, coeff_var_threshold=25.0): -def hr_calc(raw): - if SHORT_CHANNELS: - short_chans = get_short_channels(raw, max_dist=SHORT_CHANNELS_THRESHOLD) +def hr_calc(raw, seconds_to_strip_hr, l_freq, h_freq, search_min, search_max, max_low_hr, max_high_hr, smoothing_window_hr, hr_window, short_channels, short_channels_threshold, verbosity): + if short_channels: + short_chans = get_short_channels(raw, max_dist=short_channels_threshold) else: short_chans = None - sfreq, signal_trimmed, times_trimmed = short_channel_processing_for_hr(raw, short_chans) - hr_smooth_nk, mean_hr_nk = calculate_heart_rate_neurokit(sfreq, signal_trimmed) - freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy = calculate_heart_rate_scipy(sfreq, signal_trimmed) + sfreq, signal_trimmed, times_trimmed = short_channel_processing_for_hr(raw, short_chans, seconds_to_strip_hr=seconds_to_strip_hr, verbosity=verbosity) + hr_smooth_nk, mean_hr_nk = calculate_heart_rate_neurokit(sfreq, signal_trimmed, hr_low_freq=l_freq, hr_high_freq=h_freq, max_low_hr=max_low_hr, max_high_hr=max_high_hr, smoothing_window_hr=smoothing_window_hr) + freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy = calculate_heart_rate_scipy(sfreq, signal_trimmed, search_min=search_min, search_max=search_max) # HACK: This sucks but looking at the graphs I trust neurokit2 more overruled = False @@ -5112,7 +5210,7 @@ def hr_calc(raw): mean_hr_scipy = mean_hr_nk overruled = True - hr1, hr2 = plot_heart_rate(freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy, hr_smooth_nk, mean_hr_nk, times_trimmed, overruled) + hr1, hr2 = plot_heart_rate(freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy, hr_smooth_nk, mean_hr_nk, times_trimmed, overruled, hr_window=hr_window) fig = raw.plot_psd(show=False) raw_filtered = raw.copy().filter(0.5, 3, fir_design='firwin') @@ -5123,7 +5221,7 @@ def hr_calc(raw): # --- Parameters for PSD --- desired_bin_hz = 0.1 nperseg = int(sfreq / desired_bin_hz) - hr_range = (30, 180) # TODO: Should this not use the user defined values? + hr_range = (search_min, search_max) # --- Function to find strongest local peak --- def find_hr_from_psd(ch_data): @@ -5153,11 +5251,11 @@ def hr_calc(raw): return fig, hr1, hr2, low, high -def trim_participant_data(raw): +def trim_participant_data(raw, seconds_to_keep: float): if hasattr(raw, 'annotations') and len(raw.annotations) > 0: # Get time of first event first_event_time = raw.annotations.onset[0] - trim_time = max(0, first_event_time - SECONDS_TO_KEEP) # Ensure we don't go negative + trim_time = max(0, first_event_time - seconds_to_keep) # Ensure we don't go negative raw.crop(tmin=trim_time) # Shift annotation onsets to match new t=0 @@ -5180,14 +5278,14 @@ def trim_participant_data(raw): return raw, fig_trimmed -def remove_bad_channels(raw, bad_channels): +def remove_bad_channels(raw, bad_channels, max_bad_channels: int): num_bad = len(bad_channels) # Check against the threshold - if num_bad > MAX_BAD_CHANNELS: + if num_bad > max_bad_channels: raise Exception( f"Data Quality Error: {num_bad} channels flagged for removal, " - f"which exceeds the limit of {MAX_BAD_CHANNELS}. To avoid this, " + f"which exceeds the limit of {max_bad_channels}. To avoid this, " f"either lower your filtering parameters or increase MAX_BAD_CHANNELS." ) @@ -5196,9 +5294,9 @@ def remove_bad_channels(raw, bad_channels): return raw -def make_and_run_glm(raw_haemo, df_design_matrix): +def make_and_run_glm(raw_haemo, df_design_matrix, noise_model, bins, n_jobs, verbosity): - glm_est = run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbose=VERBOSITY) + glm_est = run_glm(raw_haemo, df_design_matrix, noise_model=noise_model, bins=bins, n_jobs=n_jobs, verbose=verbosity) fir_cols = [col for col in df_design_matrix.columns if "_delay_" in col] if fir_cols: @@ -5265,11 +5363,11 @@ def generate_channel_results(glm_est, file_path): return df_cha -def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path): +def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path, json_location): rois_formatted = {} try: - with open(JSON_LOCATION, 'r') as f: + with open(json_location, 'r') as f: roi_data = json.load(f) for region in roi_data.get("regions_of_interest", []): @@ -5305,7 +5403,7 @@ def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path): # --------------------------------------------------------------------- if not rois_formatted: logger.error( - f"'{JSON_LOCATION}' produced zero valid ROIs — attempting automatic " + f"'{json_location}' produced zero valid ROIs — attempting automatic " f"geometric fallback (Left/Right, then Front/Back)." ) rois_formatted = _build_geometric_fallback_rois(raw_haemo) @@ -5449,6 +5547,7 @@ def haemoglobin_concentration(raw_od, file_path, override_ppf=False, ppf_lower_w def process_participant(file_path, progress_callback=None): + """Goal is to have all of the constants be located here, and nowhere else in the code to remove any ambiguity.""" # Step 0: Setting up fig_individual: dict[str, Figure] = {} @@ -5459,7 +5558,7 @@ def process_participant(file_path, progress_callback=None): } # Step 1: Preprocessing - raw = load_snirf(file_path) + raw = load_snirf(file_path=file_path, downsample_frequency=DOWNSAMPLE_FREQUENCY, verbosity=VERBOSITY) fig_raw = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Loaded Raw", show=False) fig_individual["Loaded Raw Data"] = fig_raw if progress_callback: progress_callback(1) @@ -5467,7 +5566,7 @@ def process_participant(file_path, progress_callback=None): # Step 2: Trimming if TRIM and not FOLDING_BYP: - raw, fig_trimmed = trim_participant_data(raw) + raw, fig_trimmed = trim_participant_data(raw, seconds_to_keep=SECONDS_TO_KEEP) fig_individual["Trimmed Raw Data"] = fig_trimmed if progress_callback: progress_callback(2) logger.info("Step 2 Completed.") @@ -5485,13 +5584,28 @@ def process_participant(file_path, progress_callback=None): _short_chans = get_short_channels(raw, max_dist=SHORT_CHANNELS_THRESHOLD) # JUST FOR PLOTTING THEM SEPERATELY fig_short_chans = _short_chans.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Short Channels Only", show=False, verbose=VERBOSITY) fig_individual["Short Channels Raw Data"] = fig_short_chans - raw = get_long_channels(raw, min_dist=0, max_dist=LONG_CHANNELS_THRESHOLD) + if LONG_CHANNELS: + raw = get_long_channels(raw, min_dist=0, max_dist=LONG_CHANNELS_THRESHOLD) if progress_callback: progress_callback(4) logger.info("Step 4 Completed.") # Step 5: Heart Rate if HEART_RATE and not FOLDING_BYP: - fig, hr1, hr2, low, high = hr_calc(raw) + fig, hr1, hr2, low, high = hr_calc( + raw, + seconds_to_strip_hr=SECONDS_TO_STRIP_HR, + l_freq=HR_LOW_FREQ, + h_freq=HR_HIGH_FREQ, + search_min=HR_SEARCH_MIN, + search_max=HR_SEARCH_MAX, + max_low_hr=MAX_LOW_HR, + max_high_hr=MAX_HIGH_HR, + smoothing_window_hr=SMOOTHING_WINDOW_HR, + hr_window=HEART_RATE_WINDOW, + short_channels=SHORT_CHANNELS, + short_channels_threshold=SHORT_CHANNELS_THRESHOLD, + verbosity=VERBOSITY + ) fig_individual["Power Spectral Density"] = fig fig_individual['Heart Rate - PSD'] = hr1 fig_individual['Heart Rate - Time'] = hr2 @@ -5501,10 +5615,10 @@ def process_participant(file_path, progress_callback=None): # Step 6: Scalp Coupling Index bad_sci = [] if SCI and not FOLDING_BYP: - if HEART_RATE: - bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, low, high) + if HEART_RATE and SCI_USE_HEART_RATE_BAND: + bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, l_freq=low, h_freq=high, time_window=SCI_TIME_WINDOW, threshold=SCI_THRESHOLD) else: - bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw) + bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, l_freq=SCI_LOW_FREQ, h_freq=SCI_HIGH_FREQ, time_window=SCI_TIME_WINDOW, threshold=SCI_THRESHOLD) fig_individual["Scalp Coupling Index Heatmap"] = fig_sci_1 fig_individual["Scalp Coupling Index Binary Heatmap"] = fig_sci_2 if progress_callback: progress_callback(6) @@ -5521,7 +5635,7 @@ def process_participant(file_path, progress_callback=None): # Step 8: Peak Spectral Power bad_psp = [] if PSP and not FOLDING_BYP: - bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw) + bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=PSP_LOW_FREQ, h_freq=PSP_HIGH_FREQ) fig_individual["Peak Spectral Power Heatmap"] = fig_psp1 fig_individual["Peak Spectral Power Binary Heatmap"] = fig_psp2 if progress_callback: progress_callback(8) @@ -5546,7 +5660,7 @@ def process_participant(file_path, progress_callback=None): # Step 11: Power Spectral Density Noise bad_noise = [] if PSD_NOISE and not FOLDING_BYP: - bad_noise, fig_noise = detect_spectral_noise_spike(raw, db_limit=DB_LIMIT, freq_div=TARGET_FREQ_DIV) + bad_noise, fig_noise = detect_spectral_noise_spike(raw, db_limit=DB_LIMIT, freq_div=TARGET_FREQ_DIV, min_freq=PSD_MIN_FREQ, target_bandwith=PSD_TARGET_BANDWIDTH) fig_individual['Power Spectral Density Noise'] = fig_noise if progress_callback: progress_callback(11) logger.info("Step 11 Completed.") @@ -5567,11 +5681,11 @@ def process_participant(file_path, progress_callback=None): fig_individual["Bad Channels Data"] = fig_raw_before if bad_channels: if BAD_CHANNELS_HANDLING == "Interpolate": - raw, fig_raw_after, fig_compare = interpolate_fNIRS_bads_weighted_average(raw, max_dist=MAX_DIST, min_neighbors=MIN_NEIGHBORS) + raw, fig_raw_after, fig_compare = interpolate_fNIRS_bads_weighted_average(raw, max_dist=MAX_DIST, min_neighbors=MIN_NEIGHBORS, short_channels_threshold=SHORT_CHANNELS_THRESHOLD) fig_individual["Data after Interpolating Bad Channels"] = fig_raw_after fig_individual["Bad Channels Interpolation Results"] = fig_compare elif BAD_CHANNELS_HANDLING == "Remove": - raw = remove_bad_channels(raw, bad_channels) + raw = remove_bad_channels(raw, bad_channels, max_bad_channels=MAX_BAD_CHANNELS) if progress_callback: progress_callback(13) logger.info("Step 13 Completed.") @@ -5592,7 +5706,7 @@ def process_participant(file_path, progress_callback=None): # Step 16: Wavelet Filtering if WAVELET and not FOLDING_BYP: - raw_od, fig = calculate_and_apply_wavelet(raw_od) + raw_od, fig = calculate_and_apply_wavelet(data=raw_od, wavelet_type=WAVELET_TYPE, wavelet_level=WAVELET_LEVEL, iqr=IQR, verbosity=VERBOSITY) fig_individual["Wavelet"] = fig if progress_callback: progress_callback(16) logger.info("Step 16 Completed.") @@ -5614,36 +5728,85 @@ def process_participant(file_path, progress_callback=None): # Step 19: Filter if FILTER and not FOLDING_BYP: - raw_haemo, fig_filter, fig_raw_haemo_filter = filter_the_data(raw_haemo) + raw_haemo, fig_filter, fig_raw_haemo_filter = filter_the_data( + raw_haemo, + filter_algorithm=FILTER_ALGORITHM, + l_freq=L_FREQ, + h_freq=H_FREQ, + l_trans_bandwidth=L_TRANS_BANDWIDTH, + h_trans_bandwidth=H_TRANS_BANDWIDTH, + iir_type=IIR_TYPE, + iir_order=IIR_ORDER, + filter_length=FILTER_LENGTH, + filter_phase=FILTER_PHASE, + fir_window=FIR_WINDOW, + fir_design=FIR_DESIGN, + iir_output=IIR_OUTPUT, + passband_ripple=PASSBAND_RIPPLE, + stopband_attenuation=STOPBAND_ATTENUATION, + filter_pad=FILTER_PAD, + skip_by_annotation=SKIP_BY_ANNOTATION, + filter_n_jobs=FILTER_N_JOBS, + verbosity=VERBOSITY + ) fig_individual["Filter_1"] = fig_filter fig_individual["Filter_2"] = fig_raw_haemo_filter if progress_callback: progress_callback(19) logger.info("Step 19 Completed.") # Step 20: Extracting Events - if not FOLDING_BYP: - events, event_dict = events_from_annotations(raw_haemo) + if EVENTS and not FOLDING_BYP: + events, event_dict = events_from_annotations(raw_haemo, event_id=EVENT_ID, regexp=EVENT_REGEX, verbose=VERBOSITY) #TODO: Implement the chunk duration fig_events = plot_events(events, event_id=event_dict, sfreq=raw_haemo.info["sfreq"], show=False) fig_individual["Events"] = fig_events if progress_callback: progress_callback(20) logger.info("Step 20 Completed.") # Step 21: Epoch Calculations - if not FOLDING_BYP: - epochs, fig_epochs = epochs_calculations(raw_haemo, events, event_dict) + if EPOCHS and EVENTS and not FOLDING_BYP: + epochs, fig_epochs = epochs_calculations( + raw_haemo, + events, + event_dict, + epoch_handling=EPOCH_HANDLING, + max_shift=MAX_SHIFT, + t_min=T_MIN, + t_max=T_MAX, + baseline=(None,0), #TODO: Unhardcode this + reject_epochs=REJECT_EPOCHS, + reject_hbo_threshold=dict(hbo=REJECT_HBO_THRESHOLD) + ) for name, fig in fig_epochs: fig_individual[f"epochs_{name}"] = fig if progress_callback: progress_callback(21) logger.info("Step 21 Completed.") # Step 22: Design Matrix - raw_haemo, df_design_matrix, fig_design_matrix = make_design_matrix(raw_haemo) # Short channel is re-applied inside this method + raw_haemo, df_design_matrix, fig_design_matrix = make_design_matrix( + raw_haemo=raw_haemo, + resample=RESAMPLE, + resample_freq=RESAMPLE_FREQ, + stim_dur=STIM_DUR, + hrf_model=HRF_MODEL, + drift_model=DRIFT_MODEL, + high_pass=HIGH_PASS, + drift_order=DRIFT_ORDER, + fir_delays=FIR_DELAYS, + min_onset=MIN_ONSET, + oversampling=OVERSAMPLING, + short_channel_regression=SHORT_CHANNEL_REGRESSION, + short_channels=SHORT_CHANNELS, + long_channels=LONG_CHANNELS, + short_channels_threshold=SHORT_CHANNELS_THRESHOLD, + long_channels_threshold=LONG_CHANNELS_THRESHOLD, + folding_bypass=FOLDING_BYP + ) fig_individual["Design Matrix"] = fig_design_matrix if progress_callback: progress_callback(22) logger.info("Step 22 Completed.") # Step 23: General Linear Model - glm_est, fig_glm_topo = make_and_run_glm(raw_haemo, df_design_matrix) + glm_est, fig_glm_topo = make_and_run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbosity=VERBOSITY) fig_individual["GLM Topography"] = fig_glm_topo if progress_callback: progress_callback(23) logger.info("23") @@ -5662,7 +5825,7 @@ def process_participant(file_path, progress_callback=None): logger.info("25") # Step 26: Generate Region of Interest Results - df_roi, fig_roi = generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path) + df_roi, fig_roi = generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path, json_location=JSON_LOCATION) fig_individual["Region of Interest"] = fig_roi if progress_callback: progress_callback(26) logger.info("26") diff --git a/main.py b/main.py index 0d3d2e2..7b5acc6 100644 --- a/main.py +++ b/main.py @@ -109,7 +109,7 @@ SECTIONS = [ {"name": "SHORT_CHANNELS", "default": True, "type": bool, "advanced": False, "help": "This should be set to True if the data has a short channel present in the data. For more information about short channels, please visit the Wiki."}, {"name": "LONG_CHANNELS", "default": True, "type": bool, "advanced": True, "help": "Should channels exceeding the maximum allowed distance be removed?"}, {"name": "SHORT_CHANNELS_THRESHOLD", "default": 0.015, "type": float, "depends_on": "SHORT_CHANNELS", "advanced": False, "help": "The maximum distance the short channel can be in metres before it is no longer considered a short channel."}, - {"name": "LONG_CHANNELS_THRESHOLD", "default": 0.045, "type": float, "advanced": False, "help": "The maximum distance channels can be in metres. Any channel longer than this distance will be discarded."}, + {"name": "LONG_CHANNELS_THRESHOLD", "default": 0.045, "type": float, "depends_on": "LONG_CHANNELS", "advanced": False, "help": "The maximum distance channels can be in metres. Any channel longer than this distance will be discarded."}, ] }, { @@ -131,7 +131,7 @@ SECTIONS = [ "title": "Scalp Coupling Index", "params": [ {"name": "SCI", "default": True, "type": bool, "advanced": False, "help": "Calculate and mark channels bad based on their Scalp Coupling Index. This metric calculates the quality of the connection between the optode and the scalp."}, - {"name": "SCI_USE_HEART_RATE_BAND", "default": True, "type": bool, "depends_on": "SCI", "advanced": False, "help": "Adjust the SCI frequency band using the participant's estimated heart rate."}, + {"name": "SCI_USE_HEART_RATE_BAND", "default": True, "type": bool, "depends_on": [{"parent_name": "SCI"}, {"parent_name": "HEART_RATE"}], "advanced": False, "help": "Adjust the SCI frequency band using the participant's estimated heart rate."}, {"name": "SCI_LOW_FREQ", "default": 0.7, "type": float, "depends_on": "SCI_USE_HEART_RATE_BAND", "depends_value": False,"advanced": True, "help": "Lower frequency cutoff for SCI bandpass filtering (Hz)."}, {"name": "SCI_HIGH_FREQ", "default": 1.5, "type": float, "depends_on": "SCI_USE_HEART_RATE_BAND", "depends_value": False, "advanced": True, "help": "Upper frequency cutoff for SCI bandpass filtering (Hz)."}, {"name": "SCI_TIME_WINDOW", "default": 3, "type": int, "depends_on": "SCI", "advanced": False, "help": "Duration of each independent SCI calculation window in seconds."}, @@ -147,7 +147,6 @@ SECTIONS = [ {"name": "SNR_SIGNAL_HIGH_FREQ", "default": 0.5, "type": float, "depends_on": "SNR", "advanced": True, "help": "Upper frequency bound for the signal band used in SNR calculation (Hz)."}, {"name": "SNR_NOISE_LOW_FREQ", "default": 1.0, "type": float, "depends_on": "SNR", "advanced": True, "help": "Lower frequency bound for the noise band used in SNR calculation (Hz)."}, {"name": "SNR_NOISE_HIGH_FREQ", "default": 10.0, "type": float, "depends_on": "SNR", "advanced": True, "help": "Upper frequency bound for the noise band used in SNR calculation (Hz)."}, - {"name": "SNR_PLOT_MAX", "default": 20.0, "type": float, "depends_on": "SNR", "advanced": True, "help": "Maximum SNR value displayed on the plot scale (dB)."}, ] }, { @@ -216,9 +215,9 @@ SECTIONS = [ "title": "Wavelet filtering", "params": [ {"name": "WAVELET", "default": True, "type": bool, "advanced": False, "help": "Apply Wavelet filtering. It is a method to filter involving decomposition, threholding, and reconstruction."}, - {"name": "IQR", "default": 1.5, "type": float, "depends_on": "WAVELET", "advanced": False, "help": "Scaling factor for the Inter-Quartile Range."}, {"name": "WAVELET_TYPE", "default": "db4", "type": str, "depends_on": "WAVELET", "advanced": False, "help": "Wavelet type. Valid values are ['bior1.1', 'bior1.3', 'bior1.5', 'bior2.2', 'bior2.4', 'bior2.6', 'bior2.8', 'bior3.1', 'bior3.3', 'bior3.5', 'bior3.7', 'bior3.9', 'bior4.4', 'bior5.5', 'bior6.8', 'coif1', 'coif2', 'coif3', 'coif4', 'coif5', 'coif6', 'coif7', 'coif8', 'coif9', 'coif10', 'coif11', 'coif12', 'coif13', 'coif14', 'coif15', 'coif16', 'coif17', 'db1', 'db2', 'db3', 'db4', 'db5', 'db6', 'db7', 'db8', 'db9', 'db10', 'db11', 'db12', 'db13', 'db14', 'db15', 'db16', 'db17', 'db18', 'db19', 'db20', 'db21', 'db22', 'db23', 'db24', 'db25', 'db26', 'db27', 'db28', 'db29', 'db30', 'db31', 'db32', 'db33', 'db34', 'db35', 'db36', 'db37', 'db38', 'dmey', 'haar', 'rbio1.1', 'rbio1.3', 'rbio1.5', 'rbio2.2', 'rbio2.4', 'rbio2.6', 'rbio2.8', 'rbio3.1', 'rbio3.3', 'rbio3.5', 'rbio3.7', 'rbio3.9', 'rbio4.4', 'rbio5.5', 'rbio6.8', 'sym2', 'sym3', 'sym4', 'sym5', 'sym6', 'sym7', 'sym8', 'sym9', 'sym10', 'sym11', 'sym12', 'sym13', 'sym14', 'sym15', 'sym16', 'sym17', 'sym18', 'sym19', 'sym20']"}, {"name": "WAVELET_LEVEL", "default": 3, "type": int, "depends_on": "WAVELET", "advanced": False, "help": "Wavelet Decomposition level (must be >= 0)."}, + {"name": "IQR", "default": 1.5, "type": float, "depends_on": "WAVELET", "advanced": False, "help": "Scaling factor for the Inter-Quartile Range."}, ] }, { @@ -239,43 +238,43 @@ SECTIONS = [ "title": "Filtering", "params": [ {"name": "FILTER", "default": True, "type": bool, "advanced": False, "help": "Should the data be bandpass filtered?"}, - {"name": "FILTER_ALGORITHM", "default": ["FIR"], "type": list, "options": ["FIR", "IIR"], "exclusive": True, "advanced": False, "help": "Filtering algorithm."}, + {"name": "FILTER_ALGORITHM", "default": ["fir"], "type": list, "options": ["fir", "iir"], "exclusive": True, "advanced": False, "help": "Filtering algorithm."}, {"name": "L_FREQ", "default": 0.005, "type": float, "depends_on": "FILTER", "advanced": False, "help": "Any frequencies lower than this value will be removed."}, {"name": "H_FREQ", "default": 0.3, "type": float, "depends_on": "FILTER", "advanced": False, "help": "Any frequencies higher than this value will be removed."}, {"name": "L_TRANS_BANDWIDTH", "default": 0.002, "type": float, "depends_on": "FILTER", "advanced": True, "help": "Width of the lower transition band to prevent abrupt filter cutoff."}, {"name": "H_TRANS_BANDWIDTH", "default": 0.002, "type": float, "depends_on": "FILTER", "advanced": True, "help": "Width of the upper transition band to prevent abrupt filter cutoff."}, - {"name": "IIR_TYPE", "default": ["butterworth"], "type": list, "options": ["butterworth", "chebyshev1", "chebyshev2", "elliptic", "bessel"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "IIR filter design."}, - {"name": "IIR_ORDER", "default": 4, "type": int, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Order of the IIR filter."}, + # {"name": "IIR_TYPE", "default": ["butterworth"], "type": list, "options": ["butterworth", "chebyshev1", "chebyshev2", "elliptic", "bessel"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "IIR filter design."}, + # {"name": "IIR_ORDER", "default": 4, "type": int, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Order of the IIR filter."}, {"name": "FILTER_LENGTH", "default": "auto", "type": str, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Length of the FIR filter. 'auto' allows automatic selection."}, - {"name": "FILTER_PHASE", "default": ["zero"], "type": list, "options": ["zero", "zero-double", "minimum", "minimum-half"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Phase response of the FIR filter."}, + {"name": "FILTER_PHASE", "default": ["zero"], "type": list, "options": ["zero", "zero-double", "minimum", "minimum-half", "linear"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Phase response of the FIR filter."}, {"name": "FIR_WINDOW", "default": ["hamming"], "type": list, "options": ["hamming", "hann", "blackman"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Window function used when designing the FIR filter."}, {"name": "FIR_DESIGN", "default": ["firwin"], "type": list, "options": ["firwin", "firwin2"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Method used to design the FIR filter."}, - {"name": "IIR_OUTPUT", "default": ["sos"], "type": list, "options": ["sos", "ba", "zpk"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Representation used for IIR filter coefficients."}, - {"name": "PASSBAND_RIPPLE", "default": 1.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev1", "elliptic"], "advanced": True, "help": "Maximum allowed ripple in the passband (dB)."}, - {"name": "STOPBAND_ATTENUATION", "default": 40.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev2", "elliptic"], "advanced": True, "help": "Minimum attenuation in the stopband (dB)."}, + # {"name": "IIR_OUTPUT", "default": ["sos"], "type": list, "options": ["sos", "ba", "zpk"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Representation used for IIR filter coefficients."}, + # {"name": "PASSBAND_RIPPLE", "default": 1.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev1", "elliptic"], "advanced": True, "help": "Maximum allowed ripple in the passband (dB)."}, + # {"name": "STOPBAND_ATTENUATION", "default": 40.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev2", "elliptic"], "advanced": True, "help": "Minimum attenuation in the stopband (dB)."}, {"name": "FILTER_PAD", "default": ["reflect_limited"], "type": list, "options": ["reflect_limited", "reflect", "edge", "constant"], "exclusive": True, "depends_on": "FILTER", "advanced": True, "help": "Padding strategy used during filtering to reduce edge artifacts."}, - {"name": "SKIP_BY_ANNOTATION", "default": ["edge", "bad_acq_skip"], "type": list, "depends_on": "FILTER", "advanced": True, "help": "Annotations that should be skipped when applying the filter."}, + # {"name": "SKIP_BY_ANNOTATION", "default": ["edge", "bad_acq_skip"], "type": list, "depends_on": "FILTER", "advanced": True, "help": "Annotations that should be skipped when applying the filter."}, {"name": "FILTER_N_JOBS", "default": 1, "type": int, "advanced": True, "help": "Number of parallel jobs used during filtering. Use -1 to use all available CPUs."}, ] }, { - "title": "Extracting Events*", + "title": "Extracting Events", "params": [ - {"name": "EVENTS", "default": True, "type": bool, "advanced": False, "help": "Extract events from annotations for visualization and downstream event-based analysis."}, + {"name": "EVENTS", "default": True, "type": bool, "advanced": True, "help": "Extract events from annotations for visualization and downstream event-based analysis."}, {"name": "EVENT_ID", "default": "auto", "type": str, "advanced": True, "help": "Controls how annotation descriptions are converted into event identifiers. Use 'auto' for automatic event detection."}, - {"name": "EVENT_REGEX", "default": "^(?![Bb][Aa][Dd]|[Ee][Dd][Gg][Ee]).*$", "type": str, "advanced": True, "help": "Regular expression used to select which annotations are converted into events. By default, bad and edge annotations are ignored."}, - {"name": "EVENT_CHUNK_DURATION", "default": None, "type": float, "advanced": True, "help": "If provided, creates repeated events at this interval within longer annotations instead of only using annotation onset times."}, + {"name": "EVENT_REGEX", "default": r"^(?![Bb][Aa][Dd]|[Ee][Dd][Gg][Ee]).*$", "type": str, "advanced": True, "help": "Regular expression used to select which annotations are converted into events. By default, bad and edge annotations are ignored."}, + # {"name": "EVENT_CHUNK_DURATION", "default": 0.0, "type": float, "advanced": True, "help": "If provided, creates repeated events at this interval within longer annotations instead of only using annotation onset times."}, ] }, { "title": "Epoch Calculations", "params": [ - {"name": "EPOCHS", "default": True, "type": bool, "advanced": False, "help": "Create epochs around extracted events for condition-based analysis."}, + {"name": "EPOCHS", "default": True, "type": bool, "depends_on": "EVENTS", "advanced": True, "help": "Create epochs around extracted events for condition-based analysis."}, {"name": "EPOCH_HANDLING", "default": ["shift"], "type": list, "options": ["shift", "strict", "drop"], "exclusive": True, "advanced": False, "help": "How to handle events occurring at the same sample. Shift moves conflicting events forward, strict raises an error, and drop removes conflicting events."}, {"name": "MAX_SHIFT", "default": 5, "type": int, "depends_on": "EPOCH_HANDLING", "depends_value": "shift", "advanced": True, "help": "Maximum number of samples to shift conflicting events before failing."}, {"name": "T_MIN", "default": -5.0, "type": float, "advanced": False, "help": "Time in seconds before each event to include in the epoch."}, {"name": "T_MAX", "default": 15.0, "type": float, "advanced": False, "help": "Time in seconds after each event to include in the epoch."}, - {"name": "BASELINE", "default": ["pre_event"], "type": list, "options": ["none", "pre_event"], "exclusive": True, "advanced": False, "help": "Baseline correction applied to epochs. Pre-event uses the period before the event as baseline."}, + # {"name": "BASELINE", "default": ["pre_event"], "type": list, "options": ["none", "pre_event"], "exclusive": True, "advanced": False, "help": "Baseline correction applied to epochs. Pre-event uses the period before the event as baseline."}, {"name": "REJECT_EPOCHS", "default": True, "type": bool, "advanced": False, "help": "Automatically reject epochs containing excessively large haemoglobin amplitude changes."}, {"name": "REJECT_HBO_THRESHOLD", "default": 80e-7, "type": float, "depends_on": "REJECT_EPOCHS", "advanced": True, "help": "Maximum allowed HbO amplitude before an epoch is rejected."}, ] @@ -310,13 +309,6 @@ SECTIONS = [ {"name": "JSON_LOCATION", "default": "", "type": "json_file", "advanced": False, "help": "Location of the JSON file containing region of interest results for significance calculations."}, ] }, - { - "title": "Contrast", - "params": [ - {"name": "CONTRAST_BASELINE", "default": True, "type": bool, "advanced": True, "help": "Calculate contrasts comparing each condition against zero/baseline. This may not always represent a meaningful physiological baseline in fNIRS."}, - {"name": "CONTRAST_PAIRWISE", "default": True, "type": bool, "advanced": True, "help": "Calculate contrasts comparing each experimental condition against every other condition. The number of contrasts increases with the number of conditions."} - ] - }, { "title": "Other", "params": [ @@ -1445,7 +1437,7 @@ class MainApplication(QMainWindow): if self.advanced_parameters: self.update_sections(0) - + if hasattr(self, 'recent_files_menu'): self.update_recent_files_menu() diff --git a/src/shared/flaresbasewidget.py b/src/shared/flaresbasewidget.py index e351103..1d3c1ee 100644 --- a/src/shared/flaresbasewidget.py +++ b/src/shared/flaresbasewidget.py @@ -417,10 +417,20 @@ class ParamSection(QWidget): widget.textChanged.connect(lambda val, p=param["name"]: self.check_if_changed(p, val)) if "depends_on" in param: + deps_list = [] + + if isinstance(param["depends_on"], list): + deps_list = param["depends_on"] + + else: + deps_list.append({ + "parent_name": param["depends_on"], + "depends_value": param.get("depends_value", "True") + }) + self.dependencies.append({ "child_name": param["name"], - "parent_name": param["depends_on"], - "depends_value": param.get("depends_value", "True") + "conditions": deps_list }) widget.setToolTip(help_text) @@ -532,23 +542,44 @@ class ParamSection(QWidget): """Disables/Enables widgets based on parent selection values.""" for dep in self.dependencies: child_info = self.widgets.get(dep["child_name"]) - parent_info = self.widgets.get(dep["parent_name"]) + if not child_info: + continue + + # Default to enabled until a condition fails + all_conditions_met = True + + for cond in dep["conditions"]: + parent_name = cond.get("parent_name") or cond.get("parent") + required_value = str(cond.get("depends_value") if "depends_value" in cond else cond.get("value", "True")) + + parent_info = self.widgets.get(parent_name) + if not parent_info: + all_conditions_met = False + break - if child_info and parent_info: parent_widget = parent_info["widget"] + + # Extract current parent value based on widget type + if isinstance(parent_widget, QComboBox): + current_parent_value = parent_widget.currentText() + elif isinstance(parent_widget, QLineEdit): + current_parent_value = parent_widget.text() + elif isinstance(parent_widget, QSpinBox): + current_parent_value = str(parent_widget.value()) + else: + current_parent_value = str(parent_widget) + + # If any condition fails, flag as false + if current_parent_value != required_value: + all_conditions_met = False + break - # Get current value of parent (works for both bool-combos and list-combos) - current_parent_value = parent_widget.currentText() - - # Check if it matches the required value - is_active = (current_parent_value == dep["depends_value"]) - - # Toggle the entire row (Button, Label, and Input) - h_layout = child_info["h_layout"] - for i in range(h_layout.count()): - item = h_layout.itemAt(i).widget() - if item: - item.setEnabled(is_active) + # Toggle the entire row (Button, Label, and Input) + h_layout = child_info["h_layout"] + for i in range(h_layout.count()): + item = h_layout.itemAt(i).widget() + if item: + item.setEnabled(all_conditions_met) def _create_multiselect_dropdown(self, items): combo = FullClickComboBox()