From 24f06cfbbb7a4a3ab02931dc2a1e1f0cf17500cb Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Sat, 18 May 2019 21:43:53 -0400 Subject: [PATCH 01/12] Minor improvements on ppl minimal example script --- scripts/surface_template_full_group_pipeline_V1.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/surface_template_full_group_pipeline_V1.m b/scripts/surface_template_full_group_pipeline_V1.m index 7c4dd54e..c4cecd9a 100644 --- a/scripts/surface_template_full_group_pipeline_V1.m +++ b/scripts/surface_template_full_group_pipeline_V1.m @@ -52,12 +52,12 @@ function surface_template_full_group_pipeline_V1() %% Import data options = nst_ppl_surface_template_V1('get_options'); % get default pipeline options -[sFiles, imported] = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); +sFiles = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); % Read stimulation events from AUX channel for ifile=1:length(sFiles) - if imported(ifile) - % Read events from aux channel + evt_data = load(file_fullpath(sFiles{ifile}), 'Events'); + if ~any(strcmp({evt_data.Events.label}, 'motor')) % Insure that events were not already loaded bst_process('CallProcess', 'process_evt_read', sFiles{ifile}, [], ... 'stimchan', 'NIRS_AUX', ... 'trackmode', 3, ... % Value: detect the changes of channel value From 8e81e033a6633fcef02ebe4704e4daa501ddcf74 Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Sat, 18 May 2019 21:44:39 -0400 Subject: [PATCH 02/12] Rename utest file --- test/{SurfTplPipelineImportTest.m => SurfTplPipelineTest.m} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/{SurfTplPipelineImportTest.m => SurfTplPipelineTest.m} (100%) diff --git a/test/SurfTplPipelineImportTest.m b/test/SurfTplPipelineTest.m similarity index 100% rename from test/SurfTplPipelineImportTest.m rename to test/SurfTplPipelineTest.m From 21aafb99adc32f4cbc097062ff06703fb89bf500 Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Tue, 21 May 2019 21:27:22 -0400 Subject: [PATCH 03/12] Add wrapper functions to load/save events. To insure backward compatibility. To control default bst event format used in nirstorm. --- bst_plugin/MANIFEST | 2 + bst_plugin/nst_bst_load_events.m | 15 ++++ bst_plugin/nst_bst_save_events.m | 15 ++++ test/BstEventIOTest.m | 150 +++++++++++++++++++++++++++++++ 4 files changed, 182 insertions(+) create mode 100644 bst_plugin/nst_bst_load_events.m create mode 100644 bst_plugin/nst_bst_save_events.m create mode 100644 test/BstEventIOTest.m diff --git a/bst_plugin/MANIFEST b/bst_plugin/MANIFEST index 4dc4c6d7..1848af75 100644 --- a/bst_plugin/MANIFEST +++ b/bst_plugin/MANIFEST @@ -1,3 +1,5 @@ +nst_bst_save_events.m +nst_bst_load_events.m nst_io_fetch_sample_data.m nst_bst_set_template_anatomy.m nst_core_get_available_templates.m diff --git a/bst_plugin/nst_bst_load_events.m b/bst_plugin/nst_bst_load_events.m new file mode 100644 index 00000000..82d7cf5b --- /dev/null +++ b/bst_plugin/nst_bst_load_events.m @@ -0,0 +1,15 @@ +function bst_events = nst_bst_load_events(event_fn) +% Wrapper around bst functions to save events in a backward compatible mode +% Save according to format based on extension +wrapper.events = []; +[root, bfn, ext] = fileparts(event_fn); +switch ext + case '.mat' + % Use dummy sampling frequency high enough so that times are not rounded + wrapper.prop.sfreq = 1e6; + wrapper = import_events(wrapper, [], event_fn, 'BST'); + bst_events = wrapper.events; + otherwise + error(sprintf('Unsupported extension for event file: "%s"', ext)); +end +end \ No newline at end of file diff --git a/bst_plugin/nst_bst_save_events.m b/bst_plugin/nst_bst_save_events.m new file mode 100644 index 00000000..92f323d5 --- /dev/null +++ b/bst_plugin/nst_bst_save_events.m @@ -0,0 +1,15 @@ +function nst_bst_save_events(bst_events, event_fn) +% Wrapper around bst functions to load event in a backward compatible mode +% Load according to format based on extension + +[root, bfn, ext] = fileparts(event_fn); +switch ext + case '.mat' + wrapper.events = bst_events; + % Emulate sampling frequency high enough to avoid rounding of times + % wrapper.prop.sfreq = 2.0 / min(diff(sort([bst_events.times]))); + export_events(wrapper, [], event_fn); + otherwise + error(sprintf('Unsupported extension for event file: "%s"', ext)); +end +end \ No newline at end of file diff --git a/test/BstEventIOTest.m b/test/BstEventIOTest.m new file mode 100644 index 00000000..6de81aaf --- /dev/null +++ b/test/BstEventIOTest.m @@ -0,0 +1,150 @@ +classdef BstEventIOTest < 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_load_save_current_version(testCase) + bst_events = get_test_events(); + event_fn = fullfile(testCase.tmp_dir, 'bst_events.mat'); + nst_bst_save_events(bst_events, event_fn); + testCase.assertTrue(exist(event_fn, 'file')==2); + loaded_events = nst_bst_load_events(event_fn); + assert_struct_are_equal(bst_events, loaded_events); + end + + function test_load_save_backward_compatibility(testCase) + bst_events = get_test_events(); + % Dump event file and check that it's available online in the + % unittest nirstorm data repository + forge_current_utest_data_file(bst_events); + + % Loop through all available previous versions of dumped event + % file: load them and check that they are still consistent with + % current test events. + utest_data_dir = get_utest_data_dir(); + flist = dir(utest_data_dir); + for ifile=1:length(flist) + if ~strcmp(flist(ifile).name, '.') && ~strcmp(flist(ifile).name, '..') + loaded_events = nst_bst_load_events(fullfile(utest_data_dir, flist(ifile).name)); + assert_struct_are_equal(bst_events, loaded_events); + end + end + end + + % Error tests: non-available channels, incompatible time axes, + % Assert warning is issued when rounding of events make them + % differ from what was stored on the drive. + end +end + +function bst_event = get_test_events() +% Return a BST event struct with corner cases: +% - multiple conditions +% - some have empty notes / channels tags +% - some have non-empty notes with weird chars +% - some selected, some not +% +% Maintain this function everytime bst event structure changes +% -> to add new fields use expected default values so that backward compatibility +% is guarenteed. test_load_save_backward_compatibility should yell at +% you if not. +bst_event = db_template('event'); +bst_event(1).label = 'condition1'; +bst_event(1).color = [0., 0., 0.]; +bst_event(1).epochs = [1]; +bst_event(1).times = [10.3]; +bst_event(1).reactTimes = []; +bst_event(1).select = 1; +bst_event(1).channels = {{}}; +bst_event(1).notes = {''}; + +bst_event(2).label = 'condition2'; +bst_event(2).color = [0., 1., 0.3]; +bst_event(2).epochs = [1 1 1]; +bst_event(2).times = [1.3 7.5 12.6]; +bst_event(2).reactTimes = []; +bst_event(2).select = 0; +bst_event(2).channels = {{'S1D1', 'S2D3'}, {}, {'S121', 'S4D3'}}; +bst_event(2).notes = {sprintf('\t!:/\\,#%%*+{}[]wazaéàè?`''&$'), ... + '', '12345678910-=;,1,2||'}; + +end + +function forge_current_utest_data_file(bst_events) +% Save given bst events to file tagged with current bst version, if not +% already available. +% Insure that this file is available online. + +bst_version = bst_get('Version'); +event_bfn = sprintf('bst_events_v%s.mat', bst_version.Version); +event_fn = fullfile(get_utest_data_dir(), event_bfn); +if exist(event_fn, 'file') ~= 2 + nst_bst_save_events(bst_events, event_fn); +end +assert(exist(event_fn, 'file')==2); + +event_fn_url = [get_utest_data_url() '/' event_bfn]; +% Check if remote file exists: +jurl = java.net.URL(event_fn_url); +conn = openConnection(jurl); +conn.setConnectTimeout(15000); +status = getResponseCode(conn); +if status == 404 + error('Remote utest event data file not found: %s', event_fn_url); +end +end + +function utest_data_dir = get_utest_data_dir() +sdirs = get_utest_data_subdirs(); +utest_data_dir = fullfile(nst_get_local_user_dir(), sdirs{:}); +if ~exist(utest_data_dir, 'dir') + mkdir(utest_data_dir); +end +end + +function utest_data_url = get_utest_data_url() +sdirs = get_utest_data_subdirs(); +utest_data_url = [nst_get_repository_url() '/' strjoin(sdirs, '/')]; +end + +function sdirs = get_utest_data_subdirs() +sdirs = {'unittest', 'bst_events'}; +end + +function assert_struct_are_equal(s1, s2) + +assert(length(s1) == length(s2)); +f1 = fieldnames(s1); +f2 = fieldnames(s2); +assert(length(f1) == length(f2)); +assert(all(ismember(f1, f2))); +for i_item=1:length(s1) + for ifield=1:length(f1) + f = f1{ifield}; + if ~isequaln(s1(i_item).(f), s2(i_item).(f)) + error('Unequal values for field %s', f); + end + end +end +end \ No newline at end of file From 85f189110da3641621f89ede1c9df18df83ef9ea Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Mon, 27 May 2019 14:02:41 -0400 Subject: [PATCH 04/12] Add IO function for markings --- bst_plugin/MANIFEST | 6 ++- bst_plugin/nst_bst_export_channel_flags.m | 19 +++++++++ bst_plugin/nst_bst_export_events.m | 17 ++++++++ bst_plugin/nst_bst_import_channel_flags.m | 19 +++++++++ bst_plugin/nst_bst_import_events.m | 49 +++++++++++++++++++++++ test/BstEventIOTest.m | 40 ++---------------- 6 files changed, 111 insertions(+), 39 deletions(-) create mode 100644 bst_plugin/nst_bst_export_channel_flags.m create mode 100644 bst_plugin/nst_bst_export_events.m create mode 100644 bst_plugin/nst_bst_import_channel_flags.m create mode 100644 bst_plugin/nst_bst_import_events.m diff --git a/bst_plugin/MANIFEST b/bst_plugin/MANIFEST index 1848af75..81da6b1b 100644 --- a/bst_plugin/MANIFEST +++ b/bst_plugin/MANIFEST @@ -1,5 +1,7 @@ -nst_bst_save_events.m -nst_bst_load_events.m +nst_bst_import_channel_flags.m +nst_bst_export_channel_flags.m +nst_bst_export_events.m +nst_bst_import_events.m nst_io_fetch_sample_data.m nst_bst_set_template_anatomy.m nst_core_get_available_templates.m diff --git a/bst_plugin/nst_bst_export_channel_flags.m b/bst_plugin/nst_bst_export_channel_flags.m new file mode 100644 index 00000000..6ae88cec --- /dev/null +++ b/bst_plugin/nst_bst_export_channel_flags.m @@ -0,0 +1,19 @@ +function nst_bst_export_channel_flags(sFile, channel_flags_fn) +%NST_BST_EXPORT_CHANNEL_FLAGS Save channel flags from given sFile +% into given channel_flags_fn. +% TODO: also save list of channels +if ischar(sFile) + bst_chan_data = load(file_fullpath(sFile), 'ChannelFlag'); + bst_chan_flags = bst_chan_data.ChannelFlag; +end +[root, bfn, ext] = fileparts(channel_flags_fn); +switch ext + case '.mat' + s.channel_flags = bst_chan_flags; + bst_save(channel_flags_fn, s, 'v7'); + otherwise + error(sprintf('Unsupported extension for channel file: "%s"', ext)); +end +end + + diff --git a/bst_plugin/nst_bst_export_events.m b/bst_plugin/nst_bst_export_events.m new file mode 100644 index 00000000..da79a994 --- /dev/null +++ b/bst_plugin/nst_bst_export_events.m @@ -0,0 +1,17 @@ +function nst_bst_export_events(sFile, bst_events, event_fn) +% Wrapper around bst functions to load event in a backward compatible mode +% Load according to format based on extension + +if ischar(bst_events) + bst_evt_data = load(file_fullpath(bst_events), 'Events'); + bst_events = bst_evt_data.Events; +end +[root, bfn, ext] = fileparts(event_fn); +switch ext + case '.mat' + s.events = bst_events; + bst_save(event_fn, s, 'v7'); + otherwise + error(sprintf('Unsupported extension for event file: "%s"', ext)); +end +end \ No newline at end of file diff --git a/bst_plugin/nst_bst_import_channel_flags.m b/bst_plugin/nst_bst_import_channel_flags.m new file mode 100644 index 00000000..6f544f86 --- /dev/null +++ b/bst_plugin/nst_bst_import_channel_flags.m @@ -0,0 +1,19 @@ +function nst_bst_import_channel_flags(sFile, channel_flags_fn) +%NST_BST_EXPORT_CHANNEL_FLAGS Save channel flags from given sFile +% into given channel_flags_fn. +% TODO: also save list of channels +if ischar(sFile) + bst_chan_data = load(file_fullpath(sFile), 'ChannelFlag'); + bst_chan_flags = bst_chan_data.ChannelFlag; +end +[root, bfn, ext] = fileparts(channel_flags_fn); +switch ext + case '.mat' + s.channel_flags = bst_chan_flags; + bst_save(channel_flags_fn, s, 'v7'); + otherwise + error(sprintf('Unsupported extension for channel file: "%s"', ext)); +end +end + + diff --git a/bst_plugin/nst_bst_import_events.m b/bst_plugin/nst_bst_import_events.m new file mode 100644 index 00000000..85f1d07e --- /dev/null +++ b/bst_plugin/nst_bst_import_events.m @@ -0,0 +1,49 @@ +function sFileUpdated = nst_bst_import_events(sFile, events_fn) +% Wrapper around bst functions to import events from given file event_fn into +% given sFile (file or struct), in a backward compatible way. +% Load according to format based on extension +% +% Inputs: +% - sFile (struct or str): data to import into. +% either a brainstorm data structure or a mat file containing +% the data structure +% - events_fn (str): event file +% Event information. Suported format: .mat +% +% Outputs: +% - sFile (struct or str): same type as sFile intput +% + + + +STRUCT = 1; +FILE = 2; +if isstruct(sFile) + obj_type = STRUCT; +else + assert(ischar(sFile)); + + sFile_fn = sFile; + obj_type = FILE; + + % fileType = file_gettype( fileName ) + % See process_nst_import_csv_events for proper event importation and + % loading of sFile struct from file + + + sFile = in_fopen_bstmat(DataFile); +end + +[root, bfn, ext] = fileparts(event_fn); +switch ext + case '.mat' + sFile = import_events(sFile, [], event_fn, 'BST'); + otherwise + error(sprintf('Unsupported extension for event file: "%s"', ext)); +end + +if obj_type == FILE + bst_save(sFile_fn, sFile); +end + +end \ No newline at end of file diff --git a/test/BstEventIOTest.m b/test/BstEventIOTest.m index 6de81aaf..dc566699 100644 --- a/test/BstEventIOTest.m +++ b/test/BstEventIOTest.m @@ -22,10 +22,9 @@ function tear_down(testCase) methods(Test) - - + function test_load_save_current_version(testCase) - bst_events = get_test_events(); + bst_events = utest_get_test_bst_events(); event_fn = fullfile(testCase.tmp_dir, 'bst_events.mat'); nst_bst_save_events(bst_events, event_fn); testCase.assertTrue(exist(event_fn, 'file')==2); @@ -34,7 +33,7 @@ function test_load_save_current_version(testCase) end function test_load_save_backward_compatibility(testCase) - bst_events = get_test_events(); + bst_events = utest_get_test_bst_events(); % Dump event file and check that it's available online in the % unittest nirstorm data repository forge_current_utest_data_file(bst_events); @@ -58,39 +57,6 @@ function test_load_save_backward_compatibility(testCase) end end -function bst_event = get_test_events() -% Return a BST event struct with corner cases: -% - multiple conditions -% - some have empty notes / channels tags -% - some have non-empty notes with weird chars -% - some selected, some not -% -% Maintain this function everytime bst event structure changes -% -> to add new fields use expected default values so that backward compatibility -% is guarenteed. test_load_save_backward_compatibility should yell at -% you if not. -bst_event = db_template('event'); -bst_event(1).label = 'condition1'; -bst_event(1).color = [0., 0., 0.]; -bst_event(1).epochs = [1]; -bst_event(1).times = [10.3]; -bst_event(1).reactTimes = []; -bst_event(1).select = 1; -bst_event(1).channels = {{}}; -bst_event(1).notes = {''}; - -bst_event(2).label = 'condition2'; -bst_event(2).color = [0., 1., 0.3]; -bst_event(2).epochs = [1 1 1]; -bst_event(2).times = [1.3 7.5 12.6]; -bst_event(2).reactTimes = []; -bst_event(2).select = 0; -bst_event(2).channels = {{'S1D1', 'S2D3'}, {}, {'S121', 'S4D3'}}; -bst_event(2).notes = {sprintf('\t!:/\\,#%%*+{}[]wazaéàè?`''&$'), ... - '', '12345678910-=;,1,2||'}; - -end - function forge_current_utest_data_file(bst_events) % Save given bst events to file tagged with current bst version, if not % already available. From 7b4c336193250241db17d6afac91535671e67c15 Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Mon, 27 May 2019 15:26:30 -0400 Subject: [PATCH 05/12] Fix marking exportation & utest --- bst_plugin/nst_bst_export_events.m | 10 +- bst_plugin/nst_ppl_surface_template_V1.m | 226 ++++++++++++++++------- test/SurfTplPipelineTest.m | 167 +++++++++++++++-- test/utest_get_test_bst_events.m | 32 ++++ 4 files changed, 351 insertions(+), 84 deletions(-) create mode 100644 test/utest_get_test_bst_events.m diff --git a/bst_plugin/nst_bst_export_events.m b/bst_plugin/nst_bst_export_events.m index da79a994..d26515de 100644 --- a/bst_plugin/nst_bst_export_events.m +++ b/bst_plugin/nst_bst_export_events.m @@ -1,11 +1,15 @@ -function nst_bst_export_events(sFile, bst_events, event_fn) +function nst_bst_export_events(sFile, event_fn) % Wrapper around bst functions to load event in a backward compatible mode % Load according to format based on extension -if ischar(bst_events) - bst_evt_data = load(file_fullpath(bst_events), 'Events'); +if ischar(sFile) + bst_evt_data = load(file_fullpath(sFile), 'Events'); bst_events = bst_evt_data.Events; +else + bst_events = sFile.Events; end + + [root, bfn, ext] = fileparts(event_fn); switch ext case '.mat' diff --git a/bst_plugin/nst_ppl_surface_template_V1.m b/bst_plugin/nst_ppl_surface_template_V1.m index 4f63268c..394b926c 100644 --- a/bst_plugin/nst_ppl_surface_template_V1.m +++ b/bst_plugin/nst_ppl_surface_template_V1.m @@ -3,7 +3,7 @@ % Manage a full template- and surface-based pipeline starting from raw NIRS data % up to GLM group analysis (if enough subjects). % -% IMPORTANT: although each subject has its own optode coordinate file during importation, +% IMPORTANT: although each subject can have its own optode coordinate file during importation, % the same optode coordinates are used for all subjects. These coordinates % are the one from the first given subject. % @@ -17,8 +17,8 @@ % options = NST_PPL_SURFACE_TEMPLATE_V1('get_options'); % get default pipeline options % % % Define import options (optional): -% options.moco.export_dir = 'path/to/store/motion_events' -% options.tag_bad_channels.export_dir = 'path/to/store/bad_channels' +% options.export_dir_events = 'path/to/store/events' +% options.export_dir_channel_flags = 'path/to/store/channel_flags' % % % Import some nirs data along with event markings: % subject_names = {'subj1', 'subj2'}; @@ -36,7 +36,7 @@ % % Run the pipeline (and save user markings): % NST_PPL_SURFACE_TEMPLATE_V1('analyse', options, subject_names); % Run the full pipeline % -% % For a working example see +% % For a minimal working example see % nirstorm/script/surface_template_full_group_pipeline.m % % DEFAULT_OPTIONS = NST_PPL_SURFACE_TEMPLATE_V1('get_options') @@ -55,6 +55,30 @@ % Return: % FILES_RAW: brainstorm file pathes to imported data. % +% FILES_RAW = NST_PPL_SURFACE_TEMPLATE_V1('import', OPTIONS, NIRS_FNS, SUBJECT_NAMES) +% Import all nirs files in database and use given subjects (skip if exists). +% NIRS_FNS is a cell array of str. +% If SUBJECT_NAMES is empty or not given, then use base filename as +% subject names. If not empty, then it must be a cell array of str with the +% same length as NIRS_FNS. +% +% Used options: +% - options.import.redo +% +% Return: +% FILES_RAW: brainstorm file pathes to imported data. +% +% FILES_RAW = NST_PPL_SURFACE_TEMPLATE_V1('save_markings', OPTIONS, SUBJECT_NAMES) +% Export markings (events and bad channels) for given subjects SUBJECT_NAMES. +% +% Used options: +% - options.export_dir_events +% - options.export_dir_channel_flags +% +% Return: +% FILES_RAW: brainstorm file pathes to imported data, for which markings +% were exported. +% % NST_PPL_SURFACE_TEMPLATE_V1('analyse', OPTIONS, GROUPS | SUBJECT_NAMES) % % Apply pipeline to given group(s) of subjects. @@ -68,9 +92,9 @@ % -> head models are not recomputed for each subject. % Create a dummy subject called "full_head_model...". % Then for each subject: -% -) Export user-defined inputs [optional]: TODO -% - movement events -% - bad channels +% 0) Export user-defined markings [optional]: +% - all events (especially movement artefacts) +% - channel tags % 1) Motion correction % ASSUME: event group "movement_artefacts" exists in each FILES_RAW % and has been filled by user before calling this function. @@ -97,11 +121,12 @@ % TODO: % - handle when no contrast defined % - options documentation -% - export manual inputs % - importation of manual inputs: % - wiki page % - utest -% +% - factorize plot code +% - more comprehensive plot setup to control all options + global GlobalData; assert(ischar(action)); @@ -128,10 +153,25 @@ subject_names = cell(size(arg1)); subject_names(:) = {''}; end + if isempty(options) + options = get_options(); + end + create_dirs(options); [imported_files, redone] = import_nirs_files(arg1, subject_names, options); + import_markings(imported_files(redone==1), subject_names(redone==1), options); + % Just in case markings changed in files that were not reimported, save them: + export_markings(imported_files(redone==0), subject_names(redone==0), options); varargout{1} = imported_files; varargout{2} = redone; return; + case 'save_markings' + subject_names = arg1; + orig_cond = ['origin' get_ppl_tag()]; + files_raw = cellfun(@(s) nst_get_bst_func_files(s, orig_cond , 'Raw'), ... + subject_names, 'UniformOutput', false); + export_markings(files_raw, subject_names, options); + varargout{1} = files_raw; + return; case 'analyse' otherwise error('Unknown action: %s', action); @@ -155,6 +195,8 @@ error('"export_fig" not found. Can be installed from "https://github.com/altmany/export_fig"'); end +create_dirs(options); + force_redo = options.redo_all; protocol_info = bst_get('ProtocolInfo'); @@ -170,11 +212,6 @@ db_surface_default(0, 'Cortex', iCortex); panel_protocols('RepaintTree'); -create_dir(options.fig_dir); -create_dir(options.moco.export_dir); -create_dir(options.tag_bad_channels.export_dir); -create_dir(options.GLM_group.rois_summary.csv_export_output_dir); - % Get head model precomputed for all optode pairs % (precompute it by cloning given data if needed) file_raw1 = nst_get_bst_func_files(groups(1).subject_names{1}, ['origin' get_ppl_tag()], 'Raw'); @@ -210,6 +247,8 @@ subject_name)); %#ok end + export_markings({file_raw}, {subject_name}, options); + % Run preprocessings [sFiles_preprocessed, hb_types, redone_preprocs, preproc_folder] = preprocs(file_raw, sFile_raw_fhm, options, force_redo|fhm_redone); @@ -477,21 +516,6 @@ % Compute Scalp coupling index nst_run_bst_proc([preproc_folder 'SCI'], force_redo | options.sci.redo, 'process_nst_sci', sFile_raw); -% TODO: export motion correction tagging to external file -% sRaw = load(file_fullpath(sFile_raw)); -% sExport.Events = sRaw.Events(strcmp({sRaw.Events.label}, 'movement_artefacts')); -% export_events(sExport, [], moco_export_fn); - -% TODO: export bad channel tagging information - -% TODO: plot raw input signals -% fig_bfn = sprintf('%s_%s_signals_raw.png', SubjectName, data_tag); -% fig_fn = protect_fn_str(fullfile(options.fig_dir, fig_bfn )); -% if ~isempty(options.fig_dir) && options.make_figs && options.plot_raw_signals.do && ... -% (force_redo || options.plot_raw_signals.redo || ~exist(fig_fn, 'file')) -% plot_signals(sFile_raw, fig_fn, options); -% end - % Deglitching if options.deglitch.do redo_parent = force_redo | options.deglitch.redo; @@ -732,6 +756,22 @@ options.import.redo = 0; +options.export_dir_events = ''; % Where to export all events (mirroring), + % espcially those manually defined. + % -> will be exported everytime the pipeline + % function is run, or called with the action + % 'export_markings' + % -> will be reimported everytime reimportation + % is needed. + +options.export_dir_channel_flags = ''; % Where to export channel tags (mirroring), + % espcially those manually tagged. + % -> will be exported everytime the pipeline + % function is run, or called with the action + % 'export_markings' + % -> will be reimported everytime reimportation + % is needed. + options.head_model.surface = 'cortex_lowres'; options.sci.redo = 0; @@ -743,10 +783,7 @@ options.deglitch.agrad_std_factor = 2.5; options.moco.redo = 0; -options.moco.export_dir = ''; % Where to export motion correction events manually tagged (for backup) - % -> will be exported before running the analysis - % -> will be reimported everytime the importation - % stage is run + options.resample.redo = 0; options.resample.freq = 5; % Hz @@ -759,11 +796,6 @@ options.tag_bad_channels.redo = 0; options.tag_bad_channels.max_prop_sat_ceil = 1; % no tagging options.tag_bad_channels.max_prop_sat_floor = 1; % no tagging -options.tag_bad_channels.export_dir = ''; % Where to export bad channel taggings (backup) - % -> will be exported before - % running the analysis - % -> will be reimported everytime - % the importation stage is run options.projection.redo = 0; proj_methods = process_nst_cortical_projection('methods'); @@ -854,57 +886,113 @@ 'freq', [], ... 'baseline', []); redone_imports(ifile) = redone; - %% Manage movement event markings TODO if redone - evt_formats = bst_get('FileFilters', 'events'); - evt_format = evt_formats(strcmp('BST', evt_formats(:,3)), :); + % Create empty event group + movement_events = db_template('event'); + movement_events.label = 'movement_artefacts'; + sFile_in = bst_process('GetInputStruct', file_in); + process_nst_import_csv_events('import_events', [], sFile_in, movement_events); + end + files_in{ifile} = file_in; +end + +end + +function create_dirs(options) +create_dir(options.fig_dir); +create_dir(options.export_dir_events); +create_dir(options.export_dir_channel_flags); +create_dir(options.GLM_group.rois_summary.csv_export_output_dir); +end + + +% TODO: export motion correction tagging to external file +% sRaw = load(file_fullpath(sFile_raw)); +% sExport.Events = sRaw.Events(strcmp({sRaw.Events.label}, 'movement_artefacts')); +% export_events(sExport, [], moco_export_fn); + +% TODO: export bad channel tagging information - moco_fn = get_moco_markings_fn(subject_name, options.moco.export_dir); - if exist(moco_fn, 'file') +% TODO: plot raw input signals +% fig_bfn = sprintf('%s_%s_signals_raw.png', SubjectName, data_tag); +% fig_fn = protect_fn_str(fullfile(options.fig_dir, fig_bfn )); +% if ~isempty(options.fig_dir) && options.make_figs && options.plot_raw_signals.do && ... +% (force_redo || options.plot_raw_signals.redo || ~exist(fig_fn, 'file')) +% plot_signals(sFile_raw, fig_fn, options); +% end + + +% if exist(moco_fn, 'file') % Load event from pre-saved file % TODO: test - sFile_in = load(file_fullpath(file_in)); - [sFile_in, events] = import_events(sFile_in, [], moco_fn, evt_format); - else - % Create empty event group - movement_events = db_template('event'); - movement_events.label = 'movement_artefacts'; - sFile_in = bst_process('GetInputStruct', file_in); - process_nst_import_csv_events('import_events', [], sFile_in, movement_events); - end +% sFile_in = load(file_fullpath(file_in)); +% [sFile_in, events] = import_events(sFile_in, [], moco_fn, evt_format); +% else + +function import_markings(sFiles, subject_names, options) +io_markings(sFiles, subject_names, @nst_bst_import_events, @get_events_fn, ... + 'export_dir_events', 1, 'Loaded events', options); +io_markings(sFiles, subject_names, @nst_bst_import_channel_flags, @get_channel_flags_fn, ... + 'export_dir_channel_flags', 1, 'Loaded channel flags', options); +end + +function export_markings(sFiles, subject_names, options) +io_markings(sFiles, subject_names, @nst_bst_export_events, @get_events_fn, ... + 'export_dir_events', 0, 'Saved events', options); +io_markings(sFiles, subject_names, @nst_bst_export_channel_flags, @get_channel_flags_fn, ... + 'export_dir_channel_flags', 0, 'Saved channel flags', options); +end - bad_chans_fn = get_bad_chan_markings_fn(subject_name, options.tag_bad_channels.export_dir); - if exist(bad_chans_fn, 'file') - % TODO: load content of .mat and set channel flag - % TODO: save data -> see process_nst_tag_bad_channels +function io_markings(sFiles, subject_names, io_func, get_fn_func, dir_option_name, ... + check_exist, msg, options) +if ~isempty(options.(dir_option_name)) + assert(exist(options.(dir_option_name), 'dir')~=0); + for isubject=1:length(subject_names) + markings_fn = get_fn_func(subject_names{isubject},.... + options.(dir_option_name)); + if ~check_exist || exist(markings_fn, 'file')==2 + io_func(sFiles{isubject}, markings_fn); + write_log(sprintf('%s to: %s\n', msg, markings_fn)); end end - %% Manage bad channel markins - % TODO: update channel flags - % - files_in{ifile} = file_in; end - end -function markings_fn = get_moco_markings_fn(subject_name, export_dir) -markings_fn = ''; -if ~isempty(export_dir) - assert(exist(export_dir, 'dir')~=0); - markings_fn = fullfile(export_dir, [subject_name '_motion_events.mat']); + +% function save_markings_(sFiles, subject_names, export_func, get_fn_func, dir_option_name, msg, options) +% if ~isempty(options.(dir_option_name)) +% assert(exist(options.(dir_option_name), 'dir')~=0); +% for isubject=1:length(subject_names) +% markings_fn = get_fn_func(subject_names{isubject},.... +% options.(dir_option_name)); +% export_func(sFiles{isubject}, markings_fn); +% write_log(sprintf('%s to: %s\n', msg, markings_fn)); +% end +% end +% end + +function events_fn = get_events_fn(subject_name, export_dir) +events_fn = get_marking_fn(subject_name, export_dir, 'events'); end + +function chan_tag_fn = get_channel_flags_fn(subject_name, export_dir) +chan_tag_fn = get_marking_fn(subject_name, export_dir, 'channel_flags'); end -function markings_fn = get_bad_chan_markings_fn(subject_name, export_dir) -markings_fn = ''; +function marking_fn = get_marking_fn(subject_name, export_dir, tag) +marking_fn = ''; if ~isempty(export_dir) assert(exist(export_dir, 'dir')~=0); - markings_fn = fullfile(export_dir, [subject_name '_bad_channels.mat']); + marking_fn = fullfile(export_dir, [subject_name '_' tag '.mat']); end end %% Helper functions +function write_log(msg) +fprintf(msg); +end + function folder = create_dir(folder) % Create folder if does not exist. % Check that folder is not a subfolder of nirstorm sources (encourage good practice diff --git a/test/SurfTplPipelineTest.m b/test/SurfTplPipelineTest.m index 122f841f..a127e576 100644 --- a/test/SurfTplPipelineTest.m +++ b/test/SurfTplPipelineTest.m @@ -1,7 +1,8 @@ -classdef SurfTplPipelineImportTest < matlab.unittest.TestCase +classdef SurfTplPipelineTest < matlab.unittest.TestCase properties tmp_dir + ppl_tag end methods(TestMethodSetup) @@ -10,6 +11,15 @@ function setup(testCase) mkdir(tmpd); testCase.tmp_dir = tmpd; utest_bst_setup(); + + ProtocolName = 'nst_utest'; + %% Ensure that nst_utest protocol exists + if isempty(bst_get('Protocol', ProtocolName)) + % Create new protocol + gui_brainstorm('CreateProtocol', ProtocolName, 1, 0); %UseDefaultAnat=1, UseDefaultChannel=0. + end + + testCase.ppl_tag = '__nspst_V1'; end end @@ -22,23 +32,156 @@ function tear_down(testCase) methods(Test) - function test_import_subjects(testCase) - % TODO -% nirs_fns = simulate_nirs(); -% options = nst_ppl_surface_template_V1('setup'); -% options.moco.export_dir = fullfile(testCase.tmp_dir, 'moco'); -% sFilesRaw = nst_ppl_surface_template_V1('import', {'data1.nirs', 'data2.nirs'}, {'subj1', 'subj2'}); + % TODO: test figure output & rewriting + % TODO: check coverage +% function test_import_subjects_from_scratch(testCase) +% [nirs_fns, subject_names] = load_test_group_data(); +% sFiles = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); +% +% % Check that all origin/Raw items are there +% for isubject=1:length(subject_names) +% sFileRaw = nst_get_bst_func_files(subject_names{isubject}, ['origin' testCase.ppl_tag], 'Raw'); +% testCase.assertNotEmpty(sFileRaw); +% testCase.assertMatches(sFiles{isubject}, sFileRaw); +% % Check that event group for movement artefacts is there +% data_event = load(file_fullpath(sFileRaw), 'Events'); +% events = data_event.Events; +% testCase.assertTrue(ismember('movement_artefacts', {events.label})); +% end +% end - end +% function test_reimportation(testCase) +% % Do a first importation of all data, then remove one item and +% % check that only this one is reimported while running the +% % importation again +% [nirs_fns, subject_names] = load_test_group_data(); +% sFiles = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); +% bst_process('CallProcess', 'process_delete', sFiles{2}, [], ... +% 'target', 1); +% [sFiles, reimported] = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); +% testCase.assertTrue(all(reimported==[0 1 0])); +% sFileRaw = nst_get_bst_func_files(subject_names{2}, ['origin' testCase.ppl_tag], 'Raw'); +% testCase.assertNotEmpty(sFileRaw); +% testCase.assertMatches(sFiles{2}, sFileRaw); +% end - function test_import_events(testCase) +% + function test_save_markings(testCase) + % Do a first importation of all data, then add some events & + % tag some channels. + % Run importation again and check that user-defined markings are + % saved. + % TODO: add action 'save_markings' to explicitely save all + % events and bad channels tagging. + + % TODO: reimplement + % - utest them with MWUC from ppl function & insure backward compatibility + % - see if we still need to keep nst_bst_save_events / nst_bst_load_events + [nirs_fns, subject_names] = load_test_group_data(); + options = nst_ppl_surface_template_V1('get_options'); + event_dir = fullfile(testCase.tmp_dir, 'export_events'); + options.export_dir_events = event_dir; + chan_flags_dir = fullfile(testCase.tmp_dir, 'export_channel_flags'); + options.export_dir_channel_flags = chan_flags_dir; + + sFiles = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); + data = load(file_fullpath(sFiles{2})); + events = data.Events; + i_evt_mvt = find(strcmp('movement_artefacts', {events.label})); + testCase.assertEqual(length(i_evt_mvt), 1); + + % Add one movement artefacts + events(i_evt_mvt).times = [1.5;... + 1.87]; + events(i_evt_mvt).epochs = [1]; + events(i_evt_mvt).reactTimes = []; + events(i_evt_mvt).channels = {{}}; + events(i_evt_mvt).notes = {'participant sneezed'}; + + % Append test events + events = [events utest_get_test_bst_events()]; + + % Tag some channels as bad: + data.ChannelFlag(1) = -1; + data.ChannelFlag(end) = -1; + + % Update bst data + data.Events = events; + save(file_fullpath(sFiles{2}), '-struct', 'data'); %TODO: use bst import function instead of hacking data file + + % Save markings + sFiles_mkgs = nst_ppl_surface_template_V1('save_markings', options, subject_names); + testCase.assertMatches(sFiles_mkgs{1}, sFiles{1}); + testCase.assertMatches(sFiles_mkgs{2}, sFiles{2}); + testCase.assertMatches(sFiles_mkgs{3}, sFiles{3}); + testCase.assertEqual(length(sFiles_mkgs), length(sFiles)); + + % Check event markings + get_evt_fn = @(ii) fullfile(event_dir, [subject_names{ii} '_events.mat']); + testCase.assertTrue(exist(get_evt_fn(1), 'file')==2); + loaded_events = load(get_evt_fn(1)); + loaded_events = loaded_events.events; + testCase.assertTrue(length(loaded_events) == 1); + testCase.assertMatches(loaded_events(1).label, 'movement_artefacts'); + + event_fn = get_evt_fn(2); + testCase.assertTrue(exist(event_fn, 'file')==2); + loaded_events = load(event_fn); + loaded_events = loaded_events.events; + testCase.assertTrue(isequaln(loaded_events, events)); + + % Check bad channel markings + get_ctag_fn = @(ii) fullfile(chan_flags_dir, [subject_names{ii} '_channel_flags.mat']); + testCase.assertTrue(exist(get_ctag_fn(1), 'file')==2); + + ctag_fn = get_ctag_fn(2); + testCase.assertTrue(exist(ctag_fn, 'file')==2); + loaded_flags = load(ctag_fn); + loaded_flags = loaded_flags.channel_flags; + testCase.assertEqual(loaded_flags(1), -1); + testCase.assertEqual(loaded_flags(end), -1); + testCase.assertTrue(all(loaded_flags(2:end-1) == 1)); end +% +% function test_reimportation_with_markings(testCase) +% % Do a first importation of all data, then add some event & tag chans +% % Save event markings. Then remove item. Reimport and check +% % that user-defined event markings were restored - function test_import_bad_channels(testCase) - - end +% end +% +% function test_reimportation_with_markings_events_only(testCase) +% % Do not enable saving of bad chans tags. +% % Do a first importation of all data, then add some event & chan tags +% % Save markings. Then remove item. Reimport and check +% % that user-defined event markings were restored but not chan tags. +% end +% +% function test_reimportation_with_markings_bad_chans_only(testCase) +% % Do not enable saving of events. +% % Do a first importation of all data, then add some event & chan tags +% % Save markings. Then remove item. Reimport and check +% % that chan tags were restored but not event markings. +% end +% +% function test_default_pipeline(testCase) +% % Import data set. Run pipeline with minimal options. +% % Insure that the following outputs are produced: +% % - SCI +% % - resampled data +% % - dOD +% % - Projected data +% % - GLM 1st level with contrast effect maps +% % - GLM 2nd level with contrast t-maps +% end end end + +function [nirs_fns, subject_names] = load_test_group_data() +[nirs_fns, subject_names] = nst_io_fetch_sample_data('template_group_tapping'); +nirs_fns = nirs_fns(1:3); +subject_names = subject_names(1:3); +end \ No newline at end of file diff --git a/test/utest_get_test_bst_events.m b/test/utest_get_test_bst_events.m new file mode 100644 index 00000000..ac68167d --- /dev/null +++ b/test/utest_get_test_bst_events.m @@ -0,0 +1,32 @@ +function bst_events = utest_get_test_bst_events() +% Return a BST event struct with corner cases: +% - multiple conditions +% - some have empty notes / channels tags +% - some have non-empty notes with weird chars +% - some selected, some not +% +% Maintain this function everytime bst event structure changes +% -> to add new fields use expected default values so that backward compatibility +% is guarenteed. test_load_save_backward_compatibility should yell at +% you if not. +bst_events = db_template('event'); +bst_events(1).label = 'condition1'; +bst_events(1).color = [0., 0., 0.]; +bst_events(1).epochs = [1]; +bst_events(1).times = [10.3;10.9]; +bst_events(1).reactTimes = []; +bst_events(1).select = 1; +bst_events(1).channels = {{}}; +bst_events(1).notes = {''}; + +bst_events(2).label = 'condition2'; +bst_events(2).color = [0., 1., 0.3]; +bst_events(2).epochs = [1 1 1]; +bst_events(2).times = [1.3 7.5 12.6]; +bst_events(2).reactTimes = []; +bst_events(2).select = 0; +bst_events(2).channels = {{'S1D1', 'S2D3'}, {}, {'S121', 'S4D3'}}; +bst_events(2).notes = {sprintf('\t!:/\\,#%%*+{}[]wazaéàè?`''&$'), ... + '', '12345678910-=;,1,2||'}; + +end \ No newline at end of file From ed244e28a51f3dd1dddcaab8daa1583beb938748 Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Tue, 28 May 2019 10:26:29 -0400 Subject: [PATCH 06/12] Fix channel flag import/export. Add utest with imported data, not only raw. All utests pass for markings IO --- bst_plugin/nst_bst_export_events.m | 3 +- bst_plugin/nst_bst_import_channel_flags.m | 60 +++++- bst_plugin/nst_bst_import_events.m | 36 ++-- .../surface_template_full_group_pipeline_V1.m | 5 +- test/CSVImportTest.m | 66 +++++- test/SurfTplPipelineTest.m | 203 +++++++++++------- test/utest_import_nirs_in_bst.m | 48 +++-- 7 files changed, 306 insertions(+), 115 deletions(-) diff --git a/bst_plugin/nst_bst_export_events.m b/bst_plugin/nst_bst_export_events.m index d26515de..6bc26c09 100644 --- a/bst_plugin/nst_bst_export_events.m +++ b/bst_plugin/nst_bst_export_events.m @@ -1,6 +1,5 @@ function nst_bst_export_events(sFile, event_fn) -% Wrapper around bst functions to load event in a backward compatible mode -% Load according to format based on extension +% Wrapper around bst functions to export events. if ischar(sFile) bst_evt_data = load(file_fullpath(sFile), 'Events'); diff --git a/bst_plugin/nst_bst_import_channel_flags.m b/bst_plugin/nst_bst_import_channel_flags.m index 6f544f86..ecd8d1b2 100644 --- a/bst_plugin/nst_bst_import_channel_flags.m +++ b/bst_plugin/nst_bst_import_channel_flags.m @@ -1,19 +1,63 @@ -function nst_bst_import_channel_flags(sFile, channel_flags_fn) +function sFileUpdated = nst_bst_import_channel_flags(sFile, channel_flags_fn) %NST_BST_EXPORT_CHANNEL_FLAGS Save channel flags from given sFile % into given channel_flags_fn. -% TODO: also save list of channels -if ischar(sFile) - bst_chan_data = load(file_fullpath(sFile), 'ChannelFlag'); - bst_chan_flags = bst_chan_data.ChannelFlag; +% TODO: also save list of channels to always insure that chan flags array +% is aligned with channel info +% +% Inputs: +% - sFile (struct or str): data to import into. +% either a brainstorm data structure or a mat file containing +% the data structure +% - events_fn (str): channel flags file +% Channel flag information. Suported format: .mat +% +% Outputs: +% - sFile (struct or str): same type as sFile intput +% +% TODO: utest with struct input + +STRUCT = 1; +FILE = 2; +if isstruct(sFile) + obj_type = STRUCT; +else + assert(ischar(sFile)); + + sFile_fn = sFile; + obj_type = FILE; + + isRaw = ~isempty(strfind(sFile_fn, '_0raw')); + if isRaw + DataMat = in_bst_data(sFile_fn, 'F'); + sFile = DataMat.F; + else + DataMat = in_fopen_bstmat(sFile_fn); + sFile = DataMat; + end end + [root, bfn, ext] = fileparts(channel_flags_fn); +sFileUpdated = sFile; switch ext case '.mat' - s.channel_flags = bst_chan_flags; - bst_save(channel_flags_fn, s, 'v7'); + chan_flags_data = load(channel_flags_fn); + sFileUpdated.ChannelFlag = chan_flags_data.channel_flags; + %TODO: save channel flags properly -> should be in channel header + %file, not data file! otherwise error(sprintf('Unsupported extension for channel file: "%s"', ext)); end -end +if obj_type == FILE + % Report changes in .mat structure + DataMat.ChannelFlag = sFileUpdated.ChannelFlag; + if isRaw + DataMat.F.channelflag = chan_flags_data.channel_flags; + end + % Save file definition + bst_save(file_fullpath(sFile_fn), DataMat, 'v6', 1); + + sFileUpdated = sFile_fn; +end +end diff --git a/bst_plugin/nst_bst_import_events.m b/bst_plugin/nst_bst_import_events.m index 85f1d07e..c5f959e9 100644 --- a/bst_plugin/nst_bst_import_events.m +++ b/bst_plugin/nst_bst_import_events.m @@ -13,8 +13,7 @@ % Outputs: % - sFile (struct or str): same type as sFile intput % - - +% TODO: utest with struct input STRUCT = 1; FILE = 2; @@ -25,25 +24,36 @@ sFile_fn = sFile; obj_type = FILE; - - % fileType = file_gettype( fileName ) - % See process_nst_import_csv_events for proper event importation and - % loading of sFile struct from file - - sFile = in_fopen_bstmat(DataFile); + isRaw = ~isempty(strfind(sFile_fn, '_0raw')); + if isRaw + DataMat = in_bst_data(sFile_fn, 'F'); + sFile = DataMat.F; + else + DataMat = in_fopen_bstmat(sFile_fn); + sFile = DataMat; + end end -[root, bfn, ext] = fileparts(event_fn); +[root, bfn, ext] = fileparts(events_fn); switch ext case '.mat' - sFile = import_events(sFile, [], event_fn, 'BST'); + sFileUpdated = import_events(sFile, [], events_fn, 'BST'); otherwise error(sprintf('Unsupported extension for event file: "%s"', ext)); end -if obj_type == FILE - bst_save(sFile_fn, sFile); +if obj_type == FILE + % Report changes in .mat structure + if isRaw + DataMat.F = sFileUpdated; + else + DataMat.Events = sFileUpdated.events; + end + % Save file definition + bst_save(file_fullpath(sFile_fn), DataMat, 'v6', 1); + + sFileUpdated = sFile_fn; end -end \ No newline at end of file +end diff --git a/scripts/surface_template_full_group_pipeline_V1.m b/scripts/surface_template_full_group_pipeline_V1.m index c4cecd9a..4e1f04c1 100644 --- a/scripts/surface_template_full_group_pipeline_V1.m +++ b/scripts/surface_template_full_group_pipeline_V1.m @@ -11,7 +11,7 @@ function surface_template_full_group_pipeline_V1() % % Data are imported in brainstorm in a dedicated protocol, using specific % naming conventions specific to NST_PPL_SURFACE_TEMPLATE_V1. - % Preprocessings and processings are then run up to group-level analysis: + % Preprocessings and processings are then run, up to group-level analysis: % 1) Resampling to 5hz % 2) Bad channels detection % 3) Conversion to delta optical density @@ -76,7 +76,6 @@ function surface_template_full_group_pipeline_V1() options.GLM_1st_level.contrasts(1).label = 'motor'; options.GLM_1st_level.contrasts(1).vector = '[1 0]'; % a vector of weights, as a string -% Run the pipeline (and save user markings): +% Run the pipeline: nst_ppl_surface_template_V1('analyse', options, subject_names); % Run the full pipeline -%TODO: full reload of GUI tree at end of pipeline end \ No newline at end of file diff --git a/test/CSVImportTest.m b/test/CSVImportTest.m index 8990e96b..f034a56b 100644 --- a/test/CSVImportTest.m +++ b/test/CSVImportTest.m @@ -69,6 +69,58 @@ function test_csv_import_event_merge(testCase) assert(abs(events(i_evt).times(1, 2) - 1482333540.553) <= nirs_dt); assert(abs(events(i_evt).times(2, 2) - 1482333570.823) <= nirs_dt); end + + function test_evt_merge_imported_data(testCase) + paradigm_fn = utest_request_data({'lesca_data','lesca_task_data_block.xls'}); + paradigm_file_sel = {paradigm_fn, 'CSV'}; + nirs_fn = utest_request_data({'dummy.nirs'}); + nirs = load(nirs_fn, '-mat'); + nirs_dt = diff(nirs.t(1:2)); + + sFiles_dummy = utest_import_nirs_in_bst(nirs_fn, 1, 0); + + event_orig = db_template('event'); + event_orig.label = 'DMNirs'; + event_orig.times = [10 ; 15]; + event_orig.epochs = ones(1, size(event_orig.times, 2)); + event_orig.channels = cell(1, size(event_orig.times, 2)); + event_orig.notes = cell(1, size(event_orig.times, 2)); + process_nst_import_csv_events('import_events', [], sFiles_dummy, event_orig); + + trial_span_types = process_nst_import_csv_events('get_trial_spans_opt'); + time_units = process_nst_import_csv_events('get_time_units_opt'); + time_origin_types = process_nst_import_csv_events('get_time_origin_opt'); + + sFiles_processed = bst_process('CallProcess', ... + 'process_nst_import_csv_events', ... + sFiles_dummy, [], ... + 'evtfile', paradigm_file_sel, ... + 'delimiter', '\t', ... + 'trial_label_column', 'blockType', ... + 'trial_start_column', 'blockClockBegin', ... + 'span_type', trial_span_types.START_END, ... + 'trial_end_column', 'blockClockEnd', ... + 'time_unit', time_units.MILLISECOND, ... + 'time_origin_type', time_origin_types.OFFSET, ... + 'time_origin_offset_sec', {0.0, '', 0}, ... + 'entry_filters', 'dbCurrentSessionId=22555,blockType=DMNirs', ... + 'max_events', 1000, ... + 'confirm_importation', 0); + assert(~isempty(sFiles_processed)); + [events, DataMat] = get_events(sFiles_processed); + + % TODO: also check DataMat content aside events + assert(length(events) == 1); + assert(size(events(1).times, 1) == 2); + + i_evt = strcmp({events.label}, 'DMNirs'); + assert(events(i_evt).times(1, 1) == event_orig.times(1, 1)); + assert(events(i_evt).times(2, 1) == event_orig.times(2, 1)); + assert(abs(events(i_evt).times(1, 2) - 1482333540.553) <= nirs_dt); + assert(abs(events(i_evt).times(2, 2) - 1482333570.823) <= nirs_dt); + end + + function test_csv_import_filters(testCase) global GlobalData; @@ -497,7 +549,15 @@ function test_csv_event_import_time_unit(testCase) end end -function events = get_events(sFile) -DataMat = in_bst_data(sFile.FileName, 'F'); -events = DataMat.F.events; +function [events, DataMat, isRaw] = get_events(sFile) + +isRaw = ~isempty(strfind(sFile.FileName, '_0raw')); +if isRaw + DataMat = in_bst_data(sFile.FileName, 'F'); + events = DataMat.F.events; +else + DataMat = in_fopen_bstmat(sFile.FileName); + events = DataMat.events; +end + end diff --git a/test/SurfTplPipelineTest.m b/test/SurfTplPipelineTest.m index a127e576..fbbe4f61 100644 --- a/test/SurfTplPipelineTest.m +++ b/test/SurfTplPipelineTest.m @@ -34,49 +34,47 @@ function tear_down(testCase) % TODO: test figure output & rewriting % TODO: check coverage - -% function test_import_subjects_from_scratch(testCase) -% [nirs_fns, subject_names] = load_test_group_data(); -% sFiles = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); -% -% % Check that all origin/Raw items are there -% for isubject=1:length(subject_names) -% sFileRaw = nst_get_bst_func_files(subject_names{isubject}, ['origin' testCase.ppl_tag], 'Raw'); -% testCase.assertNotEmpty(sFileRaw); -% testCase.assertMatches(sFiles{isubject}, sFileRaw); -% % Check that event group for movement artefacts is there -% data_event = load(file_fullpath(sFileRaw), 'Events'); -% events = data_event.Events; -% testCase.assertTrue(ismember('movement_artefacts', {events.label})); -% end -% end - -% function test_reimportation(testCase) -% % Do a first importation of all data, then remove one item and -% % check that only this one is reimported while running the -% % importation again -% [nirs_fns, subject_names] = load_test_group_data(); -% sFiles = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); -% bst_process('CallProcess', 'process_delete', sFiles{2}, [], ... -% 'target', 1); -% [sFiles, reimported] = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); -% testCase.assertTrue(all(reimported==[0 1 0])); -% sFileRaw = nst_get_bst_func_files(subject_names{2}, ['origin' testCase.ppl_tag], 'Raw'); -% testCase.assertNotEmpty(sFileRaw); -% testCase.assertMatches(sFiles{2}, sFileRaw); -% end - - -% + + function test_import_subjects_from_scratch(testCase) + [nirs_fns, subject_names] = load_test_group_data(); + sFiles = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); + + % Check that all origin/Raw items are there + for isubject=1:length(subject_names) + sFileRaw = nst_get_bst_func_files(subject_names{isubject}, ['origin' testCase.ppl_tag], 'Raw'); + testCase.assertNotEmpty(sFileRaw); + testCase.assertMatches(sFiles{isubject}, sFileRaw); + % Check that event group for movement artefacts is there + data_event = load(file_fullpath(sFileRaw), 'Events'); + events = data_event.Events; + testCase.assertTrue(ismember('movement_artefacts', {events.label})); + end + end + + function test_reimportation(testCase) + % Do a first importation of all data, then remove one item and + % check that only this one is reimported while running the + % importation again + [nirs_fns, subject_names] = load_test_group_data(); + sFiles = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); + bst_process('CallProcess', 'process_delete', sFiles{2}, [], ... + 'target', 1); + [sFiles, reimported] = nst_ppl_surface_template_V1('import', [], nirs_fns, subject_names); + testCase.assertTrue(all(reimported==[0 1 0])); + sFileRaw = nst_get_bst_func_files(subject_names{2}, ['origin' testCase.ppl_tag], 'Raw'); + testCase.assertNotEmpty(sFileRaw); + testCase.assertMatches(sFiles{2}, sFileRaw); + end + function test_save_markings(testCase) % Do a first importation of all data, then add some events & % tag some channels. - % Run importation again and check that user-defined markings are + % Run importation again and check that user-defined markings are % saved. % TODO: add action 'save_markings' to explicitely save all % events and bad channels tagging. - % TODO: reimplement + % TODO: reimplement % - utest them with MWUC from ppl function & insure backward compatibility % - see if we still need to keep nst_bst_save_events / nst_bst_load_events [nirs_fns, subject_names] = load_test_group_data(); @@ -94,7 +92,7 @@ function test_save_markings(testCase) % Add one movement artefacts events(i_evt_mvt).times = [1.5;... - 1.87]; + 1.87]; events(i_evt_mvt).epochs = [1]; events(i_evt_mvt).reactTimes = []; events(i_evt_mvt).channels = {{}}; @@ -107,18 +105,18 @@ function test_save_markings(testCase) data.ChannelFlag(1) = -1; data.ChannelFlag(end) = -1; - % Update bst data + % Update bst data data.Events = events; - save(file_fullpath(sFiles{2}), '-struct', 'data'); %TODO: use bst import function instead of hacking data file + save(file_fullpath(sFiles{2}), '-struct', 'data'); %TODO: use bst import function instead of hacking data file? - % Save markings + % Save markings on drive sFiles_mkgs = nst_ppl_surface_template_V1('save_markings', options, subject_names); testCase.assertMatches(sFiles_mkgs{1}, sFiles{1}); testCase.assertMatches(sFiles_mkgs{2}, sFiles{2}); testCase.assertMatches(sFiles_mkgs{3}, sFiles{3}); testCase.assertEqual(length(sFiles_mkgs), length(sFiles)); - % Check event markings + % Check event markings that were saved on drive get_evt_fn = @(ii) fullfile(event_dir, [subject_names{ii} '_events.mat']); testCase.assertTrue(exist(get_evt_fn(1), 'file')==2); loaded_events = load(get_evt_fn(1)); @@ -144,39 +142,98 @@ function test_save_markings(testCase) testCase.assertEqual(loaded_flags(end), -1); testCase.assertTrue(all(loaded_flags(2:end-1) == 1)); end -% -% function test_reimportation_with_markings(testCase) -% % Do a first importation of all data, then add some event & tag chans -% % Save event markings. Then remove item. Reimport and check -% % that user-defined event markings were restored - -% end -% -% function test_reimportation_with_markings_events_only(testCase) -% % Do not enable saving of bad chans tags. -% % Do a first importation of all data, then add some event & chan tags -% % Save markings. Then remove item. Reimport and check -% % that user-defined event markings were restored but not chan tags. -% end -% -% function test_reimportation_with_markings_bad_chans_only(testCase) -% % Do not enable saving of events. -% % Do a first importation of all data, then add some event & chan tags -% % Save markings. Then remove item. Reimport and check -% % that chan tags were restored but not event markings. -% end -% -% function test_default_pipeline(testCase) -% % Import data set. Run pipeline with minimal options. -% % Insure that the following outputs are produced: -% % - SCI -% % - resampled data -% % - dOD -% % - Projected data -% % - GLM 1st level with contrast effect maps -% % - GLM 2nd level with contrast t-maps -% end + function test_reimportation_with_markings(testCase) + % Do a first importation of all data, then add some event & tag chans + % Save event markings. Then remove 1 item. Reimport and check + % that user-defined event markings were restored + + [nirs_fns, subject_names] = load_test_group_data(); + options = nst_ppl_surface_template_V1('get_options'); + event_dir = fullfile(testCase.tmp_dir, 'export_events'); + options.export_dir_events = event_dir; + chan_flags_dir = fullfile(testCase.tmp_dir, 'export_channel_flags'); + options.export_dir_channel_flags = chan_flags_dir; + + sFiles = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); + data = load(file_fullpath(sFiles{2})); + events = data.Events; + + % Add one movement artefacts + i_evt_mvt = find(strcmp('movement_artefacts', {events.label})); + events(i_evt_mvt).times = [1.5;... + 1.9]; + events(i_evt_mvt).epochs = [1]; + events(i_evt_mvt).reactTimes = []; + events(i_evt_mvt).channels = {{}}; + events(i_evt_mvt).notes = {'participant sneezed'}; + + % Append test events + events = [events utest_get_test_bst_events()]; + + % Tag some channels as bad: + data.ChannelFlag(1) = -1; + data.ChannelFlag(end) = -1; + + % Update bst data + data.Events = events; + bst_save(file_fullpath(sFiles{2}), data); %TODO: use bst import function instead of hacking data file? + + % Save markings + sFiles_mkgs = nst_ppl_surface_template_V1('save_markings', options, subject_names); + testCase.assertMatches(sFiles_mkgs{1}, sFiles{1}); + testCase.assertMatches(sFiles_mkgs{2}, sFiles{2}); + testCase.assertMatches(sFiles_mkgs{3}, sFiles{3}); + testCase.assertEqual(length(sFiles_mkgs), length(sFiles)); + + % Delete one item + bst_process('CallProcess', 'process_delete', sFiles{2}, [], ... + 'target', 1); + + % Import again + [files_raw, import_flags] = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); + testCase.assertTrue(all(import_flags == [0 1 0])); + + % Check event markings in DB + sFile1 = in_fopen(files_raw{1}, 'BST-DATA'); + testCase.assertTrue(length(sFile1.events) == 1); + testCase.assertMatches(sFile1.events(1).label, 'movement_artefacts'); + + sFile2 = in_fopen(files_raw{2}, 'BST-DATA'); + testCase.assertTrue(isequaln(sFile2.events, events)); + + % Check bad channel markings + testCase.assertTrue(all(sFile1.channelflag == 1)); + testCase.assertEqual(sFile2.channelflag(1), -1); + testCase.assertEqual(sFile2.channelflag(end), -1); + testCase.assertTrue(all(sFile2.channelflag(2:end-1) == 1)); + end + + % + % function test_reimportation_with_markings_events_only(testCase) + % % Do not enable saving of bad chans tags. + % % Do a first importation of all data, then add some event & chan tags + % % Save markings. Then remove item. Reimport and check + % % that user-defined event markings were restored but not chan tags. + % end + % + % function test_reimportation_with_markings_bad_chans_only(testCase) + % % Do not enable saving of events. + % % Do a first importation of all data, then add some event & chan tags + % % Save markings. Then remove item. Reimport and check + % % that chan tags were restored but not event markings. + % end + % + % function test_default_pipeline(testCase) + % % Import data set. Run pipeline with minimal options. + % % Insure that the following outputs are produced: + % % - SCI + % % - resampled data + % % - dOD + % % - Projected data + % % - GLM 1st level with contrast effect maps + % % - GLM 2nd level with contrast t-maps + % end end end diff --git a/test/utest_import_nirs_in_bst.m b/test/utest_import_nirs_in_bst.m index 593770b3..7b57d85d 100644 --- a/test/utest_import_nirs_in_bst.m +++ b/test/utest_import_nirs_in_bst.m @@ -1,9 +1,16 @@ -function sFile = utest_import_nirs_in_bst(nirs_fn, clean) +function sFile = utest_import_nirs_in_bst(nirs_fn, clean, raw_importation) + +% RAW=1 -> import as link to raw file +% RAW=0 -> import data into brainstorm DB. if nargin < 2 clean = 1; % always clean by default end +if nargin < 3 + raw_importation = 1; +end + ProtocolName = 'nst_utest'; %% Clean nst_utest protocol if needed @@ -24,21 +31,36 @@ gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); end -%% Import data as raw file +%% Resolve file format [rr, bfn, ext] = fileparts(nirs_fn); switch(ext) case '.nirs' %Homer - sFile = bst_process('CallProcess', 'process_import_data_raw', [], [], ... - 'subjectname', 'test_subject', ... - 'datafile', {nirs_fn, 'NIRS-BRS'}, ... - 'channelreplace', 1, ... - 'channelalign', 0); + format = 'NIRS-BRS'; case '.txt' %Artinis TODO: enable importation in Brainstorm - sFile = bst_process('CallProcess', 'process_import_data_raw', [], [], ... - 'subjectname', 'test_subject', ... - 'datafile', {nirs_fn, 'NIRS-ARTINIS'}, ... - 'channelreplace', 1, ... - 'channelalign', 0); + format = 'NIRS-ARTINIS'; end - % [sSubject, iSubject] = bst_get('Subject', sFile.SubjectName); + +%% Import data as raw or data file +subject_name = 'test_subject'; +if raw_importation + sFile = bst_process('CallProcess', 'process_import_data_raw', [], [], ... + 'subjectname', subject_name, ... + 'datafile', {nirs_fn, format}, ... + 'channelreplace', 1, ... + 'channelalign', 0); +else + sFile = bst_process('CallProcess', 'process_import_data_time', [], [], ... + 'subjectname', subject_name, ... + 'condition', 'test', ... + 'datafile', {nirs_fn, format}, ... + 'timewindow', [], ... + 'split', 0, ... + 'ignoreshort', 1, ... + 'channelalign', 0, ... + 'usectfcomp', 0, ... + 'usessp', 0, ... + 'freq', [], ... + 'baseline', []); +end + end From dec65c7aa1055e68771857fd9a7815b733aaef47 Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Fri, 14 Jun 2019 13:44:25 -0400 Subject: [PATCH 07/12] Fix usage of .samples in event --- bst_plugin/process_nst_motion_correction.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/bst_plugin/process_nst_motion_correction.m b/bst_plugin/process_nst_motion_correction.m index c0b61ee0..bc97df9c 100644 --- a/bst_plugin/process_nst_motion_correction.m +++ b/bst_plugin/process_nst_motion_correction.m @@ -84,6 +84,7 @@ if strcmp(events(ievt).label, event_name) event = events(ievt); ievt_mvt = ievt; + break; end end if isempty(event) @@ -95,7 +96,8 @@ if isempty(event.times) % no marked event data_corr = sDataIn.F'; else - data_corr = Compute(sDataIn.F', sDataIn.Time', event.samples'); + samples = round((event.times' - sDataIn.Time(1)) ./ diff(sDataIn.Time(1:2))) + 1; + data_corr = Compute(sDataIn.F', sDataIn.Time', samples); end if 0 nirs_data_full = in_bst(sInputs(iInput).FileName, [], 1, 0, 'no'); From 5596c24d09460705e15ab515829df6e1fd9eb99f Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Sat, 6 Jul 2019 10:31:20 -0400 Subject: [PATCH 08/12] Add import/export of manual markings --- bst_plugin/nst_ppl_surface_template_V1.m | 146 ++++++++++++++---- ...template_full_group_pipeline_V1_all_opts.m | 118 ++++++++++++++ test/SurfTplPipelineTest.m | 74 ++++++++- 3 files changed, 303 insertions(+), 35 deletions(-) create mode 100644 scripts/surface_template_full_group_pipeline_V1_all_opts.m diff --git a/bst_plugin/nst_ppl_surface_template_V1.m b/bst_plugin/nst_ppl_surface_template_V1.m index 394b926c..7a7a834e 100644 --- a/bst_plugin/nst_ppl_surface_template_V1.m +++ b/bst_plugin/nst_ppl_surface_template_V1.m @@ -11,34 +11,50 @@ % such as movement events and bad channels. This allows to safely flush all % brainstorm data while keeping markings. % +%% Overview +% % This function is intended to be called from batch scripts where the user % can add some custom steps. Here is the workflow: % +% %% Importation script +% % options = NST_PPL_SURFACE_TEMPLATE_V1('get_options'); % get default pipeline options % -% % Define import options (optional): -% options.export_dir_events = 'path/to/store/events' -% options.export_dir_channel_flags = 'path/to/store/channel_flags' +% % Define options (optional): +% options.export_dir_events = 'path/to/export/events'; +% options.export_dir_channel_flags = 'path/to/export/channel_flags'; % -% % Import some nirs data along with event markings: +% % Import some nirs data: % subject_names = {'subj1', 'subj2'}; % sFilesRaw = NST_PPL_SURFACE_TEMPLATE_V1('import', options, {'data1.nirs', 'data2.nirs'}, subject_names); % for ifile=1:length(sFilesRaw) % % Tweak sFilesRaw{ifile} here, eg import stimulation event. % end % -% % User can manually tag motion events and bad channels here +% % After the importation script, motion events and bad channels can be manually tagged. +% +% %% Markings export script +% +% NST_PPL_SURFACE_TEMPLATE_V1('save_markings', options, subject_names); +% +% %% Analysis script +% options = NST_PPL_SURFACE_TEMPLATE_V1('get_options'); +% % Define export options again (redundant so this could be factorized) +% options.export_dir_events = 'path/to/export/events'; +% options.export_dir_channel_flags = 'path/to/export/channel_flags'; % % % Customize options: % options.GLM_1st_level.contrasts(1).name = 'my_contrast1'; % options.GLM_1st_level.contrasts(1).vector = [0 1 -1 0]; % % % Run the pipeline (and save user markings): -% NST_PPL_SURFACE_TEMPLATE_V1('analyse', options, subject_names); % Run the full pipeline +% NST_PPL_SURFACE_TEMPLATE_V1('analyse', options); % Run the full pipeline % % % For a minimal working example see % nirstorm/script/surface_template_full_group_pipeline.m % +%% Setup and importation +% % DEFAULT_OPTIONS = NST_PPL_SURFACE_TEMPLATE_V1('get_options') % Return default options % @@ -79,7 +95,12 @@ % FILES_RAW: brainstorm file pathes to imported data, for which markings % were exported. % -% NST_PPL_SURFACE_TEMPLATE_V1('analyse', OPTIONS, GROUPS | SUBJECT_NAMES) +% FILES_RAW = NST_PPL_SURFACE_TEMPLATE_V1('save_markings', OPTIONS) +% Export markings (events and bad channels) for all subjects in current protocol. +% +%% Analysis +% +% NST_PPL_SURFACE_TEMPLATE_V1('analyse', OPTIONS, GROUPS | SUBJECT_NAMES) % % Apply pipeline to given group(s) of subjects. % ASSUME: all subjects in a given protocol have the same template @@ -95,19 +116,22 @@ % 0) Export user-defined markings [optional]: % - all events (especially movement artefacts) % - channel tags -% 1) Motion correction -% ASSUME: event group "movement_artefacts" exists in each FILES_RAW -% and has been filled by user before calling this function. -% Note that NST_PPL_SURFACE_TEMPLATE_V1('import') -% creates this event group if necessary. -% 2) Resampling: +% WARNING: the function does run the analyses again if markings have changed. +% Either use the redo option or delete specific analysis items +% in the brainstorm DB to force their recomputation. +% 1) Deglitching +% 2) Motion correction +% NOTE: NST_PPL_SURFACE_TEMPLATE_V1('import') creates the event group +% "movement_artefacts". The user should fill it to +% specify motion events. +% 3) Resampling: % TODO: check interpolation errors when there are spikes -% 3) Detect bad channels -% 4) Convert to delta optical density -% 5) High pass filter -% 6) Compute head model (from "full_head_model...") -% 7) Project on the cortical surface -% 8) 1st level GLM: +% 4) Detect bad channels +% 5) Convert to delta optical density +% 6) High pass filter +% 7) Compute head model (from "full_head_model...") +% 8) Project on the cortical surface +% 9) 1st level GLM: % - build design matrix from stimulation events % - OLS fit with pre-coloring % - compute contrasts @@ -118,10 +142,16 @@ % - MFX contrast t-maps % 2) Extract group-masked subject-level maps [optional] % -% TODO: +% TODO:% - importation of manual inputs: +% - find a way to invalidate analysis if markings have changed +% -> compare if markings about to be saved differ from disk version +% -> save only if they differ and store modification date (use ISO with +% time zone) +% -> store analysis date +% -> compare dates before running analysis, redo if analysis date < marking date +% -> for now, must be handled by the user. % - handle when no contrast defined % - options documentation -% - importation of manual inputs: % - wiki page % - utest % - factorize plot code @@ -165,10 +195,29 @@ varargout{2} = redone; return; case 'save_markings' - subject_names = arg1; + if nargin >= 3 + subject_names = arg1; + else + % Get all subjects in current protocol (ignore Group_analysis + % and subject holding full head model) + if isempty(GlobalData.DataBase.iProtocol) || (GlobalData.DataBase.iProtocol == 0) + error('Cannot find current protocol'); + end + sSubjects = GlobalData.DataBase.ProtocolSubjects(GlobalData.DataBase.iProtocol); + subject_names = ignore_forged_subjects({sSubjects.Subject.Name}); + end orig_cond = ['origin' get_ppl_tag()]; files_raw = cellfun(@(s) nst_get_bst_func_files(s, orig_cond , 'Raw'), ... subject_names, 'UniformOutput', false); + missing = cellfun(@isempty, files_raw); + if any(missing) + warning(sprintf('Missing origin/raw data files for subjects (will be ignored):\n%s\n', ... + strjoin(cellfun(@(s) sprintf(' - %s', s), subject_names(missing), ... + 'UniformOutput', false), '\n'))); + end + files_raw = files_raw(~missing); + subject_names = subject_names(~missing); + create_dirs(options); export_markings(files_raw, subject_names, options); varargout{1} = files_raw; return; @@ -227,7 +276,7 @@ all_sFiles_con = cell(1, nb_groups); for igroup=1:nb_groups redo_group = options.GLM_group.redo; - subject_names = groups(igroup).subject_names; + subject_names = ignore_forged_subjects(groups(igroup).subject_names); group_label = groups(igroup).label; %% Within-subject analyses @@ -392,12 +441,13 @@ [sFile_subj_zmat, redone] = nst_run_bst_proc(['Group analysis/' target_group_condition_names{igroup} '/' hb_types{ihb} ' | con ' contrasts(icon).label ' |' group_comment_tags{igroup} ' masked z-scores'], ... redone | options.GLM_group.rois_summary.redo, ... 'process_nst_glm_group_subjs_zmat', ... - all_sFiles_con{igroup}(ihb, icon, :), sFile_gp_masks{igroup, ihb, icon}); + all_sFiles_con{igroup}(ihb, icon, :), sFile_gp_masks{igroup, ihb, icon}, ... + 'keep_only_first_roi', options.GLM_group.rois_summary.keep_only_first_roi); all_sFiles_subj_zmat{ihb, icon} = sFile_subj_zmat; if redone % Set contrast name as prefix for each ROI column bst_process('CallProcess', 'process_nst_prefix_matrix', ... - sFile_subj_zmat, [], 'col_prefixes', [contrasts(icon).label '_']); + sFile_subj_zmat, [], 'col_prefixes', [protect_con_str(contrasts(icon).label) '_']); end end @@ -533,7 +583,7 @@ 'process_nst_motion_correction', sFile_deglitched, [], ... 'option_event_name', 'movement_artefacts'); -% Resample to 5Hz (save some space) +% Resample (save some space) redo_parent = redo_parent | options.resample.redo; [sFileMocoResampled, redo_parent] = nst_run_bst_proc([preproc_folder 'Motion-corrected | Resampled'], redo_parent, ... 'process_resample', sFileMoco, [], ... @@ -571,11 +621,18 @@ % Project and convert to d[HbX] redo_parent = redo_parent | options.projection.redo; proj_method = options.projection.method; -[sFilesHbProj, redo_parent] = nst_run_bst_proc({[preproc_folder 'dHbO_cortex'], [preproc_folder 'dHbR_cortex']}, redo_parent, ... +hb_types = process_nst_cortical_projection('get_hb_types', options.projection.compute_hbt); +% TODO: align output order with hb_types +if options.projection.compute_hbt + proj_outputs = {[preproc_folder 'dHbO_cortex'], [preproc_folder 'dHbR_cortex'], [preproc_folder 'dHbT_cortex']}; +else + proj_outputs = {[preproc_folder 'dHbO_cortex'], [preproc_folder 'dHbR_cortex']}; +end +[sFilesHbProj, redo_parent] = nst_run_bst_proc(proj_outputs, redo_parent, ... 'process_nst_cortical_projection', sFile_dOD_filtered, [], ... 'method', proj_method, ... + 'compute_hbt', options.projection.compute_hbt, ... 'sparse_storage', options.projection.sparse_storage); -hb_types = process_nst_cortical_projection('get_hb_types'); redone = redo_parent; end @@ -679,10 +736,14 @@ end end +function fm_subject = get_full_head_model_subject_name() + fm_subject = ['full_head_model' get_ppl_tag()]; +end + function [file_raw_fm, fm_subject, redone] = get_sFile_for_full_head_model(sfile_raw, options, force_redo) redone = 0; -fm_subject = ['full_head_model' get_ppl_tag()]; +fm_subject = get_full_head_model_subject_name(); subject_name = fileparts(sfile_raw); % Check if given subject is dummy one created to hold full head model if strcmp(subject_name, fm_subject) @@ -747,9 +808,16 @@ 'use_closest_wl', 1, 'use_all_pairs', 1, ... 'force_median_spread', 0, ... 'normalize_fluence', 1, ... + 'sensitivity_threshold_pct', 0.5, ... 'smoothing_fwhm', 0); end +function subject_names = ignore_forged_subjects(subject_names) +ignore_subjects = {get_full_head_model_subject_name(), 'Group_analysis'}; +subject_names = subject_names(~ismember(subject_names, ignore_subjects)); +end + + function options = get_options() options.redo_all = 0; @@ -784,6 +852,8 @@ options.moco.redo = 0; +% TODO: disable resampling by default +% Maybe issue warning if sampling > 5Hz or nb_samples > 10000 options.resample.redo = 0; options.resample.freq = 5; % Hz @@ -800,7 +870,7 @@ options.projection.redo = 0; proj_methods = process_nst_cortical_projection('methods'); options.projection.method = proj_methods.Sensitivity_based_interpolation; % proj_methods.MNE; - +options.projection.compute_hbt = 0; options.projection.sparse_storage = 0; options.clean_preprocessings = 0; @@ -830,6 +900,7 @@ options.GLM_group.rois_summary.do = 0; options.GLM_group.rois_summary.atlas = 'MarsAtlas'; options.GLM_group.rois_summary.matrix_col_prefix = ''; +options.GLM_group.rois_summary.keep_only_first_roi = 0; options.GLM_group.rois_summary.csv_export_output_dir = ''; mask_combinations = get_mask_combinations(); options.GLM_group.rois_summary.group_masks_combination = mask_combinations.none; % mask_combinations.intersection, mask_combinations.union @@ -1042,6 +1113,17 @@ function write_log(msg) end +function con_str = protect_con_str(con_str) +con_str = strrep(con_str, ' - ', '_minus_'); +con_str = strrep(con_str, '-', '_minus_'); +con_str = strrep(con_str, ' + ', '_plus_'); +con_str = strrep(con_str, '+', '_plus_'); +con_str = strrep(con_str, ' * ', '_x_'); +con_str = strrep(con_str, '*', '_x_'); +con_str = strrep(con_str, ':', '_'); +con_str = strrep(con_str, ' ', '_'); +end + function sfn = protect_fn_str(s) sfn = strrep(s, ' | ', '--'); sfn = strrep(s, ' : ', '--'); @@ -1051,9 +1133,11 @@ function write_log(msg) function plot_stat(sFile_ttest, fig_fn, options, show_colbar, save_colbar, sSubjectDefault) - % TODO: set colormap % TODO: set t-stat cmap boundaries +% TODO: control window size / figure size. +% TODO: contolr displayed scouts +global GlobalData; hFigSurfData = view_surface_data(sSubjectDefault.Surface(sSubjectDefault.iCortex).FileName, ... sFile_ttest, 'NIRS', 'NewFigure'); diff --git a/scripts/surface_template_full_group_pipeline_V1_all_opts.m b/scripts/surface_template_full_group_pipeline_V1_all_opts.m new file mode 100644 index 00000000..cc106d91 --- /dev/null +++ b/scripts/surface_template_full_group_pipeline_V1_all_opts.m @@ -0,0 +1,118 @@ +function surface_template_full_group_pipeline_V1_all_opts() + % Example for the template- and surface-based full pipeline, using the + % function NST_PPL_SURFACE_TEMPLATE_V1. + % + % This script downloads some sample data of 10 "dummy" subjects (27 Mb), + % as well as the Colin27_4NIRS template (19 Mb) if not available. + % For the analysis part, precomputed fluence data are also downloaded. + % Total max amount of data to download: 50 Mb, the user is asked for download + % confirmation. + % All downloaded data will be stored in .brainstorm/defaults/nirstorm + % + % Since several files will be written outside of the brainstorm database, + % a temporary folder will be used (see output_dir). + % + % Data are imported in brainstorm in a dedicated protocol, using specific + % naming conventions handled by NST_PPL_SURFACE_TEMPLATE_V1. + % Then, this script emulates user-defined movement artefacts markings as + % well as bad channel taggings. + % These markings will be saved in output_dir. Everytime NST_PPL_SURFACE_TEMPLATE_V1 + % is run, markings in the brainstorm DB are saved in the output_dir. + % When some nirs data is not available in the DB and need (re)importation, + % these markings are loaded at the same time. + % + % Preprocessings and processings are then run up to group-level analysis: + % 0) Figure outputs of raw time-series and SCI maps + % 1) Motion-correction + % 2) Resampling to 5hz + % 3) Bad channels detection + % 4) Conversion to delta optical density + % 5) High pass filter + % 6) Figure output of preprocessed time-series + % 5) Projection on the cortical surface using head model + % 6) 1st level GLM with pre-coloring + % 7) Subject-level t-maps with figure outputs + % 8) group-level GLM with MFX contrast t-maps + % + % This script illustrates a minimal fully functional analysis pipeline that can + % serve as a basis for another custom study. + % + % For a more detailed description, see the wiki page: + % https://github.com/Nirstorm/nirstorm/wiki/%5BWIP%5D-GLM-pipeline:-surface-and-template-based#example-script-with-all-options + % + % For a minimal example script, see: + % surface_template_full_group_pipeline_V1.m + % + +%% Setup brainstorm +if ~brainstorm('status') + % Start brainstorm without the GUI if not already running + brainstorm nogui +end + +%% Check Protocol +protocol_name = 'TestSurfaceTemplateGroupPipelineV1'; +if isempty(bst_get('Protocol', protocol_name)) + gui_brainstorm('CreateProtocol', protocol_name, 1, 0); % UseDefaultAnat=1, UseDefaultChannel=0 +end + +% Set template for default anatomy +nst_bst_set_template_anatomy('Colin27_4NIRS_Jan19'); + +%% Fetch data +% Get list of local nirs files for the group data. +% The function nst_io_fetch_sample_data takes care of downloading data to +% .brainstorm/defaults/nirstorm/sample_data if necessary +[nirs_fns, subject_names] = nst_io_fetch_sample_data('template_group_tapping'); + +%% Set root output dir (to store markings, figures and results) +tmp_folder = tempdir(); +if ~exist(tmp_folder, 'dir') + error('Cannot locate temporary folder'); +end +output_dir = fullfile(tmp_folder, 'nst_ppl_surface_template_V1_example'); +if ~exist(output_dir, 'dir') + mkdir(output_dir); +end + +%% Get default options +options = nst_ppl_surface_template_V1('get_options'); % get default pipeline options + +%% Set output directories for markings +options.moco.export_dir = fullfile(output_dir, 'markings_moco'); +options.tag_bad_channels.export_dir = fullfile(output_dir, 'markings_bad_channels'); + +%% Import data +sFiles = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); + +% Read stimulation events from AUX channel +for ifile=1:length(sFiles) + evt_data = load(file_fullpath(sFiles{ifile}), 'Events'); + if ~any(strcmp({evt_data.Events.label}, 'motor')) % Insure that events were not already loaded + bst_process('CallProcess', 'process_evt_read', sFiles{ifile}, [], ... + 'stimchan', 'NIRS_AUX', ... + 'trackmode', 3, ... % Value: detect the changes of channel value + 'zero', 0); + % Rename event AUX1 -> motor + bst_process('CallProcess', 'process_evt_rename', sFiles{ifile}, [], ... + 'src', 'AUX1', 'dest', 'motor'); + % Convert to extended event-> add duration of 30 sec to all motor events + bst_process('CallProcess', 'process_evt_extended', sFiles{ifile}, [], ... + 'eventname', 'motor', 'timewindow', [0, 30]); + end +end + +%% Simulate user-inputs +% Add tagging of movement artefacts + +% Add bad channel tagging + +%% Run pipeline +options.GLM_1st_level.stimulation_events = {'motor'}; +options.GLM_1st_level.contrasts(1).label = 'motor'; +options.GLM_1st_level.contrasts(1).vector = '[1 0]'; % a vector of weights, as a string + +% Run the pipeline (and save user markings): +nst_ppl_surface_template_V1('analyse', options, subject_names); % Run the full pipeline +%TODO: full reload of GUI tree at end of pipeline +end diff --git a/test/SurfTplPipelineTest.m b/test/SurfTplPipelineTest.m index fbbe4f61..a7349be5 100644 --- a/test/SurfTplPipelineTest.m +++ b/test/SurfTplPipelineTest.m @@ -32,6 +32,7 @@ function tear_down(testCase) methods(Test) + % TODO: test full default pipeline % TODO: test figure output & rewriting % TODO: check coverage @@ -71,11 +72,9 @@ function test_save_markings(testCase) % tag some channels. % Run importation again and check that user-defined markings are % saved. - % TODO: add action 'save_markings' to explicitely save all - % events and bad channels tagging. - % TODO: reimplement - % - utest them with MWUC from ppl function & insure backward compatibility + % TODO: + % - insure backward compatibility % - see if we still need to keep nst_bst_save_events / nst_bst_load_events [nirs_fns, subject_names] = load_test_group_data(); options = nst_ppl_surface_template_V1('get_options'); @@ -208,6 +207,73 @@ function test_reimportation_with_markings(testCase) testCase.assertEqual(sFile2.channelflag(end), -1); testCase.assertTrue(all(sFile2.channelflag(2:end-1) == 1)); end + + function test_reimportation_with_markings_all_subjects(testCase) + % Do a first importation of all data, then add some event & tag chans + % Save event markings. Then remove 1 item. Reimport and check + % that user-defined event markings were restored + + [nirs_fns, subject_names] = load_test_group_data(); + options = nst_ppl_surface_template_V1('get_options'); + event_dir = fullfile(testCase.tmp_dir, 'export_events'); + options.export_dir_events = event_dir; + chan_flags_dir = fullfile(testCase.tmp_dir, 'export_channel_flags'); + options.export_dir_channel_flags = chan_flags_dir; + + sFiles = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); + data = load(file_fullpath(sFiles{2})); + events = data.Events; + + % Add one movement artefacts + i_evt_mvt = find(strcmp('movement_artefacts', {events.label})); + events(i_evt_mvt).times = [1.5;... + 1.9]; + events(i_evt_mvt).epochs = [1]; + events(i_evt_mvt).reactTimes = []; + events(i_evt_mvt).channels = {{}}; + events(i_evt_mvt).notes = {'participant sneezed'}; + + % Append test events + events = [events utest_get_test_bst_events()]; + + % Tag some channels as bad: + data.ChannelFlag(1) = -1; + data.ChannelFlag(end) = -1; + + % Update bst data + data.Events = events; + bst_save(file_fullpath(sFiles{2}), data); %TODO: use bst import function instead of hacking data file? + + % Save markings + sFiles_mkgs = nst_ppl_surface_template_V1('save_markings', options); + testCase.assertMatches(sFiles_mkgs{1}, sFiles{1}); + testCase.assertMatches(sFiles_mkgs{2}, sFiles{2}); + testCase.assertMatches(sFiles_mkgs{3}, sFiles{3}); + testCase.assertEqual(length(sFiles_mkgs), length(sFiles)); + + % Delete one item + bst_process('CallProcess', 'process_delete', sFiles{2}, [], ... + 'target', 1); + + % Import again + [files_raw, import_flags] = nst_ppl_surface_template_V1('import', options, nirs_fns, subject_names); + testCase.assertTrue(all(import_flags == [0 1 0])); + + % Check event markings in DB + sFile1 = in_fopen(files_raw{1}, 'BST-DATA'); + testCase.assertTrue(length(sFile1.events) == 1); + testCase.assertMatches(sFile1.events(1).label, 'movement_artefacts'); + + sFile2 = in_fopen(files_raw{2}, 'BST-DATA'); + testCase.assertTrue(isequaln(sFile2.events, events)); + + % Check bad channel markings + testCase.assertTrue(all(sFile1.channelflag == 1)); + testCase.assertEqual(sFile2.channelflag(1), -1); + testCase.assertEqual(sFile2.channelflag(end), -1); + testCase.assertTrue(all(sFile2.channelflag(2:end-1) == 1)); + end + % % function test_reimportation_with_markings_events_only(testCase) From 5f5a73f1a1df7d62685941d86d8be30d011d4f1d Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Sat, 6 Jul 2019 10:32:06 -0400 Subject: [PATCH 09/12] Minor fix when some outputs already existing --- bst_plugin/nst_run_bst_proc.m | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/bst_plugin/nst_run_bst_proc.m b/bst_plugin/nst_run_bst_proc.m index 9c60f92f..916fbd28 100644 --- a/bst_plugin/nst_run_bst_proc.m +++ b/bst_plugin/nst_run_bst_proc.m @@ -158,7 +158,7 @@ sStudy = bst_get('Study', iStudy); end outputs(i_item).sStudy = sStudy; - outputs(1).iStudy = iStudy; + outputs(i_item).iStudy = iStudy; % Check if output exists in specified condition [selected_files, file_type] = nst_get_bst_func_files(subject_name, outputs(i_item).condition, outputs(i_item).comment); @@ -173,11 +173,11 @@ end if ~isempty(file_type) sFilesOut_types{i_item} = file_type; + sFilesOut{i_item} = selected_files; else sFilesOut_types{i_item} = ''; + sFilesOut{i_item} = ''; end - sFilesOut{i_item} = selected_files; - else outputs(i_item).sStudy = []; outputs(1).iStudy = []; @@ -188,9 +188,9 @@ end if ~isempty(duplicates) - bst_error(sprintf('Cannot safely manage unique outputs. Found duplicate items: %s', strjoin(duplicates, ', '))); - sFilesOut = {}; - return; + error(sprintf('Cannot safely manage unique outputs. Found duplicate items: %s', strjoin(duplicates, ', '))); +% sFilesOut = {}; +% return; end existing = cellfun(@(s) ~isempty(s), sFilesOut); @@ -204,11 +204,11 @@ prev_iHeadModel = strcmp({outputs(1).sStudy.HeadModel.FileName}, sFilesOut{1}); outputs(1).sStudy = delete_head_model(outputs(1).sStudy, outputs(1).iStudy, prev_iHeadModel); else - bst_process('CallProcess', 'process_delete', sFilesOut, [], ... + bst_process('CallProcess', 'process_delete', sFilesOut(existing), [], ... 'target', 1); end bst_report('Info', ProcessName, sFiles, ... - sprintf('Force redo - removed previous result(s): %s', strjoin(sFilesOut, ', ')) ); + sprintf('Force redo - removed previous result(s): %s', strjoin(sFilesOut(existing), ', ')) ); end % Special case for head model which is not returned in sFilesOut From 56b0d37795bb8505700c6ca25933bf86f62a84b4 Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Sat, 6 Jul 2019 10:32:40 -0400 Subject: [PATCH 10/12] Add computation of HbT --- bst_plugin/process_nst_cortical_projection.m | 27 +++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/bst_plugin/process_nst_cortical_projection.m b/bst_plugin/process_nst_cortical_projection.m index aa3221f5..5bea3126 100644 --- a/bst_plugin/process_nst_cortical_projection.m +++ b/bst_plugin/process_nst_cortical_projection.m @@ -44,6 +44,10 @@ choices = methods(); sProcess.options.method.Value = {choices.Sensitivity_based_interpolation, ... fieldnames(choices)}; + +sProcess.options.compute_hbt.Comment = 'Compute HbT'; +sProcess.options.compute_hbt.Type = 'checkbox'; +sProcess.options.compute_hbt.Value = 0; sProcess.options.sparse_storage.Comment = 'Sparse storage'; sProcess.options.sparse_storage.Type = 'checkbox'; @@ -125,7 +129,7 @@ end nb_vertices = size(sensitivity_surf, 3); -hb_types = get_hb_types(); +hb_types = get_hb_types(sProcess.options.compute_hbt.Value); pinv_ext = pinv(ext_coeffs); pdata_hbo = zeros(nb_vertices, nb_samples); @@ -149,6 +153,16 @@ extra); OutputFiles{end+1} = ResultFile; +if sProcess.options.compute_hbt.Value + pdata_hbt = pdata_hbo + pdata_hbr; + [sStudy, ResultFile] = nst_bst_add_surf_data(pdata_hbt, sDataIn.Time, ... + head_model, 'hbt_proj', sprintf('HbT cortex (%s)', mtag), ... + sInputs, sStudy, 'Projected HbT signals', [], sProcess.options.sparse_storage.Value, ... + extra); + OutputFiles{end+1} = ResultFile; + clear pdata_hbt; +end + clear pdata_hbo pdata_hbr; if sProcess.options.save_mixing_mat.Value @@ -185,9 +199,16 @@ end -function hb_types = get_hb_types() +function hb_types = get_hb_types(enable_hbt) +if nargin < 1 + enable_hbt = 0; +end % Return Hb type order as used for the cols of the matrix of exctinction coeffs -hb_types = {'HbO', 'HbR'}; +if enable_hbt + hb_types = {'HbO', 'HbR', 'HbT'}; +else + hb_types = {'HbO', 'HbR'}; +end end function save_chan_data(chan_data, sDataIn, iStudy, sStudy, comment, tag) From 9e43193cb0efe8fa2b9b5e4b7a1e4d7e3a867dee Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Sat, 6 Jul 2019 10:34:12 -0400 Subject: [PATCH 11/12] Add option to better handle scenario of single channel recording --- bst_plugin/process_nst_glm_group_subjs_zmat.m | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/bst_plugin/process_nst_glm_group_subjs_zmat.m b/bst_plugin/process_nst_glm_group_subjs_zmat.m index 81994bd3..159e540d 100644 --- a/bst_plugin/process_nst_glm_group_subjs_zmat.m +++ b/bst_plugin/process_nst_glm_group_subjs_zmat.m @@ -37,6 +37,13 @@ sProcess.InputTypes = {'results'}; sProcess.OutputTypes = {'matrix'}; + % Useful for single-channel experiment where signal is duplicated while + % projecting -> discard regional aspects + sProcess.options.keep_only_first_roi.Comment = 'Keep only first ROI'; + sProcess.options.keep_only_first_roi.Type = 'checkbox'; + sProcess.options.keep_only_first_roi.Hidden = 1; + sProcess.options.keep_only_first_roi.Value = 0; + sProcess.nInputs = 2; sProcess.nMinFiles = 1; sProcess.isPaired = 0; @@ -80,10 +87,15 @@ roi_indexes = roi_indexes(2:end); end end + + if sProcess.options.keep_only_first_roi.Value + roi_indexes = roi_indexes(1); + end + if ~isempty(mask_data.from_atlas) roi_names = {mask_data.from_atlas.Scouts(roi_indexes).Label}; else - roi_names = arrayfun(@(n) num2str(n), roi_indexes, 'UniformOutput', 0); + roi_names = arrayfun(@(n) ['roi_' num2str(n)], roi_indexes, 'UniformOutput', 0); end % Loop over full list of ROIs within FOV and fill matrix with mean of From 2d34b504e2f7774cbf4278fe3f8163fe20a6041e Mon Sep 17 00:00:00 2001 From: tvvincent <20100thomas@gmail.com> Date: Sat, 6 Jul 2019 10:34:58 -0400 Subject: [PATCH 12/12] Fix negative signals produced while doing moco --- bst_plugin/process_nst_motion_correction.m | 23 ++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/bst_plugin/process_nst_motion_correction.m b/bst_plugin/process_nst_motion_correction.m index bc97df9c..9910faa6 100644 --- a/bst_plugin/process_nst_motion_correction.m +++ b/bst_plugin/process_nst_motion_correction.m @@ -92,12 +92,29 @@ OutputFiles = {}; return end - + if isempty(event.times) % no marked event data_corr = sDataIn.F'; else + % Process only NIRS channels + channels = in_bst_channel(sInputs(iInput).ChannelFile); + nirs_ichans = channel_find(channels.Channel, 'NIRS'); + data_nirs = sDataIn.F(nirs_ichans, :)'; + prev_negs = any(data_nirs <= 0, 1); + % Keep track of channels that contained neg values samples = round((event.times' - sDataIn.Time(1)) ./ diff(sDataIn.Time(1:2))) + 1; - data_corr = Compute(sDataIn.F', sDataIn.Time', samples); + data_corr = Compute(data_nirs, sDataIn.Time', samples); + % Fix negative values created by moco + new_negs = any(data_corr <= 0, 1) & ~prev_negs; + if any(new_negs) + warning('Motion correction introduced negative values. Will be fixed by global offset'); + % TODO: use pair-specific offset see /home/tom/Projects/Research/Software/mfipcode/sandbox/nirstorm/mfip_correct_negChan.m + offset = 2*abs(min(data_corr(:))); + data_corr = data_corr + offset; + end + data_corr_full = sDataIn.F'; + data_corr_full(:, nirs_ichans) = data_corr; + data_corr = data_corr_full; end if 0 nirs_data_full = in_bst(sInputs(iInput).FileName, [], 1, 0, 'no'); @@ -114,6 +131,8 @@ % OutputFiles = {sInputs(iInput).FileName}; end + + % Save time-series data sDataOut = db_template('data'); sDataOut.F = data_corr';