Skip to content

Commit

Permalink
Merge pull request #85 from Nirstorm/surface_glm
Browse files Browse the repository at this point in the history 
Surface glm pipeline - glitch removal
  • Loading branch information
@thomas-vincent
thomas-vincent committed Feb 21, 2019
2 parents cdcfec2 + 4f4f1da commit 4fcaf1b
Show file tree
Hide file tree
Showing 5 changed files with 279 additions and 4 deletions.
1 change: 1 addition & 0 deletions bst_plugin/MANIFEST
View file
@@ -1,3 +1,4 @@
process_nst_deglitch.m
process_nst_prefix_matrix.m
nst_save_table_in_bst.m
nst_parse_bst_item_name.m
Expand Down
23 changes: 19 additions & 4 deletions bst_plugin/nst_ppl_surface_template_V1.m
View file
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
121 changes: 121 additions & 0 deletions bst_plugin/process_nst_deglitch.m
View file
@@ -0,0 +1,121 @@
function varargout = process_nst_deglitch( varargin )

% @=============================================================================
% This software is part of the Brainstorm software:
% http://neuroimage.usc.edu/brainstorm
%
% Copyright (c)2000-2013 Brainstorm by the University of Southern California
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Thomas Vincent (2019)

eval(macro_method);
end


%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>

% Description the process
sProcess.Comment = 'Remove glitches';
sProcess.Category = 'Filter';
sProcess.SubGroup = 'NIRS';
sProcess.Index = 1403;
sProcess.isSeparator = 0;
sProcess.Description = 'https://github.com/Nirstorm/nirstorm/wiki/Remove-glitches';

% Definition of the input accepted by this process
sProcess.InputTypes = {'raw', 'data'};
sProcess.OutputTypes = {'raw', 'data'};

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 =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
Comment = sProcess.Comment;
end

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

%% ===== Compute =====
function [nirs_deglitched, glitch_flags] = Compute(nirs_sig, std_factor_thresh)

if nargin < 2
std_factor_thresh = 2.5;
end

[nb_positions, nb_samples] = size(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, :)), ...
1, sum(glitch_flags(ipos, :)));
if glitch_flags(ipos, end)
nirs_deglitched(ipos, end) = mean(nirs_sig(ipos, nb_samples-2:nb_samples-1));
glitch_flags(ipos, end) = 0;
end
if glitch_flags(ipos, 1)
nirs_deglitched(ipos, 1) = mean(nirs_sig(ipos, 2:3));
glitch_flags(ipos, 1) = 0;
end
% Maybe the following can be vectorized but there was something strange with
% find function on 2D array
glitch_i = find(glitch_flags(ipos, :));
nirs_deglitched(ipos, glitch_i) = (nirs_sig(ipos, glitch_i-1) + nirs_sig(ipos, glitch_i+1)) ./ 2;
end
end

function glitch_flags = detect_glitches(nirs_sig, std_factor_thresh)

% nirs_sig: (nb_positions, nb_samples)

[nb_pos, nb_samples] = size(nirs_sig);

if nargin < 2
std_factor_thresh = 2.5;
end

signal_tmp = [nirs_sig(:, 2) nirs_sig nirs_sig(:, end-1)]; % mirror edges
grad = diff(signal_tmp, 1, 2);
abs_grad = abs(grad);

% 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
52 changes: 52 additions & 0 deletions scripts/deglitch_example.m
View file
@@ -0,0 +1,52 @@
function deglitch_example()
%% Generate signal with glitches
nb_samples = 500;
nb_glitches = 10;
nb_draws = 2;
signal = randn(nb_samples, nb_draws);
gen_spikes = [ones(1, nb_draws) ; randi([3,nb_samples-2], nb_glitches-2, nb_draws) ; zeros(1, nb_draws) + nb_samples];
for idraw=1:nb_draws
signal(gen_spikes(:, idraw), idraw) = rand(nb_glitches, 1) * 5 + 5;
end

%% Detect glitches
std_factor_thresh = 2.5;
signal_tmp = [signal(2, :) ; signal ; signal(end-1, :)]; % mirror edges
grad = diff(signal_tmp);
abs_grad = abs(grad);
std_agrad = std(abs_grad);
glitch_canditates = [zeros(1, nb_draws) ; (abs_grad > std_factor_thresh * std_agrad) .* sign(grad)];
glitch_flags = [abs(diff(glitch_canditates))==2 ; false(1, nb_draws)] & glitch_canditates;
glitch_flags = glitch_flags(2:end-1, :); % remove mirrored edges

%% Clean glitches
signal_cleaned = signal;
for idraw=1:nb_draws
signal_cleaned(glitch_flags(:, idraw), idraw) = repmat(mean(signal(:,idraw)), sum(glitch_flags(:, idraw)), 1);
if glitch_flags(end, idraw)
signal_cleaned(end, idraw) = mean(signal(nb_samples-2:nb_samples-1, idraw));
glitch_flags(end, idraw) = 0;
end
if glitch_flags(1, idraw)
signal_cleaned(1, idraw) = mean(signal(2:3, idraw));
glitch_flags(1, idraw) = 0;
end
% Maybe the following can be vectorized but there was smthg strange with
% find on 2D array
glitch_i = find(glitch_flags(:, idraw));
signal_cleaned(glitch_i, idraw) = (signal(glitch_i-1, idraw) + signal(glitch_i+1, idraw)) ./ 2;
end

%% Plot results
figure(); hold on;
plot(signal(:, 1), 'b', 'Linewidth', 3);
plot(find(glitch_flags(:,1)), signal(glitch_flags(:,1), 1), 'Xg', 'MarkerSize', 10);
plot(signal_cleaned(:, 1), 'r', 'Linewidth', 2);
plot(gen_spikes(:, 1), signal(gen_spikes(:, 1), 1), 'ok', 'MarkerSize', 10);

figure(); hold on;
plot(signal(:, 2), 'b', 'Linewidth', 3);
plot(find(glitch_flags(:,2)), signal(glitch_flags(:,2), 2), 'Xg', 'MarkerSize', 10);
plot(signal_cleaned(:, 2), 'r', 'Linewidth', 2);
plot(gen_spikes(:, 2), signal(gen_spikes(:, 2), 2), 'ok', 'MarkerSize', 10);
end
86 changes: 86 additions & 0 deletions test/DeglitchTest.m
View file
@@ -0,0 +1,86 @@
classdef DeglitchTest < matlab.unittest.TestCase

properties
tmp_dir
end

methods(TestMethodSetup)
function setup(testCase)
tmpd = tempname;
mkdir(tmpd);
testCase.tmp_dir = tmpd;
utest_bst_setup();
end
end

methods(TestMethodTeardown)
function tear_down(testCase)
rmdir(testCase.tmp_dir, 's');
utest_clean_bst();
end
end

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
deglitch_fn = utest_request_data({'artifacts','glitch_data.mat'});
deglitch_data = load(deglitch_fn);

% Import in brainstorm
% subject_name = bst_create_test_subject('');
deglitched_signal = process_nst_deglitch('Compute', deglitch_data.nirs_signal);
if 0
for ipos=1:size(deglitch_data.nirs_signal, 1)
figure(); hold on;
plot(deglitch_data.nirs_signal(ipos, :), 'b', 'LineWidth', 3);
plot(deglitched_signal(ipos, :), 'r');
end
end

% Check result
nb_samples = size(deglitch_data.nirs_signal, 2);
for ipos=1:size(deglitch_data.nirs_signal, 1)
gflags = deglitch_data.glitch_flags(ipos, :);
gflags([1 end]) = 0;
signal_clean = deglitch_data.nirs_signal(ipos, :);
signal_clean(1) = mean(signal_clean(2:3));
signal_clean(end) = mean(signal_clean(nb_samples-2:nb_samples-1));
glitches_idx = find(gflags);
signal_clean(gflags) = (signal_clean(glitches_idx-1) + signal_clean(glitches_idx+1)) / 2;

if 0
figure(); hold on;
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
end
end

0 comments on commit 4fcaf1b

Please sign in to comment.