Skip to content

Commit

Permalink
Fix and improve deglitching. Add it in the surface- and template-base…
Browse files Browse the repository at this point in the history 
…d full pipeline, as an optional step.
  • Loading branch information
@thomas-vincent
thomas-vincent committed Feb 21, 2019
1 parent e911370 commit 4f4f1da
Show file tree
Hide file tree
Showing 3 changed files with 73 additions and 28 deletions.
23 changes: 19 additions & 4 deletions bst_plugin/nst_ppl_surface_template_V1.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,11 +200,11 @@

% Run preprocessings
[sFiles_preprocessed, hb_types, redone_preprocs, preproc_folder] = preprocs(file_raw, sFile_raw_fhm, options, force_redo|fhm_redone);

% Run 1st level GLM
[sFiles_GLM, sFiles_con, redone_any_contrast, glm_folder] = glm_1st_level(sFiles_preprocessed, options, ...
redone_preprocs | force_redo);

if isubject==1
all_sFiles_con = cell(size(sFiles_con, 1), size(sFiles_con, 2), length(subject_names));
end
Expand Down Expand Up @@ -445,10 +445,21 @@

% TODO: export bad channel tagging information

% Deglitching
if options.deglitch.do
redo_parent = force_redo | options.deglitch.redo;
sFile_deglitched = nst_run_bst_proc([preproc_folder 'Deglitched'], redo_parent, ...
'process_nst_deglitch', sFile_raw, [], ...
'factor_std_grad', options.deglitch.agrad_std_factor);
else
redo_parent = force_redo;
sFile_deglitched = sFile_raw;
end

% Motion correction
redo_parent = force_redo | options.moco.redo;
redo_parent = redo_parent | options.moco.redo;
[sFileMoco, redo_parent] = nst_run_bst_proc([preproc_folder 'Motion-corrected'], redo_parent, ...
'process_nst_motion_correction', sFile_raw, [], ...
'process_nst_motion_correction', sFile_deglitched, [], ...
'option_event_name', 'movement_artefacts');

% Resample to 5Hz (save some space)
Expand Down Expand Up @@ -675,6 +686,10 @@

options.head_model.redo = 0;

options.deglitch.do = 0;
options.deglitch.redo = 0;
options.deglitch.agrad_std_factor = 2.5;

options.moco.redo = 0;
options.moco.export_dir = fullfile('.', 'moco_marking');

Expand Down
38 changes: 30 additions & 8 deletions bst_plugin/process_nst_deglitch.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,15 @@
sProcess.InputTypes = {'raw', 'data'};
sProcess.OutputTypes = {'raw', 'data'};

% TODO: option to only tag as events
sProcess.options.factor_std_grad.Comment = 'Variation threshold: ';
sProcess.options.factor_std_grad.Type = 'value';
sProcess.options.factor_std_grad.Value = {2.5, '* std(gradient)', 2}; % last if nb of subunit digits, use 0 for integer

% sProcess.options.tag_only.Comment = 'Only mark glitches as events';
% sProcess.options.tag_only.Type = 'checkbox';
% sProcess.options.tag_only.Value = 0;

%TODO: apply only to NIRS channels
end

%% ===== FORMAT COMMENT =====
Expand All @@ -49,16 +57,18 @@

%% ===== RUN =====
function sInputs = Run(sProcess, sInputs) %#ok<DEFNU>
sInputs.A = Compute(sInputs.A, sProcess.options.factor_std_grad.Value{1});
end

sInputs.A = Compute(sInputs.A);
%% ===== Compute =====
function [nirs_deglitched, glitch_flags] = Compute(nirs_sig, std_factor_thresh)

% TODO: manage detecting only
if nargin < 2
std_factor_thresh = 2.5;
end

%% ===== Compute =====
function [nirs_deglitched, glitch_flags] = Compute(nirs_sig)
[nb_positions, nb_samples] = size(nirs_sig);
glitch_flags = detect_glitches(nirs_sig);
glitch_flags = detect_glitches(nirs_sig, std_factor_thresh);
nirs_deglitched = nirs_sig;
for ipos=1:nb_positions
nirs_deglitched(ipos, glitch_flags(ipos, :)) = repmat(mean(nirs_sig(ipos, :)), ...
Expand Down Expand Up @@ -91,9 +101,21 @@
signal_tmp = [nirs_sig(:, 2) nirs_sig nirs_sig(:, end-1)]; % mirror edges
grad = diff(signal_tmp, 1, 2);
abs_grad = abs(grad);
std_agrad = std(abs_grad, 0, 2);

% Robust standard deviation of absolute gradient
% -> clip values within 1% to 99% of total range, to remove extreme values
% then compute reference std of absolute gradient over these clipped values

sorted_signal = sort(signal_tmp, 2);
thresh_low = sorted_signal(:, round(nb_samples*0.01));
thresh_high = sorted_signal(:, round(nb_samples*0.99));
for ipos=1:size(signal_tmp, 1)
signal_tmp(ipos, signal_tmp(ipos, :) < thresh_low(ipos)) = thresh_low(ipos);
signal_tmp(ipos, signal_tmp(ipos, :) > thresh_high(ipos)) = thresh_high(ipos);
end
std_agrad = std(abs(diff(signal_tmp, 1, 2)), 0, 2);

glitch_canditates = [zeros(nb_pos, 1) (abs_grad > std_factor_thresh * std_agrad) .* sign(grad)];
glitch_flags = [abs(diff(glitch_canditates, 1, 2))==2 false(nb_pos, 1)] & glitch_canditates;
glitch_flags = glitch_flags(:, 2:end-1); % remove mirrored edges

end
40 changes: 24 additions & 16 deletions test/DeglitchTest.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,27 @@ function tear_down(testCase)

methods(Test)

function test_detection(testCase)

% Request utest data -> time-series + ground-truth marking
deglitch_fn = utest_request_data({'artifacts','glitch_data.mat'});
deglitch_data = load(deglitch_fn);

glitch_flags = process_nst_deglitch('detect_glitches', deglitch_data.nirs_signal);

if 0
for ipos=1:size(deglitch_data.nirs_signal, 1)
figure(); hold on;
plot(deglitch_data.nirs_signal(ipos, :), 'g', 'LineWidth', 3);
glitches_idx = find(glitch_flags(ipos, :));
plot(glitches_idx, deglitch_data.nirs_signal(ipos,glitches_idx), 'or', 'MarkerSize', 8);
end
end

% Check against ground-truth markings
testCase.assertTrue(all(glitch_flags(:)==deglitch_data.glitch_flags(:)));
end

function test_deglitch(testCase)

% Request utest data -> time-series + ground-truth marking
Expand Down Expand Up @@ -52,27 +73,14 @@ function test_deglitch(testCase)

if 0
figure(); hold on;
plot(deglitch_data.nirs_signal(ipos, :), 'k', 'LineWidth', 3);
plot(deglitch_data.nirs_signal(ipos, :), 'g', 'LineWidth', 3);
plot(signal_clean, 'b', 'LineWidth', 2);
plot(deglitched_signal(ipos, :), 'r');
plot(glitches_idx, deglitch_data.nirs_signal(ipos,glitches_idx), 'or', 'MarkerSize', 8);
end

testCase.assertTrue(all_close(signal_clean, deglitched_signal(ipos, :)));
end
end

function test_detection(testCase)

% Request utest data -> time-series + ground-truth marking
deglitch_fn = utest_request_data({'artifacts','glitch_data.mat'});
deglitch_data = load(deglitch_fn);

glitch_flags = process_nst_deglitch('detect_glitches', deglitch_data.nirs_signal);

% Check against ground-truth markings
testCase.assertTrue(all(glitch_flags(:)==deglitch_data.glitch_flags(:)));
end


end
end
end

0 comments on commit 4f4f1da

Please sign in to comment.