From 1c268bb55f5cfdefae2aede7d215328ff8497330 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 29 Apr 2026 06:50:33 -0700 Subject: [PATCH 01/55] Add `subtightplot` external module - Adds `subtightplot.m`, A wrapper function for `subplot`. Adds the ability to define the gap between neighbouring subplots. - `subplot` prior to R2019b lacked this functionality, and the gap between subplots can reach 40% of figure area, which is pretty lavish. - In 2019b, MATLAB introduced `tiledLayout` and `nexttile` that can achieve the same but for backward compatibility it would be good to just stick to this external module. --- external/subtightplot/license.txt | 26 +++++++++++ external/subtightplot/subtightplot.m | 67 ++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 external/subtightplot/license.txt create mode 100644 external/subtightplot/subtightplot.m diff --git a/external/subtightplot/license.txt b/external/subtightplot/license.txt new file mode 100644 index 0000000000..795a7dad01 --- /dev/null +++ b/external/subtightplot/license.txt @@ -0,0 +1,26 @@ +Copyright (c) 2012, Felipe G. Nievinski +Copyright (c) 2010, Pekka Kumpulainen +Copyright (c) 2011, Nikolay S. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/external/subtightplot/subtightplot.m b/external/subtightplot/subtightplot.m new file mode 100644 index 0000000000..ee01794961 --- /dev/null +++ b/external/subtightplot/subtightplot.m @@ -0,0 +1,67 @@ +function h=subtightplot(m,n,p,gap,marg_h,marg_w,varargin) +%function h=subtightplot(m,n,p,gap,marg_h,marg_w,varargin) +% +% Functional purpose: A wrapper function for Matlab function subplot. Adds the ability to define the gap between +% neighbouring subplots. Unfotrtunately Matlab subplot function lacks this functionality, and the gap between +% subplots can reach 40% of figure area, which is pretty lavish. +% +% Input arguments (defaults exist): +% gap- two elements vector [vertical,horizontal] defining the gap between neighbouring axes. Default value +% is 0.01. Note this vale will cause titles legends and labels to collide with the subplots, while presenting +% relatively large axis. +% marg_h margins in height in normalized units (0...1) +% or [lower uppper] for different lower and upper margins +% marg_w margins in width in normalized units (0...1) +% or [left right] for different left and right margins +% +% Output arguments: same as subplot- none, or axes handle according to function call. +% +% Issues & Comments: Note that if additional elements are used in order to be passed to subplot, gap parameter must +% be defined. For default gap value use empty element- []. +% +% Usage example: h=subtightplot((2,3,1:2,[0.5,0.2]) + +if (nargin<4) || isempty(gap), gap=0.01; end +if (nargin<5) || isempty(marg_h), marg_h=0.05; end +if (nargin<5) || isempty(marg_w), marg_w=marg_h; end +if isscalar(gap), gap(2)=gap; end +if isscalar(marg_h), marg_h(2)=marg_h; end +if isscalar(marg_w), marg_w(2)=marg_w; end +gap_vert = gap(1); +gap_horz = gap(2); +marg_lower = marg_h(1); +marg_upper = marg_h(2); +marg_left = marg_w(1); +marg_right = marg_w(2); + +%note n and m are switched as Matlab indexing is column-wise, while subplot indexing is row-wise :( +[subplot_col,subplot_row]=ind2sub([n,m],p); + +% note subplot suppors vector p inputs- so a merged subplot of higher dimentions will be created +subplot_cols=1+max(subplot_col)-min(subplot_col); % number of column elements in merged subplot +subplot_rows=1+max(subplot_row)-min(subplot_row); % number of row elements in merged subplot + +% single subplot dimensions: +%height=(1-(m+1)*gap_vert)/m; +%axh = (1-sum(marg_h)-(Nh-1)*gap(1))/Nh; +height=(1-(marg_lower+marg_upper)-(m-1)*gap_vert)/m; +%width =(1-(n+1)*gap_horz)/n; +%axw = (1-sum(marg_w)-(Nw-1)*gap(2))/Nw; +width =(1-(marg_left+marg_right)-(n-1)*gap_horz)/n; + +% merged subplot dimensions: +merged_height=subplot_rows*( height+gap_vert )- gap_vert; +merged_width= subplot_cols*( width +gap_horz )- gap_horz; + +% merged subplot position: +merged_bottom=(m-max(subplot_row))*(height+gap_vert) +marg_lower; +merged_left=(min(subplot_col)-1)*(width+gap_horz) +marg_left; +pos_vec=[merged_left merged_bottom merged_width merged_height]; + +% h_subplot=subplot(m,n,p,varargin{:},'Position',pos_vec); +% Above line doesn't work as subplot tends to ignore 'position' when same mnp is utilized +h=subplot('Position',pos_vec,varargin{:}); + +if (nargout < 1), clear h; end + +end From 81f6f3cfe83bb49b84abaf0256a89dfd580927a2 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 29 Apr 2026 06:53:08 -0700 Subject: [PATCH 02/55] Add `ChanTable` to `export_channel_atlas` output Required for plotting Fastgraph as it uses the anatomical labels based on the atlas chosen by the user --- toolbox/io/export_channel_atlas.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/toolbox/io/export_channel_atlas.m b/toolbox/io/export_channel_atlas.m index 1b6ea178df..8dac6eb169 100644 --- a/toolbox/io/export_channel_atlas.m +++ b/toolbox/io/export_channel_atlas.m @@ -1,4 +1,4 @@ -function TsvFile = export_channel_atlas(ChannelFile, Modality, TsvFile, Radius, isProba, isInteractive) +function [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, Modality, TsvFile, Radius, isProba, isInteractive) % EXPORT_CHANNEL_ATLAS: Compute anatomical labels for SEEG/ECOG contacts from volume and surface parcellations % % USAGE: TsvFile = export_channel_atlas(ChannelFile, Modality='ECOG+SEEG', TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1) @@ -78,7 +78,7 @@ % ===== SELECT OUTPUT FILE ===== -if isempty(TsvFile) +if isempty(TsvFile) && isInteractive % Get default directories and formats LastUsedDirs = bst_get('LastUsedDirs'); % Default output filename From 85ae038f2fd9c481e35ff85e609dbecdb14e1738 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 1 May 2026 01:52:58 -0700 Subject: [PATCH 03/55] Add `process_fastgraph` for plotting Fastgraph --- toolbox/process/functions/process_fastgraph.m | 846 ++++++++++++++++++ 1 file changed, 846 insertions(+) create mode 100644 toolbox/process/functions/process_fastgraph.m diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m new file mode 100644 index 0000000000..98901cb2f6 --- /dev/null +++ b/toolbox/process/functions/process_fastgraph.m @@ -0,0 +1,846 @@ +function varargout = process_fastgraph( varargin ) +% PROCESS_FASTGRAPH: Plot fastgraph for one or more SEEG recordings. +% For each stimulation pair, channels are split by hemisphere, sorted +% by a user-selected metric, filtered by atlas region or scout label, and +% plotted as stacked area plots +% +% USAGE: +% OutputFiles = process_fastgraph('Run', sProcess, sInputs) + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c) University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% 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: Kenneth N. Taylor, 2020 +% John C. Mosher, 2020 +% Chinmay Chinara, 2026 + +eval(macro_method); +end + +%% ===== GET DESCRIPTION ===== +function sProcess = GetDescription() %#ok +% Describe the process and its UI options +sProcess.Comment = 'Plot Fastgraphs'; +sProcess.Category = 'Custom'; +sProcess.SubGroup = 'Stimulation'; +sProcess.Index = 1100; +% Definition of the input accepted by this process +sProcess.InputTypes = {'data'}; +sProcess.OutputTypes = {'data'}; +sProcess.nInputs = 1; +sProcess.nMinFiles = 1; +% Atlas to use for plotting Fastgraph +sProcess.options.atlas.Comment = 'Atlas to plot: '; +sProcess.options.atlas.Type = 'atlas'; +sProcess.options.atlas.Value = []; +% Color Fastgraph by region or by label +sProcess.options.label2.Comment = 'Color Fastgraph by region or by label ?'; +sProcess.options.label2.Type = 'label'; +sProcess.options.colorscheme.Comment = {'Region', 'Label'; 'Region', 'Label'}; +sProcess.options.colorscheme.Type = 'radio_label'; +sProcess.options.colorscheme.Value = 'Region'; +% Select regions to include +sProcess.options.label3.Comment = 'Select region(s) to include:'; +sProcess.options.label3.Type = 'label'; +sProcess.options.regionprefrontal.Comment = '1: Prefrontal'; +sProcess.options.regionprefrontal.Type = 'checkbox'; +sProcess.options.regionprefrontal.Value = 1; +sProcess.options.regionfrontal.Comment = '2: Frontal'; +sProcess.options.regionfrontal.Type = 'checkbox'; +sProcess.options.regionfrontal.Value = 1; +sProcess.options.regioncentral.Comment = '3: Central'; +sProcess.options.regioncentral.Type = 'checkbox'; +sProcess.options.regioncentral.Value = 1; +sProcess.options.regionparietal.Comment = '4: Parietal'; +sProcess.options.regionparietal.Type = 'checkbox'; +sProcess.options.regionparietal.Value = 1; +sProcess.options.regiontemporal.Comment = '5: Temporal'; +sProcess.options.regiontemporal.Type = 'checkbox'; +sProcess.options.regiontemporal.Value = 1; +sProcess.options.regionoccipital.Comment = '6: Occipital'; +sProcess.options.regionoccipital.Type = 'checkbox'; +sProcess.options.regionoccipital.Value = 1; +sProcess.options.regionlimbic.Comment = '7: Limbic'; +sProcess.options.regionlimbic.Type = 'checkbox'; +sProcess.options.regionlimbic.Value = 1; +% Atlas scout labels to plot +sProcess.options.label4.Comment = 'For multiple labels: separate them with commas'; +sProcess.options.label4.Type = 'label'; +sProcess.options.atlasscoutlabels.Comment = 'Atlas scout labels to plot: '; +sProcess.options.atlasscoutlabels.Type = 'text'; +sProcess.options.atlasscoutlabels.Value = ''; +% Add separator +sProcess.options.separator1.Type = 'separator'; +% Method for sorting the data +sProcess.options.label5.Comment = 'Select method to sort the data:'; +sProcess.options.label5.Type = 'label'; +sProcess.options.sortmethod.Comment = {'Root Mean Square', 'Max Absolute'}; +sProcess.options.sortmethod.Type = 'radio'; +sProcess.options.sortmethod.Value = 1; +% Sort window +sProcess.options.label6.Comment = 'Choose range to sort over:'; +sProcess.options.label6.Type = 'label'; +sProcess.options.label7.Comment = ['' ... + 'Early latency:    0-60 ms
' ... + 'Middle latency: 60-250 ms
' ... + 'Late latency:     250-600 ms']; +sProcess.options.label7.Type = 'label'; +sProcess.options.sortwindow.Comment = 'Sort range: '; +sProcess.options.sortwindow.Type = 'timewindow'; +sProcess.options.sortwindow.Value = []; +% Add separator +sProcess.options.separator2.Type = 'separator'; +% Plot window +sProcess.options.plotwindow.Comment = 'Plot range: '; +sProcess.options.plotwindow.Type = 'timewindow'; +sProcess.options.plotwindow.Value = []; +% Edge transparency of plot +sProcess.options.edgealpha.Comment = 'Edge transparency of plot: '; +sProcess.options.edgealpha.Type = 'value'; +sProcess.options.edgealpha.Value = {0.05,' ', 2}; +% Exclude contacts within a certain distance from the stimulation sites +sProcess.options.label8.Comment = ['' ... + 'Exclude analysis of contacts within this distance from the stimulation site']; +sProcess.options.label8.Type = 'label'; +sProcess.options.excluderadius.Comment = 'Exclusion zone radius: '; +sProcess.options.excluderadius.Type = 'value'; +sProcess.options.excluderadius.Value = {20,'mm', 0}; +end + +%% ===== FORMAT COMMENT ===== +function Comment = FormatComment(sProcess) %#ok + Comment = sProcess.Comment; +end + +%% ===== GET OPTIONS ===== +function OPTIONS = GetOptions(sProcess) + OPTIONS = struct(); + % Atlas to use for plotting Fastgraph + OPTIONS.Atlas = sProcess.options.atlas.Value; + % Color figure by region or by label + OPTIONS.ColorScheme = sProcess.options.colorscheme.Value; + % Select regions to include + OPTIONS.Region = logical([sProcess.options.regionprefrontal.Value + sProcess.options.regionfrontal.Value + sProcess.options.regioncentral.Value + sProcess.options.regionparietal.Value + sProcess.options.regiontemporal.Value + sProcess.options.regionoccipital.Value + sProcess.options.regionlimbic.Value]); + % Atlas scout labels to plot + OPTIONS.AtlasScoutLabels = strtrim(strsplit(sProcess.options.atlasscoutlabels.Value,',')); + % Method for sorting the data + OPTIONS.SortMethod = sProcess.options.sortmethod.Value; + % Sort window + if isfield(sProcess.options, 'sortwindow') && isfield(sProcess.options.sortwindow, 'Value') && iscell(sProcess.options.sortwindow.Value) && ~isempty(sProcess.options.sortwindow.Value) + OPTIONS.SortWindow = round((sProcess.options.sortwindow.Value{1} * 1000)) + 101; + else + OPTIONS.SortWindow = []; + end + % Plot window + if isfield(sProcess.options, 'plotwindow') && isfield(sProcess.options.plotwindow, 'Value') && iscell(sProcess.options.plotwindow.Value) && ~isempty(sProcess.options.plotwindow.Value) + OPTIONS.PlotWindow = round((sProcess.options.plotwindow.Value{1} * 1000)); + else + OPTIONS.PlotWindow = []; + end + % Edge transparency of plot + OPTIONS.EdgeAlpha = sProcess.options.edgealpha.Value{1}; + % Exclude contacts within a certain distance of stimulation sites + OPTIONS.ExcludeRadius = sProcess.options.excluderadius.Value{1}; +end + +%% ===== RUN ===== +function OutputFiles = Run(sProcess, sInputs) %#ok + % Initialize output + OutputFiles = {}; + + % Get options + OPTIONS = GetOptions(sProcess); + + % Early exit if no region is selected + if ~any(OPTIONS.Region) + bst_report('Error', sProcess, [], 'No region selected. Select at least one region to run the analysis.'); + return; + end + + % Get subject + sSubject = bst_get('Subject', sInputs(1).SubjectName); + CortexFile = sSubject.Surface(sSubject.iCortex).FileName; + sCortex = bst_memory('LoadSurface', CortexFile); + + % Get the last used atlas if atlas not selected + if isempty(OPTIONS.Atlas) + OPTIONS.Atlas = sCortex.Atlas(sCortex.iAtlas).Name; + end + % Find the atlas selected by the user + iAtlas = find(strcmpi({sCortex.Atlas.Name}, OPTIONS.Atlas), 1); + + % Early exit if any entered atlas scout label does not exist + allAtlasScoutLabels = {sCortex.Atlas(iAtlas).Scouts.Label}; + enteredLabels = OPTIONS.AtlasScoutLabels(~cellfun(@isempty, OPTIONS.AtlasScoutLabels)); + if ~isempty(enteredLabels) && ~all(ismember(enteredLabels, allAtlasScoutLabels)) + bst_report('Error', sProcess, [], 'One or more scout labels entered are not present in the selected atlas'); + return; + end + + % Load the channel file + ChannelFile = file_fullpath(sInputs(1).ChannelFile); + ChannelMat = load(ChannelFile); + % Get indices of SEEG channels + iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); + % Get the midpoint location of each stimulation pair from channel + stimLocs = GetStimLocs(sInputs, ChannelMat); + % Sort fastgraphs by stimulation-site location for LAPRAP style display + sSortedFastgraphLocIdxs = SortLAPRAP(stimLocs); + + % Load SEEG recordings after applying Fastgraph sorting + [seegData, excludedContacts] = GetSeegData(sInputs, sSortedFastgraphLocIdxs, stimLocs, ChannelMat, OPTIONS); + % Split SEEG contacts into left and right hemisphere groups + sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat); + % Compute anatomical labels for the contacts from volume/surface parcellations + [~, chanTableWithAtlas] = export_channel_atlas(ChannelFile, 'SEEG', [], 10, 0, 0); + % Locate atlas related columns from channel table above + hit = cellfun(@(x) ischar(x) && (~isempty(strfind(OPTIONS.Atlas, x)) || ~isempty(strfind(x, OPTIONS.Atlas))), chanTableWithAtlas(1,:)); + % Columns whose header matches the atlas name + cols = find(any(hit, 1)); + % Extract SEEG channel names and their atlas scout labels + chanNamesSeeg = chanTableWithAtlas(2:end, 1); + atlasScoutLabelsSeeg = chanTableWithAtlas(2:end, cols); + + % Create figure for Fastgraph + figure; + % Maximize figure + set(gcf, 'Position', get(0,'Screensize')); + % Shared y-axis limits across Fastgraph subplots + commonAxisLimits = []; + % Reserve one extra subplot for the legend + nSubplots = length(sInputs)+1; + % Define the plot parameters + % Subplot grid dimensions + nRows = floor(sqrt(nSubplots/1.5)); + nCols = ceil(nSubplots/floor(sqrt(nSubplots/1.5))); + % Subplot spacing and margins + gap = [0.075 0.0175]; + horzMargin = 0.03; + vertMargin = 0.015; + % Generate one fastgraph per selected input + bst_progress('start', 'Process', 'Plotting Fastgraphs...', 0, 100); + for iSubplot = 1:nSubplots-1 + % Show progress + progressPrc = round(100 .* iSubplot ./ (nSubplots-1)); + bst_progress('set', progressPrc); + % Data to be plotted for the current subplot + subplotData = struct(); + Fout = seegData{iSubplot}.F(iSeeg, :); + % Keep only left-hemisphere channels if present + if any(sContactGroupLocIdxs.Left) + subplotData.leftData = Fout(sContactGroupLocIdxs.Left,:); + end + % Keep only left-hemisphere channels if present + if any(sContactGroupLocIdxs.Right) + subplotData.rightData = Fout(sContactGroupLocIdxs.Right,:); + end + % Sort channels within each hemisphere using the selected metric and time window + sSubplotDataSorted = ApplyDataSorting(subplotData, seegData, OPTIONS); + % Create the subplot with custom spacing + subtightplot(nRows, nCols, iSubplot, gap, horzMargin, vertMargin); + % Plot the Fastgraph for the current stimulation pair + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInputs, stimLocs, iSubplot, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); + % Tighten axes to the plotted data and store the current axis handle + axis tight + axisLimits = axis; + axSubplots(iSubplot) = gca; + % Update the shared y-axis limits so all Fastgraph subplots can + % use the same vertical range for visual comparison + if iSubplot == 1 + commonAxisLimits = axisLimits; + else + commonAxisLimits(3) = min(commonAxisLimits(3), axisLimits(3)); + commonAxisLimits(4) = max(commonAxisLimits(4), axisLimits(4)); + end + % Apply edge transparency to the subplot + if exist('hLeftAreaPLot','var') + set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); + end + if exist('hRightAreaPLot','var') + set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); + end + % Add the stimulation pair and atlas scout label as the subplot title + AddFastgraphTitle(sInputs, sSortedFastgraphLocIdxs.All(iSubplot), chanNamesSeeg, atlasScoutLabelsSeeg); + end + % Apply the shared y-axis limits to all Fastgraph subplots + for iSubplot = 1:nSubplots-1 + axSubplots(iSubplot).YLim = commonAxisLimits(3:4); + end + % Link subplot axes so that zooming stays synchronized + linkaxes(axSubplots) + set(gcf,'units','normalized','outerposition',[0 0 1 1]) + zoom on + + % === Use the final subplot to display legend === + bst_progress('text', 'Plotting legend...'); + % Generate a cortex snapshot with atlas scout for display + imgCortex = GenerateCortexSnapshot(sSubject, OPTIONS); + % Create the legend subplot with the same spacing settings + subtightplot(nRows, nCols, iSubplot+1, gap, horzMargin, vertMargin); + % Plot the reference panel with the cortex snapshot and axis labels + axSubplots(iSubplot+1) = gca; + PlotLegend(axSubplots(iSubplot+1), imgCortex, round(axSubplots(1).XLim), [0 1], 'Time (ms)', 'Voltage (mV)'); + + % Close progress + bst_progress('stop'); +end + +%% ===== GET STIMULATION SITE CONTACT LOCATION ===== +% Get the midpoint location of each stimulation pair from channel +function stimLocs = GetStimLocs(sInputs, ChannelMat) + % Preallocate one [x y z] midpoint per stimulation pair + stimLocs = zeros(numel(sInputs), 3); + % Get channel names once for lookup + chanNames = {ChannelMat.Channel.Name}; + + % Loop over all stimulation entries + for k = 1:numel(sInputs) + % Split the comment into the two parts + parts = strsplit(sInputs(k).Comment, '-'); + if numel(parts) ~= 2 + continue; + end + % Clean extracted comment + contact1 = parts{1}; + contact2 = parts{2}; + % Get the contact names + contact1Parts = strsplit(contact1); + contact2Parts = strsplit(contact2); + contact1 = contact1Parts{end}; + contact2 = contact2Parts{1}; + % Find the channel indices + iContact1 = find(strcmp(chanNames, contact1), 1); + iContact2 = find(strcmp(chanNames, contact2), 1); + % Compute midpoint only if both contacts exist + if ~isempty(iContact1) && ~isempty(iContact2) + loc1 = ChannelMat.Channel(iContact1).Loc(:)'; + loc2 = ChannelMat.Channel(iContact2).Loc(:)'; + stimLocs(k, :) = (loc1 + loc2) / 2; + end + end +end + +%% ===== LAPRAP STYLE LOCATION SORTING ===== +% Get indices of location sorted in (L)eft side (A)nterior to (P)osterior (LAP), +% then (R)ight side (A)nterior to (P)osterior (RAP) style given the contact locations +% +% Contacts are first separated into left and right hemispheres using the +% y coordinate (left: y >= 0, right: y < 0). Within each hemisphere, contacts +% are ordered by x-coordinate in descending order. +% +% Repeated locations (when there are multiple recordings from the same stimulation site) +% are handled safely by using the original row index as a secondary sorting key. +% This keeps identical locations grouped together while preserving their original input order. +% +% Contacts exactly on the midline (y == 0) are assigned to the left hemisphere. +function sSortedLocIdxs = SortLAPRAP(contactLocs) + % Initialize output structure + sSortedLocIdxs = struct(); + % Original row index of each location + contactIdxs = (1:size(contactLocs, 1))'; + % Append original row indices so duplicate coordinates keep input order + contactLocsWithIdx = [contactLocs, contactIdxs]; + % Identify contacts in the left and right hemispheres + isLeftHemisphere = contactLocsWithIdx(:, 2) >= 0; + isRightHemisphere = ~isLeftHemisphere; + % Extract contact locations for each hemisphere + leftContactLocs = contactLocsWithIdx(isLeftHemisphere, :); + rightContactLocs = contactLocsWithIdx(isRightHemisphere, :); + % Sort left and right hemisphere contacts by x-coordinate in descending order (-xCoordColumn). + % Use original index (idxColumn) as a secondary key so repeated locations remain grouped + % and keep their original input order + xCoordColumn = 1; + idxColumn = 4; + leftContactLocs = sortrows(leftContactLocs, [xCoordColumn, idxColumn], {'descend' 'ascend'}); + rightContactLocs = sortrows(rightContactLocs, [xCoordColumn, idxColumn], {'descend' 'ascend'}); + % Store sorted original indices for each hemisphere + sSortedLocIdxs.Left = leftContactLocs(:, 4)'; + sSortedLocIdxs.Right = rightContactLocs(:, 4)'; + % Combined sorted indices + sSortedLocIdxs.All = [sSortedLocIdxs.Left, sSortedLocIdxs.Right]; +end + +%% ===== LOAD AND FILTER SEEG DATA ===== +% Load each selected SEEG block and optionally exclude contacts based on +% distance from the stimulation site +function [seegData, excludedContacts] = GetSeegData(sInputs, sSortedFastgraphLocIdxs, stimLocs, ChannelMat, OPTIONS) + % Intialize output + seegData = cell(numel(sInputs), 1); + excludedContacts = cell(numel(sInputs), 1); + % Get index of SEEG channel types + iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); + for k = 1:numel(sInputs) + % Load current file + data = load(file_fullpath(sInputs(sSortedFastgraphLocIdxs.All(k)).FileName)); + % Mark bad channels as NaN + data.F(data.ChannelFlag<0, :) = NaN; + if ~isempty(stimLocs) + % Current stimulation center + stimCenter = stimLocs(sSortedFastgraphLocIdxs.All(k), :); + % Compute distance from stimulation site to each SEEG contact (mm) + contactDist = zeros(1, numel(ChannelMat.Channel)); + for j = iSeeg + contactDist(j) = norm(stimCenter - ChannelMat.Channel(j).Loc', 2) * 1000; + end + % Exclude stimulation contacts themselves + iStimContacts = (contactDist > 0) & (contactDist <= 2); + % Exclude contacts within user-provided distance from the stimulation sites + % iExcluded = contactDist > OPTIONS.ExcludeRadius; + iExcluded = (contactDist > 2) & (contactDist <= OPTIONS.ExcludeRadius); + % Keep only valid SEEG contacts + isSeeg = strcmp('SEEG',{ChannelMat.Channel.Type}); + validContacts = isSeeg & ~iExcluded & ~iStimContacts; + excludedContacts{k} = ~validContacts; + % Report excluded contacts + fprintf('Contacts excluded for being within the %d mm exclusion zone "%s":\n', OPTIONS.ExcludeRadius, sInputs(k).Comment); + fprintf('%s %s ', ChannelMat.Channel(iStimContacts).Name, ChannelMat.Channel(iExcluded).Name); + fprintf('\n\n'); + % Remove excluded channels + data.F(excludedContacts{k}, :) = NaN; + else + % If no stimulation locations are available, keep only SEEG channels + excludedContacts{k} = ~iSeeg; + end + % Store SEEG data for the current block + seegData{k} = data; + end +end + +%% ===== SPLIT CONTACTS TO LEFT/RIGHT HEMISPHERE ===== +% Split SEEG contacts into left and right hemisphere groups +function sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat) + % Initialize output structure + sContactGroupLocIdxs = struct(); + % Get index of SEEG channel type + iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); + % Check if valid location data is available + noLocations = isempty(stimLocs) || ~any(stimLocs(:)); + if noLocations + % Use channel group names to assign hemisphere + sContactGroupLocIdxs.Left = zeros(1, length(iSeeg)); + for i = 1:length(iSeeg) + % Left groups start with an apostrophe + sContactGroupLocIdxs.Left(i) = strcmp(ChannelMat.Channel(iSeeg(i)).Group(1), ''''); + end + % Remaining contacts belong to the right hemisphere + sContactGroupLocIdxs.Right = ~sContactGroupLocIdxs.Left; + else + % Store SEEG contact coordinates + contactLocs = zeros(length(iSeeg), 3); + for i = 1:length(iSeeg) + contactLocs(i, :) = ChannelMat.Channel(iSeeg(i)).Loc'; + end + % Use coordinates to split contacts by hemisphere + sContactGroupLocIdxs = SortLAPRAP(contactLocs); + end +end + +%% ===== WITHIN-HEMISPHERE DATA SORTING ===== +% Sort left and right hemisphere channel data within a selected time +% window using either RMS amplitude or maximum absolute amplitude +function sSorted = ApplyDataSorting(subplotData, seegData, OPTIONS) + % Initialize output structure + sSorted = struct(); + % Get sample indices used for sorting + if isempty(OPTIONS.SortWindow) + sortWindowIdx = 1:size(seegData{1}.F,2); + else + sortWindowIdx = OPTIONS.SortWindow(1):OPTIONS.SortWindow(2); + end + % Sort channels within each hemisphere using the selected metric + switch OPTIONS.SortMethod + case 1 % Root Mean Square + if ~isempty(subplotData.leftData) + leftDataRms = sqrt(sum(subplotData.leftData(:,sortWindowIdx).^2, 2)); + leftDataRms(isnan(leftDataRms)) = -Inf; + [sSorted.Vals.Left, sSorted.Idxs.Left] = sort(leftDataRms,'ascend'); + end + if ~isempty(subplotData.rightData) + rightDataRms = sqrt(sum(subplotData.rightData(:,sortWindowIdx).^2, 2)); + rightDataRms(isnan(rightDataRms)) = -Inf; + [sSorted.Vals.Right, sSorted.Idxs.Right] = sort(rightDataRms,'ascend'); + end + case 2 % Max Absolute + if ~isempty(subplotData.leftData) + leftDataMax = max(abs(subplotData.leftData(:,sortWindowIdx)),[],2); + leftDataMax(isnan(leftDataMax)) = -Inf; + [sSorted.Vals.Left, sSorted.Idxs.Left] = sort(leftDataMax,1,'ascend'); + end + if ~isempty(subplotData.rightData) + rightDataMax = max(abs(subplotData.rightData(:,sortWindowIdx)),[],2); + rightDataMax(isnan(rightDataMax)) = -Inf; + [sSorted.Vals.Right, sSorted.Idxs.Right] = sort(rightDataMax,1,'ascend'); + end + end +end + +%% ===== PLOT FASTGRAPH ===== +% Create one fastgraph subplot. +% Left-hemisphere SEEG channels are plotted as positive stacked areas +% Right-hemisphere SEEG channels are plotted as negative stacked areas +function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInputs, stimLocs, iSubplot, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) + % Initialize output handles + hLeftAreaPlot = []; + hRightAreaPlot = []; + + % Get cortex to be used for region/color lookup + sSubject = bst_get('Subject', sInputs(1).SubjectName); + CortexFile = sSubject.Surface(sSubject.iCortex).FileName; + sCortex = bst_memory('LoadSurface', CortexFile); + % Resolve selected scouts + selectedScoutLabels = ResolveScoutSelection(sCortex, OPTIONS); + % Get indices of all SEEG channels + iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); + % Match channel names against atlas table names + [~, iChanLocs] = ismember({ChannelMat.Channel.Name}, chanNamesSeeg); + % Check whether stimulation locations are available + hasStimLocs = any(stimLocs(:)); + % Select the time samples to display + plotWindowIdx = OPTIONS.PlotWindow(1) + 101 : OPTIONS.PlotWindow(2) + 101; + timeMs = seegData{iSubplot}.Time(plotWindowIdx) * 1000; + + fprintf('\n===== Fastgraph %d/%d: Stimulation site "%s" =====\n', iSubplot, numel(sInputs), sInputs(iSubplot).Comment) + % Extract SEEG data once for this subplot + Fout = seegData{iSubplot}.F(iSeeg, :); + + % Loop over left and right hemispheres + for iSide = 1:2 + if iSide == 1 + % Left hemisphere settings + if isempty(subplotData.leftData) + continue; + end + sideName = 'Left'; + groupLocIdxs = sContactGroupLocIdxs.Left; + sortedIdxs = sSubplotDataSorted.Idxs.Left; + signFactor = 1; + else + % Right hemisphere settings + if isempty(subplotData.rightData) + continue; + end + sideName = 'Right'; + groupLocIdxs = sContactGroupLocIdxs.Right; + sortedIdxs = sSubplotDataSorted.Idxs.Right; + signFactor = -1; + end + + % Reorder SEEG channels for the current hemisphere + contactIdxs = groupLocIdxs(sortedIdxs); + plotLocs = iSeeg(contactIdxs); + hemiData = abs(Fout(contactIdxs, :)); + % Get atlas scout labels for these channels + channelScoutLabels = cell(1, numel(plotLocs)); + for i = 1:numel(plotLocs) + channelScoutLabels{i} = atlasScoutLabelsSeeg{iChanLocs(plotLocs(i))}; + end + % Filter channels using resolved scout selection + if hasStimLocs + toPlot = ismember(channelScoutLabels, selectedScoutLabels); + else + toPlot = true(1, numel(plotLocs)); + end + % Keep track of number of channel before filtering + nChannelsBeforeFilter = numel(plotLocs); + % Keep only channels that pass the filters + plotLocs = plotLocs(toPlot); + hemiData = hemiData(toPlot, :); + channelScoutLabels = channelScoutLabels(toPlot); + + % Skip plotting if no channels remain + if isempty(plotLocs) + fprintf('\n%s contacts and atlas scout labels:\n', sideName); + if nChannelsBeforeFilter > 0 + fprintf('Nothing to plot. All contacts were excluded by the atlas/scout selection.\n'); + else + fprintf('Nothing to plot. No contacts are available for this hemisphere.\n'); + end + continue; + end + + % Plot stacked area traces for the current hemisphere + hAreaPlot = area(timeMs, signFactor * hemiData(:, plotWindowIdx)'); + + % Print labels and assign colors + fprintf('\n%s contacts and atlas scout labels:\n', sideName); + isAllContactsExcluded = 1; + for i = 1:numel(plotLocs) + atlasScoutLabelSeeg = channelScoutLabels{i}; + if ~excludedContacts{iSubplot}(plotLocs(i)) + fprintf('%s - %s\n', ChannelMat.Channel(plotLocs(i)).Name, atlasScoutLabelSeeg); + isAllContactsExcluded = 0; + end + region = GetRegionFromScouts(sCortex, atlasScoutLabelSeeg, OPTIONS); + hAreaPlot(i).FaceColor = region.Color; + end + if isAllContactsExcluded + fprintf('Nothing plotted. All contacts lie in the excluded region.\n'); + end + % Store handles in the correct output variable + if iSide == 1 + hLeftAreaPlot = hAreaPlot; + % Keep current plot so right side plot can be added + hold on; + else + hRightAreaPlot = hAreaPlot; + % Release the hold state after plotting both sides + hold off; + end + end +end + +%% ===== ATLAS REGION FROM SCOUTS ===== +% Map an atlas scout label to a Brainstorm region code and plot color +function region = GetRegionFromScouts(sCortex, inputAtlasScoutLabel, OPTIONS) + % Default output if no matching scout is found + region.Name = '?'; + region.Color = [0.5 0.5 0.5]; + % Find the atlas selected by the user + iAtlas = find(strcmpi({sCortex.Atlas.Name}, OPTIONS.Atlas), 1); + if isempty(iAtlas) + return; + end + % Get the selected atlas + atlas = sCortex.Atlas(iAtlas); + % Match the input atlas scout label against atlas scouts + for iScout = 1:numel(atlas.Scouts) + atlasScoutLabel = atlas.Scouts(iScout).Label(1:end-2); + if ~isempty(strfind(lower(inputAtlasScoutLabel), lower(atlasScoutLabel))) + % Matching scout found: assign region name + region.Name = atlas.Scouts(iScout).Region(2:end); + % Assign color based on the selected color scheme + if strcmp(OPTIONS.ColorScheme, 'Region') + region.Color = panel_scout('GetRegionColor', atlas.Scouts(iScout).Region); + else + region.Color = atlas.Scouts(iScout).Color; + end + return; + end + end +end + +%% ===== FASTGRAPH TITLE ===== +% Build the title shown above each subplot using the stimulation pair and +% the atlas label associated with the first contact +function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutLabelsSeeg) + % Split the comment into the two parts + parts = strsplit(sInputs(iSortedFastgraph).Comment, '-'); + % Clean extracted comment + contact1 = strtrim(parts{1}); + % Get the contact names + contact1Parts = strsplit(contact1); + contact1 = contact1Parts{end}; + % Look up atlas label for the first contact + iContact1 = find(strcmp(chanNamesSeeg, contact1), 1); + if ~isempty(iContact1) + contact1AtlasScoutLabel = atlasScoutLabelsSeeg{iContact1}; + else + contact1AtlasScoutLabel = '?'; + end + title(sprintf('%s\n%s', sInputs(iSortedFastgraph).Comment, contact1AtlasScoutLabel),'fontsize', 8); +end + +%% ===== RESOLVE SELECTED SCOUTS ===== +% Resolve which atlas scouts should be used based on either: +% 1) explicit scout labels entered by the user, or +% 2) selected anatomical regions from the checkboxes +function [selectedScoutLabels, iSelectedScouts, iAtlas] = ResolveScoutSelection(sCortex, OPTIONS) + % Default outputs + selectedScoutLabels = {}; + iSelectedScouts = []; + iAtlas = []; + % Find selected atlas + iAtlas = find(strcmpi({sCortex.Atlas.Name}, OPTIONS.Atlas), 1); + if isempty(iAtlas) + return; + end + atlas = sCortex.Atlas(iAtlas); + % Keep only non-empty scout labels entered in the GUI + enteredLabels = OPTIONS.AtlasScoutLabels(~cellfun(@isempty, OPTIONS.AtlasScoutLabels)); + if ~isempty(enteredLabels) + % Explicit scout-label filtering + isKeep = ismember({atlas.Scouts.Label}, enteredLabels); + else + % Region-based filtering + allRegionCodes = {'PF','F','C','P','T','O','L'}; + selectedRegions = allRegionCodes(OPTIONS.Region); + % Remove the leading character from Brainstorm scout region code + scoutRegions = cellfun(@(x) x(2:end), {atlas.Scouts.Region}, 'UniformOutput', false); + isKeep = ismember(scoutRegions, selectedRegions); + end + % Return selected scout indices and labels + iSelectedScouts = find(isKeep); + selectedScoutLabels = {atlas.Scouts(iSelectedScouts).Label}; +end + +%% ===== GENERATE IMAGE FOR LEGEND ===== +% Render the cortex surface with only the scouts selected from the GUI and +% color them either by region or by label +function imgCortex = GenerateCortexSnapshot(sSubject, OPTIONS) + % Default output + imgCortex = []; + % Load cortex + CortexFile = sSubject.Surface(sSubject.iCortex).FileName; + sCortex = bst_memory('LoadSurface', CortexFile); + % Resolve selected scouts from GUI options + [~, iSelectedScouts, iAtlas] = ResolveScoutSelection(sCortex, OPTIONS); + if isempty(iAtlas) || isempty(iSelectedScouts) + return; + end + % Open cortex figure + hFigSurf = view_surface(CortexFile); + figure_3d('SetStandardView', hFigSurf, 'left'); + bst_figures('SetBackgroundColor', hFigSurf, [1 1 1]); + % Select atlas + panel_scout('SetCurrentAtlas', iAtlas); + % Set options + switch(OPTIONS.ColorScheme) + case 'Region' + panel_scout('SetScoutsOptions', 0, 0, 1, 'select', 0, 1, 0, 1); + case 'Label' + panel_scout('SetScoutsOptions', 0, 0, 1, 'select', 0, 1, 0, 0); + end + % Show only selected scouts + panel_scout('SetSelectedScouts', iSelectedScouts); + % Capture image + img = out_figure_image(hFigSurf); + % Crop background + bgColor = img(1,1,:); + mask = (img(:,:,1) == bgColor(1)) & ... + (img(:,:,2) == bgColor(2)) & ... + (img(:,:,3) == bgColor(3)); + goodRows = any(~mask, 2); + goodCols = any(~mask, 1); + imgCortex = img(goodRows, goodCols, :); + % Close figure + close(hFigSurf); +end + +%% ===== PLOT LEGEND ===== +% Shows the legend for the Fastgraph plots as in the paper +function PlotLegend(axSubplotLegend, brainImg, xRange, yRange, xLabel, yLabel) + % === Prepare the plot area === + % Set the visible x- and y-axis limits + set(axSubplotLegend, 'XLim', xRange, 'YLim', yRange); + % Add x-axis label + axSubplotLegend.XLabel.String = xLabel; + % Move x-axis label closer to the axis (slightly upward) + axSubplotLegend.XLabel.Position = [mean(axSubplotLegend.XLim), axSubplotLegend.YLim(1) - 0.01, 0]; + % Add y-axis label + axSubplotLegend.YLabel.String = yLabel; + % Move the y-axis label closer to the axis (slightly right) + axSubplotLegend.YLabel.Position = [axSubplotLegend.XLim(1) - 5, mean(axSubplotLegend.YLim), 0]; + % Show ticks only at the minimum and maximum values of each axis + axSubplotLegend.XTick = [xRange(1), xRange(2)]; + axSubplotLegend.YTick = [yRange(1), yRange(2)]; + + % === Create overlay axes for the brain atlas image === + axImg = axes('Parent', ancestor(axSubplotLegend, 'figure'), ... + 'Units', 'pixels', ... + 'Color', 'none'); + % Display the brain image inside the overlay axes + hImg = imshow(brainImg, 'Parent', axImg); + % Hide the overlay axes so only the image is visible + axis(axImg, 'off'); + % Keep the original axes limits fixed so the image does not alter them + axis(axSubplotLegend, 'manual'); + % Initial placement + UpdateLegendImage(axSubplotLegend, axImg, brainImg); + % Update placement whenever the figure is resized/moved + hFig = ancestor(axSubplotLegend, 'figure'); + hFig.SizeChangedFcn = @(~,~)UpdateLegendImage(axSubplotLegend, axImg, brainImg); + + % Add left/right hemisphere labels with pixel-based spacing + AddLegendHemisphereLabels(axSubplotLegend, xRange, yRange); +end + +%% ===== ADD 'L/R' HEMISPHERE LABELS IN THE LEGEND ===== +% Add 'L/R' hemisphere labels to the legend axes +function AddLegendHemisphereLabels(axSubplotLegend, xRange, yRange) + % Position labels near the right side of the legend axes + xSpan = diff(xRange); + ySpan = diff(yRange); + xLR = xRange(2) - 0.08 * xSpan; + + % Get axes height in pixels + oldUnits = axSubplotLegend.Units; + axSubplotLegend.Units = 'pixels'; + axPos = axSubplotLegend.Position; + axSubplotLegend.Units = oldUnits; + + % Convert a fixed pixel gap into data units + pixelsPerDataY = axPos(4) / ySpan; + gapPx = max(14, axSubplotLegend.FontSize + 4); + gapData = gapPx / pixelsPerDataY; + + % Place labels above and below the x-axis + yAxisLevel = yRange(1); + yL = yAxisLevel + gapData; + yR = yAxisLevel - gapData; + + % Draw the labels + text(axSubplotLegend, xLR, yL, 'L', ... + 'FontSize', 8, ... + 'FontWeight', 'bold', ... + 'HorizontalAlignment', 'right', ... + 'VerticalAlignment', 'middle', ... + 'Clipping', 'off', ... + 'Margin', 1); + + text(axSubplotLegend, xLR, yR, 'R', ... + 'FontSize', 8, ... + 'FontWeight', 'bold', ... + 'HorizontalAlignment', 'right', ... + 'VerticalAlignment', 'middle', ... + 'Clipping', 'off', ... + 'Margin', 1); +end + +%% ===== UPDATE LEGEND IMAGE ===== +% Update the overlay image position so it stays centered inside the +% legend subplot when the figure is resized or moved across screens +function UpdateLegendImage(axSubplotLegend, axImg, brainImg) + % Get original image size in pixels + imgH = size(brainImg, 1); + imgW = size(brainImg, 2); + % Read the legend subplot position in pixel units + oldUnits = axSubplotLegend.Units; + axSubplotLegend.Units = 'pixels'; + % Get the axes position in pixel units: [left, bottom, width, height] + pos = axSubplotLegend.Position; + axSubplotLegend.Units = oldUnits; + % Available subplot width and height in pixels + boxW = pos(3); + boxH = pos(4); + % Scale the image to fit inside the subplot while preserving aspect ratio + scale = min(boxW / imgW, boxH / imgH) * 0.75; + newW = imgW * scale; + newH = imgH * scale; + % Center the image inside the legend subplot + xLeft = pos(1) + (boxW - newW) / 2; + yBottom = pos(2) + (boxH - newH) / 2; + % Update the overlay axes position in pixel coordinates + axImg.Units = 'pixels'; + axImg.Position = [xLeft, yBottom, newW, newH]; +end \ No newline at end of file From fe624c21bfeba32a64aa75d29cb61d6c8a7b77ef Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 1 May 2026 01:53:38 -0700 Subject: [PATCH 04/55] Add Fastgraph tutorial script --- toolbox/script/tutorial_fastgraph.m | 251 ++++++++++++++++++++++++++++ 1 file changed, 251 insertions(+) create mode 100644 toolbox/script/tutorial_fastgraph.m diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m new file mode 100644 index 0000000000..5d803188f6 --- /dev/null +++ b/toolbox/script/tutorial_fastgraph.m @@ -0,0 +1,251 @@ +function tutorial_fastgraph(tutorial_dir, reports_dir) +% TUTORIAL_FASTGRAPH: Script that reproduces the results of the online tutorial "Fastgraph". +% +% CORRESPONDING ONLINE TUTORIALS: +% https://neuroimage.usc.edu/brainstorm/Tutorials/FastGraph +% +% INPUTS: +% - tutorial_dir : Directory where the tutorial_fastgraph.zip file has been unzipped +% - reports_dir : Directory where to save the execution report (instead of displaying it) + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c) University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% 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: Chinmay Chinara, 2026 +% John C. Mosher, 2026 + +%% ===== PARSE INPUTS ===== +% Output folder for reports +if (nargin < 2) || isempty(reports_dir) || ~isfolder(reports_dir) + reports_dir = []; +end +% You have to specify the folder in which the tutorial dataset is unzipped +if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir) + error('The first argument must be the full path to the tutorial dataset folder.'); +end +% Subject name +SubjectName = 'Subject01'; + +%% ===== FILES TO IMPORT ===== +% Build the path of the files to import +tutorial_dir = bst_fullfile(tutorial_dir, 'tutorial_fastgraph'); +MriFilePre = bst_fullfile(tutorial_dir, 'anatomy', 'pre_T1.nii.gz'); +MriCat12Path = fullfile(tutorial_dir, 'anatomy', 'cat12'); +BaselineFile = bst_fullfile(tutorial_dir, 'recordings', 'Baseline.edf'); +ElecPosFile = bst_fullfile(tutorial_dir, 'recordings', 'Subject01_electrodes_mm.tsv'); +% Check if the folder contains the required files +if ~file_exist(BaselineFile) + error(['The folder ' tutorial_dir ' does not contain the folder from the file tutorial_fastgraph.zip.']); +end +isMriSegmented = file_exist(bst_fullfile(MriCat12Path, 'Subject01.nii')); + +%% ===== CREATE PROTOCOL ===== +% The protocol name has to be a valid folder name (no spaces, no weird characters...) +ProtocolName = 'TutorialFastgraph'; +% Start brainstorm without the GUI +if ~brainstorm('status') + brainstorm nogui +end +% Delete existing protocol +gui_brainstorm('DeleteProtocol', ProtocolName); +% Create new protocol +gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); +% Start a new report +bst_report('Start'); + +%% ===== IMPORT MRI AND CT VOLUMES ===== +if ~isMriSegmented + % Process: Import MRI + bst_process('CallProcess', 'process_import_mri', [], [], ... + 'subjectname', SubjectName, ... + 'voltype', 'mri', ... + 'comment', 'pre_T1', ... + 'mrifile', {MriFilePre, 'ALL'}, ... + 'nas', [107, 176, 105], ... + 'lpa', [ 34, 89, 74], ... + 'rpa', [175, 89, 74]); + % Process: Segment MRI with CAT12 + bst_process('CallProcess', 'process_segment_cat12', [], [], ... + 'subjectname', SubjectName, ... + 'nvertices', 15000, ... + 'tpmnii', {'', 'Nifti1'}, ... + 'sphreg', 1, ... % Use spherical registration + 'vol', 0, ... % No volume parcellations + 'extramaps', 0, ... % No additional cortical maps + 'cerebellum', 0); +else + % Process: Import anatomy folder + bst_process('CallProcess', 'process_import_anatomy', [], [], ... + 'subjectname', SubjectName, ... + 'mrifile', {MriCat12Path, 'CAT12'}, ... + 'nvertices', 15000, ... + 'nas', [107, 176, 105], ... + 'lpa', [ 34, 89, 74], ... + 'rpa', [175, 89, 74]); +end +% Get filename for imported volumes +[sSubject, iSubject] = bst_get('Subject', SubjectName); +% Reference MRI +DbMriFilePre = sSubject.Anatomy(sSubject.iAnatomy).FileName; + +%% ===== CREATE SEEG CONTACT IMPLANTATION ===== +iStudyImplantation = db_add_condition(SubjectName, 'Implantation'); +% Import locations and convert to subject coordinate system (SCS) +ImplantationChannelFile = import_channel(iStudyImplantation, ElecPosFile, 'BIDS-SCANRAS-MM', 1, 0, 1, 0, 2, DbMriFilePre); +% Snapshot: SEEG electrodes in MRI slices +hFigMri3d = view_channels_3d(ImplantationChannelFile, 'SEEG', 'anatomy', 1, 0); +bst_report('Snapshot', hFigMri3d, ImplantationChannelFile, 'SEEG electrodes in 3D MRI slices'); +close(hFigMri3d); + +%% ===== ACCESS THE RECORDINGS ===== +% Process: Create link to raw file +sFileRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ... + 'subjectname', SubjectName, ... + 'datafile', {BaselineFile, 'EEG-EDF'}, ... + 'channelreplace', 0, ... + 'channelalign', 0); +% Process: Add EEG positions +bst_process('CallProcess', 'process_channel_addloc', sFileRaw, [], ... + 'channelfile', {ImplantationChannelFile, 'BST'}, ... + 'fixunits', 0, ... % No automatic fixing of distance units required + 'vox2ras', 0); % Do not use the voxel=>subject transformation, already in SCS +% Process: Add EEG positions +bst_process('CallProcess', 'process_channel_addloc', sFileRaw, [], ... + 'channelfile', {ImplantationChannelFile, 'BST'}, ... + 'fixunits', 0, ... % No automatic fixing of distance units required + 'vox2ras', 0); % Do not use the voxel=>subject transformation, already in SCS + +% Process: Customize SPES +bst_process('CallProcess', 'process_customize_spes_nk', sFileRaw, [], ... + 'stimstartlabel', 'SB', ... + 'stimstoplabel', 'SE', ... + 'stimchan', 'DC10', ... + 'stimlabel', 'STIM', ... + 'buffertime', 2, ... % in s + 'offset', -0.001, ... % in ms + 'evtaddoddeven', 1); + +% Process: Load the Stim Start blocks +sFilesStimStart = bst_process('CallProcess', 'process_import_data_event', sFileRaw, [], ... + 'subjectname', SubjectName, ... + 'condition', '', ... + 'eventname', 'SB', ... + 'epochtime', [-2 32], ... % in s + 'createcond', 0, ... + 'ignoreshort', 1, ... + 'usectfcomp', 1, ... + 'usessp', 1, ... + 'freq', [], ... + 'baseline', []); + +% Process: Remove SPES artifacts +sFilesStimStartClean = bst_process('CallProcess', 'process_remove_spes_artifacts', sFilesStimStart, [], ... + 'stimevent', 'STIM', ... + 'cutoff', 2, ... + 'timeart', 0.005, ... % in ms + 'timespline', 0.003); % in ms + +% Process: Load the ODD events +sFilesOdd = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartClean, [], ... + 'subjectname', SubjectName, ... + 'condition', '', ... + 'eventname', 'ODD', ... + 'timewindow', [-2 32], ... % in s + 'epochtime', [-0.100 0.900], ... % in ms + 'createcond', 0, ... + 'ignoreshort', 1, ... + 'usectfcomp', 1, ... + 'usessp', 1, ... + 'freq', [], ... + 'baseline', []); + +% Process: Load the EVEN events +sFilesEven = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartClean, [], ... + 'subjectname', SubjectName, ... + 'condition', '', ... + 'eventname', 'EVEN', ... + 'timewindow', [-2 32], ... % in s + 'epochtime', [-0.100 0.900], ... % in ms + 'createcond', 0, ... + 'ignoreshort', 1, ... + 'usectfcomp', 1, ... + 'usessp', 1, ... + 'freq', [], ... + 'baseline', []); + +% Process: Get ODD average (by trial group) +sFilesAvgOdd = bst_process('CallProcess', 'process_average', sFilesOdd, [], ... + 'avgtype', 5, ... % Trial group (folder average) + 'avg_func', 1, ... % Arithmetic average: mean(x) + 'weighted', 0, ... + 'keepevents', 0); + +% Process: Get ODD average (by trial group) +sFilesAvgEven = bst_process('CallProcess', 'process_average', sFilesEven, [], ... + 'avgtype', 5, ... % Trial group (folder average) + 'avg_func', 1, ... % Arithmetic average: mean(x) + 'weighted', 0, ... + 'keepevents', 0); + +% Process: Average the ODD and EVEN +sFilesFastgraph = {}; +for i = 1:length(sFilesAvgOdd) + sFileAvg = bst_process('CallProcess', 'process_average', [sFilesAvgOdd(i), sFilesAvgEven(i)], [], ... + 'avgtype', 1, ... % Everything + 'avg_func', 1, ... % Arithmetic average: mean(x) + 'weighted', 0, ... + 'keepevents', 0); + CommentMat = in_bst_data(sFileAvg.FileName, 'Comment'); + % Get the stim site info + stimSiteInfo = strtrim(strrep(sFilesAvgOdd(i).Comment,'Avg: ODD','')); + stimSiteInfo = strtrim(stimSiteInfo(1:(end-10))); % remove file count + CommentMat.Comment = ['Avg: ' stimSiteInfo]; + % Save changes + bst_save(file_fullpath(sFileAvg.FileName), CommentMat,'v7', 1); + % Register output + sFilesFastgraph{end+1} = sFileAvg.FileName; +end +db_reload_conditions(iSubject); + +% Process: Plot Fastgraph +bst_process('CallProcess', 'process_fastgraph', sFilesFastgraph, [], ... + 'atlas', 'Desikan-Killiany', ... + 'colorscheme', 'Region', ... % Color figures by region + 'regionprefrontal', 1, ... + 'regionfrontal', 1, ... + 'regioncentral', 1, ... + 'regionparietal', 1, ... + 'regiontemporal', 1, ... + 'regionoccipital', 1, ... + 'regionlimbic', 1, ... + 'atlasscoutlabels', '', ... + 'sortmethod', 1, ... % "Root Mean Square" to sort data + 'sortwindow', [0.060, 0.250], ... % Range to sort the data (in ms) + 'plotwindow', [-0.100, 0.900], ... % Plot window + 'edgealpha', 0.05, ... % Edge transparency of plot + 'excluderadius', 20); % Exclusion zone radius + +%% ===== SAVE AND DISPLAY REPORT ===== +ReportFile = bst_report('Save', []); +if ~isempty(reports_dir) && ~isempty(ReportFile) + bst_report('Export', ReportFile, reports_dir); +else + bst_report('Open', ReportFile); +end + +disp([10 'DEMO> Fastgraph tutorial completed' 10]); \ No newline at end of file From 1bbc413a77cde0559d5676fcc12c2da3290dbfd6 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 1 May 2026 10:27:47 -0700 Subject: [PATCH 05/55] Tutorial: Use `regexp` for getting stim site info --- toolbox/script/tutorial_fastgraph.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index 5d803188f6..179b3090b2 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -212,9 +212,9 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'keepevents', 0); CommentMat = in_bst_data(sFileAvg.FileName, 'Comment'); % Get the stim site info - stimSiteInfo = strtrim(strrep(sFilesAvgOdd(i).Comment,'Avg: ODD','')); - stimSiteInfo = strtrim(stimSiteInfo(1:(end-10))); % remove file count - CommentMat.Comment = ['Avg: ' stimSiteInfo]; + % Example: "Avg: ODD A'2-A'3 4.0 #2 (8 files)" > "A'2-A'3 4.0 #2" + stimSiteInfo = regexp(sFilesAvgOdd(i).Comment, '^Avg:\s+\w+\s+(.+?)\s*\(.*\)$', 'tokens', 'once'); + CommentMat.Comment = ['Avg: ' stimSiteInfo{1}]; % Save changes bst_save(file_fullpath(sFileAvg.FileName), CommentMat,'v7', 1); % Register output From 62a04579d102bce09704e5bd993d27196a180a0d Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 1 May 2026 10:40:07 -0700 Subject: [PATCH 06/55] Tutorial: Add comments --- toolbox/script/tutorial_fastgraph.m | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index 179b3090b2..42ef319373 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -156,7 +156,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) % Process: Remove SPES artifacts sFilesStimStartClean = bst_process('CallProcess', 'process_remove_spes_artifacts', sFilesStimStart, [], ... 'stimevent', 'STIM', ... - 'cutoff', 2, ... + 'cutoff', 2, ... % in Hz 'timeart', 0.005, ... % in ms 'timespline', 0.003); % in ms @@ -188,21 +188,21 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'freq', [], ... 'baseline', []); -% Process: Get ODD average (by trial group) +% Process: Average of only ODDs (by trial group) sFilesAvgOdd = bst_process('CallProcess', 'process_average', sFilesOdd, [], ... 'avgtype', 5, ... % Trial group (folder average) 'avg_func', 1, ... % Arithmetic average: mean(x) 'weighted', 0, ... 'keepevents', 0); -% Process: Get ODD average (by trial group) +% Process: Average of only evens (by trial group) sFilesAvgEven = bst_process('CallProcess', 'process_average', sFilesEven, [], ... 'avgtype', 5, ... % Trial group (folder average) 'avg_func', 1, ... % Arithmetic average: mean(x) 'weighted', 0, ... 'keepevents', 0); -% Process: Average the ODD and EVEN +% Process: Average the ODD and EVEN per stimulation site per session sFilesFastgraph = {}; for i = 1:length(sFilesAvgOdd) sFileAvg = bst_process('CallProcess', 'process_average', [sFilesAvgOdd(i), sFilesAvgEven(i)], [], ... @@ -235,8 +235,8 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'regionlimbic', 1, ... 'atlasscoutlabels', '', ... 'sortmethod', 1, ... % "Root Mean Square" to sort data - 'sortwindow', [0.060, 0.250], ... % Range to sort the data (in ms) - 'plotwindow', [-0.100, 0.900], ... % Plot window + 'sortwindow', [0.060, 0.250], ... % Range (middle latency) to sort the data (in ms) + 'plotwindow', [-0.100, 0.900], ... % Plot window (in ms) 'edgealpha', 0.05, ... % Edge transparency of plot 'excluderadius', 20); % Exclusion zone radius From 0142d57ddcfca6e5dc5016b14920f6cf53b4a484 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 11 Jun 2026 13:13:38 -0700 Subject: [PATCH 07/55] Move to subgroup `FAST graph` --- toolbox/process/functions/process_fastgraph.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 98901cb2f6..b5f8a8bbfd 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -37,8 +37,9 @@ % Describe the process and its UI options sProcess.Comment = 'Plot Fastgraphs'; sProcess.Category = 'Custom'; -sProcess.SubGroup = 'Stimulation'; -sProcess.Index = 1100; +sProcess.SubGroup = 'FAST graph'; +sProcess.Index = 1303; +sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/FastGraph'; % Definition of the input accepted by this process sProcess.InputTypes = {'data'}; sProcess.OutputTypes = {'data'}; From 005acb634fa64f0a77d2fe23d711b0a31d46322d Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 17 Jun 2026 00:01:32 -0700 Subject: [PATCH 08/55] Use type as `scout` instead of `atlas` --- toolbox/process/functions/process_fastgraph.m | 22 +++++++------------ 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index b5f8a8bbfd..8c253424b5 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -45,10 +45,10 @@ sProcess.OutputTypes = {'data'}; sProcess.nInputs = 1; sProcess.nMinFiles = 1; -% Atlas to use for plotting Fastgraph -sProcess.options.atlas.Comment = 'Atlas to plot: '; -sProcess.options.atlas.Type = 'atlas'; -sProcess.options.atlas.Value = []; +% Scouts to use for plotting Fastgraph +sProcess.options.scouts.Comment = ''; +sProcess.options.scouts.Type = 'scout'; +sProcess.options.scouts.Value = {}; % Color Fastgraph by region or by label sProcess.options.label2.Comment = 'Color Fastgraph by region or by label ?'; sProcess.options.label2.Type = 'label'; @@ -79,12 +79,6 @@ sProcess.options.regionlimbic.Comment = '7: Limbic'; sProcess.options.regionlimbic.Type = 'checkbox'; sProcess.options.regionlimbic.Value = 1; -% Atlas scout labels to plot -sProcess.options.label4.Comment = 'For multiple labels: separate them with commas'; -sProcess.options.label4.Type = 'label'; -sProcess.options.atlasscoutlabels.Comment = 'Atlas scout labels to plot: '; -sProcess.options.atlasscoutlabels.Type = 'text'; -sProcess.options.atlasscoutlabels.Value = ''; % Add separator sProcess.options.separator1.Type = 'separator'; % Method for sorting the data @@ -131,8 +125,10 @@ %% ===== GET OPTIONS ===== function OPTIONS = GetOptions(sProcess) OPTIONS = struct(); - % Atlas to use for plotting Fastgraph - OPTIONS.Atlas = sProcess.options.atlas.Value; + % Atlas and scouts to use for plotting Fastgraph + ScoutsList = sProcess.options.scouts.Value; + OPTIONS.Atlas = ScoutsList{1,1}; + OPTIONS.AtlasScoutLabels = ScoutsList{1,2}; % Color figure by region or by label OPTIONS.ColorScheme = sProcess.options.colorscheme.Value; % Select regions to include @@ -143,8 +139,6 @@ sProcess.options.regiontemporal.Value sProcess.options.regionoccipital.Value sProcess.options.regionlimbic.Value]); - % Atlas scout labels to plot - OPTIONS.AtlasScoutLabels = strtrim(strsplit(sProcess.options.atlasscoutlabels.Value,',')); % Method for sorting the data OPTIONS.SortMethod = sProcess.options.sortmethod.Value; % Sort window From 8b3f6287fe697605c09251b36197493418eebc62 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 17 Jun 2026 00:02:03 -0700 Subject: [PATCH 09/55] Clean up --- toolbox/process/functions/process_fastgraph.m | 44 +++++-------------- 1 file changed, 11 insertions(+), 33 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 8c253424b5..838f6389a9 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -173,26 +173,6 @@ return; end - % Get subject - sSubject = bst_get('Subject', sInputs(1).SubjectName); - CortexFile = sSubject.Surface(sSubject.iCortex).FileName; - sCortex = bst_memory('LoadSurface', CortexFile); - - % Get the last used atlas if atlas not selected - if isempty(OPTIONS.Atlas) - OPTIONS.Atlas = sCortex.Atlas(sCortex.iAtlas).Name; - end - % Find the atlas selected by the user - iAtlas = find(strcmpi({sCortex.Atlas.Name}, OPTIONS.Atlas), 1); - - % Early exit if any entered atlas scout label does not exist - allAtlasScoutLabels = {sCortex.Atlas(iAtlas).Scouts.Label}; - enteredLabels = OPTIONS.AtlasScoutLabels(~cellfun(@isempty, OPTIONS.AtlasScoutLabels)); - if ~isempty(enteredLabels) && ~all(ismember(enteredLabels, allAtlasScoutLabels)) - bst_report('Error', sProcess, [], 'One or more scout labels entered are not present in the selected atlas'); - return; - end - % Load the channel file ChannelFile = file_fullpath(sInputs(1).ChannelFile); ChannelMat = load(ChannelFile); @@ -290,7 +270,7 @@ % === Use the final subplot to display legend === bst_progress('text', 'Plotting legend...'); % Generate a cortex snapshot with atlas scout for display - imgCortex = GenerateCortexSnapshot(sSubject, OPTIONS); + imgCortex = GenerateCortexSnapshot(sInputs, OPTIONS); % Create the legend subplot with the same spacing settings subtightplot(nRows, nCols, iSubplot+1, gap, horzMargin, vertMargin); % Plot the reference panel with the cortex snapshot and axis labels @@ -563,11 +543,11 @@ hemiData = hemiData(toPlot, :); channelScoutLabels = channelScoutLabels(toPlot); - % Skip plotting if no channels remain - if isempty(plotLocs) - fprintf('\n%s contacts and atlas scout labels:\n', sideName); + % Skip plotting if no channels remain after atlas/scout filtering + fprintf('\n%s contacts and atlas scout labels:\n', sideName); + if isempty(plotLocs) if nChannelsBeforeFilter > 0 - fprintf('Nothing to plot. All contacts were excluded by the atlas/scout selection.\n'); + fprintf('Nothing to plot. All contacts were filtered out by the selected atlas/scout regions.\n'); else fprintf('Nothing to plot. No contacts are available for this hemisphere.\n'); end @@ -578,7 +558,6 @@ hAreaPlot = area(timeMs, signFactor * hemiData(:, plotWindowIdx)'); % Print labels and assign colors - fprintf('\n%s contacts and atlas scout labels:\n', sideName); isAllContactsExcluded = 1; for i = 1:numel(plotLocs) atlasScoutLabelSeeg = channelScoutLabels{i}; @@ -590,7 +569,7 @@ hAreaPlot(i).FaceColor = region.Color; end if isAllContactsExcluded - fprintf('Nothing plotted. All contacts lie in the excluded region.\n'); + fprintf('Nothing plotted. All contacts lie within the stimulation-site exclusion zone.\n'); end % Store handles in the correct output variable if iSide == 1 @@ -658,7 +637,7 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL %% ===== RESOLVE SELECTED SCOUTS ===== % Resolve which atlas scouts should be used based on either: -% 1) explicit scout labels entered by the user, or +% 1) explicit scout labels selected by the user, or % 2) selected anatomical regions from the checkboxes function [selectedScoutLabels, iSelectedScouts, iAtlas] = ResolveScoutSelection(sCortex, OPTIONS) % Default outputs @@ -671,11 +650,9 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL return; end atlas = sCortex.Atlas(iAtlas); - % Keep only non-empty scout labels entered in the GUI - enteredLabels = OPTIONS.AtlasScoutLabels(~cellfun(@isempty, OPTIONS.AtlasScoutLabels)); - if ~isempty(enteredLabels) + if ~isempty(OPTIONS.AtlasScoutLabels) % Explicit scout-label filtering - isKeep = ismember({atlas.Scouts.Label}, enteredLabels); + isKeep = ismember({atlas.Scouts.Label}, OPTIONS.AtlasScoutLabels); else % Region-based filtering allRegionCodes = {'PF','F','C','P','T','O','L'}; @@ -692,10 +669,11 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL %% ===== GENERATE IMAGE FOR LEGEND ===== % Render the cortex surface with only the scouts selected from the GUI and % color them either by region or by label -function imgCortex = GenerateCortexSnapshot(sSubject, OPTIONS) +function imgCortex = GenerateCortexSnapshot(sInputs, OPTIONS) % Default output imgCortex = []; % Load cortex + sSubject = bst_get('Subject', sInputs(1).SubjectName); CortexFile = sSubject.Surface(sSubject.iCortex).FileName; sCortex = bst_memory('LoadSurface', CortexFile); % Resolve selected scouts from GUI options From f052c0f43fae043ba67e36982cc5b075bc7c5449 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 17 Jun 2026 11:40:16 -0700 Subject: [PATCH 10/55] Remove duplicate `process_channel_addloc` call --- toolbox/script/tutorial_fastgraph.m | 5 ----- 1 file changed, 5 deletions(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index 42ef319373..4fc24ef3c2 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -120,11 +120,6 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'channelreplace', 0, ... 'channelalign', 0); % Process: Add EEG positions -bst_process('CallProcess', 'process_channel_addloc', sFileRaw, [], ... - 'channelfile', {ImplantationChannelFile, 'BST'}, ... - 'fixunits', 0, ... % No automatic fixing of distance units required - 'vox2ras', 0); % Do not use the voxel=>subject transformation, already in SCS -% Process: Add EEG positions bst_process('CallProcess', 'process_channel_addloc', sFileRaw, [], ... 'channelfile', {ImplantationChannelFile, 'BST'}, ... 'fixunits', 0, ... % No automatic fixing of distance units required From 46fbdf8f3ea68e059539208d76b442078221929b Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 17 Jun 2026 11:51:41 -0700 Subject: [PATCH 11/55] `channel_add_loc`: Explicitly check for channel type in the imported external file --- toolbox/sensors/channel_add_loc.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/toolbox/sensors/channel_add_loc.m b/toolbox/sensors/channel_add_loc.m index 223faea1de..83a8e48b18 100644 --- a/toolbox/sensors/channel_add_loc.m +++ b/toolbox/sensors/channel_add_loc.m @@ -151,7 +151,9 @@ function channel_add_loc(iStudies, LocChannelFile, isInteractive, isMni) % If the channel is already considered as EEG, do not change its type, otherwise set it to EEG if ~ismember(ChannelMat.Channel(ic).Type, {'EEG','SEEG','ECOG'}) ChannelMat.Channel(ic).Type = 'EEG'; - elseif ismember(LocChannelMat.Channel(idef).Type, {'SEEG','ECOG'}) + end + % Check for type in the imported external file + if ismember(LocChannelMat.Channel(idef).Type, {'SEEG','ECOG'}) ChannelMat.Channel(ic).Type = LocChannelMat.Channel(idef).Type; end ChannelMat.Channel(ic).Loc = LocChannelMat.Channel(idef).Loc; From e95a1dac16874b987a4b5bad58d4003810411c3d Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 17 Jun 2026 14:01:57 -0700 Subject: [PATCH 12/55] Refactor tutorial script - Rename processes as per the changes in PR #911 and #912 - Adding `ODD` and `EVEN` events not part of the tutorial (left at user's discretion with a note added to the tutorial about it). As recommended by @jcmosher. --- toolbox/script/tutorial_fastgraph.m | 65 ++++++----------------------- 1 file changed, 13 insertions(+), 52 deletions(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index 4fc24ef3c2..16f65147bb 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -126,14 +126,14 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'vox2ras', 0); % Do not use the voxel=>subject transformation, already in SCS % Process: Customize SPES -bst_process('CallProcess', 'process_customize_spes_nk', sFileRaw, [], ... +bst_process('CallProcess', 'process_evt_detect_spes', sFileRaw, [], ... 'stimstartlabel', 'SB', ... 'stimstoplabel', 'SE', ... 'stimchan', 'DC10', ... 'stimlabel', 'STIM', ... 'buffertime', 2, ... % in s 'offset', -0.001, ... % in ms - 'evtaddoddeven', 1); + 'evtaddoddeven', 0); % Process: Load the Stim Start blocks sFilesStimStart = bst_process('CallProcess', 'process_import_data_event', sFileRaw, [], ... @@ -149,31 +149,20 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'baseline', []); % Process: Remove SPES artifacts -sFilesStimStartClean = bst_process('CallProcess', 'process_remove_spes_artifacts', sFilesStimStart, [], ... +sFilesStimStartRmSpes = bst_process('CallProcess', 'process_remove_spes_artifacts', sFilesStimStart, [], ... 'stimevent', 'STIM', ... - 'cutoff', 2, ... % in Hz 'timeart', 0.005, ... % in ms 'timespline', 0.003); % in ms -% Process: Load the ODD events -sFilesOdd = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartClean, [], ... - 'subjectname', SubjectName, ... - 'condition', '', ... - 'eventname', 'ODD', ... - 'timewindow', [-2 32], ... % in s - 'epochtime', [-0.100 0.900], ... % in ms - 'createcond', 0, ... - 'ignoreshort', 1, ... - 'usectfcomp', 1, ... - 'usessp', 1, ... - 'freq', [], ... - 'baseline', []); +% Process: Remove drift EMD +sFilesStimStartEmd = bst_process('CallProcess', 'process_remove_drift_emd', sFilesStimStartRmSpes, [], ... + 'cutoff', 2); % in ms -% Process: Load the EVEN events -sFilesEven = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartClean, [], ... +% Process: Load the STIM events +sFilesStim = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartEmd, [], ... 'subjectname', SubjectName, ... 'condition', '', ... - 'eventname', 'EVEN', ... + 'eventname', 'STIM', ... 'timewindow', [-2 32], ... % in s 'epochtime', [-0.100 0.900], ... % in ms 'createcond', 0, ... @@ -183,43 +172,16 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'freq', [], ... 'baseline', []); -% Process: Average of only ODDs (by trial group) -sFilesAvgOdd = bst_process('CallProcess', 'process_average', sFilesOdd, [], ... - 'avgtype', 5, ... % Trial group (folder average) - 'avg_func', 1, ... % Arithmetic average: mean(x) - 'weighted', 0, ... - 'keepevents', 0); - -% Process: Average of only evens (by trial group) -sFilesAvgEven = bst_process('CallProcess', 'process_average', sFilesEven, [], ... +% Process: Average (by trial group) +sFilesAvg = bst_process('CallProcess', 'process_average', sFilesStim, [], ... 'avgtype', 5, ... % Trial group (folder average) 'avg_func', 1, ... % Arithmetic average: mean(x) 'weighted', 0, ... 'keepevents', 0); -% Process: Average the ODD and EVEN per stimulation site per session -sFilesFastgraph = {}; -for i = 1:length(sFilesAvgOdd) - sFileAvg = bst_process('CallProcess', 'process_average', [sFilesAvgOdd(i), sFilesAvgEven(i)], [], ... - 'avgtype', 1, ... % Everything - 'avg_func', 1, ... % Arithmetic average: mean(x) - 'weighted', 0, ... - 'keepevents', 0); - CommentMat = in_bst_data(sFileAvg.FileName, 'Comment'); - % Get the stim site info - % Example: "Avg: ODD A'2-A'3 4.0 #2 (8 files)" > "A'2-A'3 4.0 #2" - stimSiteInfo = regexp(sFilesAvgOdd(i).Comment, '^Avg:\s+\w+\s+(.+?)\s*\(.*\)$', 'tokens', 'once'); - CommentMat.Comment = ['Avg: ' stimSiteInfo{1}]; - % Save changes - bst_save(file_fullpath(sFileAvg.FileName), CommentMat,'v7', 1); - % Register output - sFilesFastgraph{end+1} = sFileAvg.FileName; -end -db_reload_conditions(iSubject); - % Process: Plot Fastgraph -bst_process('CallProcess', 'process_fastgraph', sFilesFastgraph, [], ... - 'atlas', 'Desikan-Killiany', ... +bst_process('CallProcess', 'process_fastgraph', sFilesAvg, [], ... + 'scouts', {'Desikan-Killiany', {}}, ... 'colorscheme', 'Region', ... % Color figures by region 'regionprefrontal', 1, ... 'regionfrontal', 1, ... @@ -228,7 +190,6 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'regiontemporal', 1, ... 'regionoccipital', 1, ... 'regionlimbic', 1, ... - 'atlasscoutlabels', '', ... 'sortmethod', 1, ... % "Root Mean Square" to sort data 'sortwindow', [0.060, 0.250], ... % Range (middle latency) to sort the data (in ms) 'plotwindow', [-0.100, 0.900], ... % Plot window (in ms) From 90d73429146243ab6a04299148c1032e04e28557 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 18 Jun 2026 14:50:24 -0700 Subject: [PATCH 13/55] Typo: Should be `Hz` --- toolbox/script/tutorial_fastgraph.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index 16f65147bb..698170a413 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -156,7 +156,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) % Process: Remove drift EMD sFilesStimStartEmd = bst_process('CallProcess', 'process_remove_drift_emd', sFilesStimStartRmSpes, [], ... - 'cutoff', 2); % in ms + 'cutoff', 2); % in Hz % Process: Load the STIM events sFilesStim = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartEmd, [], ... From 7f05e98b85cb151ed6263362626cb8b4ca7c978b Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 18 Jun 2026 22:55:25 -0700 Subject: [PATCH 14/55] Use `radio_label` for `SortMethod` --- toolbox/process/functions/process_fastgraph.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 838f6389a9..a3872e4b43 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -84,9 +84,9 @@ % Method for sorting the data sProcess.options.label5.Comment = 'Select method to sort the data:'; sProcess.options.label5.Type = 'label'; -sProcess.options.sortmethod.Comment = {'Root Mean Square', 'Max Absolute'}; -sProcess.options.sortmethod.Type = 'radio'; -sProcess.options.sortmethod.Value = 1; +sProcess.options.sortmethod.Comment = {'Root Mean Square', 'Max Absolute'; 'Root Mean Square', 'Max Absolute'}; +sProcess.options.sortmethod.Type = 'radio_label'; +sProcess.options.sortmethod.Value = 'Root Mean Square'; % Sort window sProcess.options.label6.Comment = 'Choose range to sort over:'; sProcess.options.label6.Type = 'label'; @@ -445,7 +445,7 @@ end % Sort channels within each hemisphere using the selected metric switch OPTIONS.SortMethod - case 1 % Root Mean Square + case 'Root Mean Square' if ~isempty(subplotData.leftData) leftDataRms = sqrt(sum(subplotData.leftData(:,sortWindowIdx).^2, 2)); leftDataRms(isnan(leftDataRms)) = -Inf; @@ -456,7 +456,7 @@ rightDataRms(isnan(rightDataRms)) = -Inf; [sSorted.Vals.Right, sSorted.Idxs.Right] = sort(rightDataRms,'ascend'); end - case 2 % Max Absolute + case 'Max Absolute' if ~isempty(subplotData.leftData) leftDataMax = max(abs(subplotData.leftData(:,sortWindowIdx)),[],2); leftDataMax(isnan(leftDataMax)) = -Inf; From 7e55eea2e4d3a10a265f77e4a74b5f924f10837f Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 18 Jun 2026 23:14:49 -0700 Subject: [PATCH 15/55] Use `list_horizontal` type for regions list --- toolbox/process/functions/process_fastgraph.m | 34 +++---------------- 1 file changed, 4 insertions(+), 30 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index a3872e4b43..a02446f880 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -56,29 +56,9 @@ sProcess.options.colorscheme.Type = 'radio_label'; sProcess.options.colorscheme.Value = 'Region'; % Select regions to include -sProcess.options.label3.Comment = 'Select region(s) to include:'; -sProcess.options.label3.Type = 'label'; -sProcess.options.regionprefrontal.Comment = '1: Prefrontal'; -sProcess.options.regionprefrontal.Type = 'checkbox'; -sProcess.options.regionprefrontal.Value = 1; -sProcess.options.regionfrontal.Comment = '2: Frontal'; -sProcess.options.regionfrontal.Type = 'checkbox'; -sProcess.options.regionfrontal.Value = 1; -sProcess.options.regioncentral.Comment = '3: Central'; -sProcess.options.regioncentral.Type = 'checkbox'; -sProcess.options.regioncentral.Value = 1; -sProcess.options.regionparietal.Comment = '4: Parietal'; -sProcess.options.regionparietal.Type = 'checkbox'; -sProcess.options.regionparietal.Value = 1; -sProcess.options.regiontemporal.Comment = '5: Temporal'; -sProcess.options.regiontemporal.Type = 'checkbox'; -sProcess.options.regiontemporal.Value = 1; -sProcess.options.regionoccipital.Comment = '6: Occipital'; -sProcess.options.regionoccipital.Type = 'checkbox'; -sProcess.options.regionoccipital.Value = 1; -sProcess.options.regionlimbic.Comment = '7: Limbic'; -sProcess.options.regionlimbic.Type = 'checkbox'; -sProcess.options.regionlimbic.Value = 1; +sProcess.options.region.Comment = [{'Prefrontal', 'Frontal', 'Central', 'Parietal', 'Temporal', 'Occipital', 'Limbic'}, {'Select region(s) to include:'}]; +sProcess.options.region.Type = 'list_horizontal'; +sProcess.options.region.Value = ''; % Add separator sProcess.options.separator1.Type = 'separator'; % Method for sorting the data @@ -132,13 +112,7 @@ % Color figure by region or by label OPTIONS.ColorScheme = sProcess.options.colorscheme.Value; % Select regions to include - OPTIONS.Region = logical([sProcess.options.regionprefrontal.Value - sProcess.options.regionfrontal.Value - sProcess.options.regioncentral.Value - sProcess.options.regionparietal.Value - sProcess.options.regiontemporal.Value - sProcess.options.regionoccipital.Value - sProcess.options.regionlimbic.Value]); + OPTIONS.Region = sProcess.options.region.Value; % Method for sorting the data OPTIONS.SortMethod = sProcess.options.sortmethod.Value; % Sort window From 09a461796d62c7622e21298321acd5aa552e91bd Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 19 Jun 2026 00:09:51 -0700 Subject: [PATCH 16/55] Use `radio_linelabel` type for `colorscheme` --- toolbox/process/functions/process_fastgraph.m | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index a02446f880..a20525ae6f 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -50,10 +50,9 @@ sProcess.options.scouts.Type = 'scout'; sProcess.options.scouts.Value = {}; % Color Fastgraph by region or by label -sProcess.options.label2.Comment = 'Color Fastgraph by region or by label ?'; -sProcess.options.label2.Type = 'label'; -sProcess.options.colorscheme.Comment = {'Region', 'Label'; 'Region', 'Label'}; -sProcess.options.colorscheme.Type = 'radio_label'; +sProcess.options.colorscheme.Comment = {'Region', 'Label', 'FastGraph color:'; ... + 'Region', 'Label', ''}; +sProcess.options.colorscheme.Type = 'radio_linelabel'; sProcess.options.colorscheme.Value = 'Region'; % Select regions to include sProcess.options.region.Comment = [{'Prefrontal', 'Frontal', 'Central', 'Parietal', 'Temporal', 'Occipital', 'Limbic'}, {'Select region(s) to include:'}]; From e7d60851eac72c5e42c839b81435a4619bc25cf8 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 19 Jun 2026 00:13:56 -0700 Subject: [PATCH 17/55] Merge the consecutive labels --- toolbox/process/functions/process_fastgraph.m | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index a20525ae6f..3ae4a1f574 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -67,13 +67,12 @@ sProcess.options.sortmethod.Type = 'radio_label'; sProcess.options.sortmethod.Value = 'Root Mean Square'; % Sort window -sProcess.options.label6.Comment = 'Choose range to sort over:'; -sProcess.options.label6.Type = 'label'; -sProcess.options.label7.Comment = ['' ... +sProcess.options.label6.Comment = ['Choose range to sort over:' ... + '' ... 'Early latency:    0-60 ms
' ... 'Middle latency: 60-250 ms
' ... - 'Late latency:     250-600 ms']; -sProcess.options.label7.Type = 'label'; + 'Late latency:     250-600 ms']; +sProcess.options.label6.Type = 'label'; sProcess.options.sortwindow.Comment = 'Sort range: '; sProcess.options.sortwindow.Type = 'timewindow'; sProcess.options.sortwindow.Value = []; @@ -88,9 +87,9 @@ sProcess.options.edgealpha.Type = 'value'; sProcess.options.edgealpha.Value = {0.05,' ', 2}; % Exclude contacts within a certain distance from the stimulation sites -sProcess.options.label8.Comment = ['' ... +sProcess.options.label7.Comment = ['' ... 'Exclude analysis of contacts within this distance from the stimulation site']; -sProcess.options.label8.Type = 'label'; +sProcess.options.label7.Type = 'label'; sProcess.options.excluderadius.Comment = 'Exclusion zone radius: '; sProcess.options.excluderadius.Type = 'value'; sProcess.options.excluderadius.Value = {20,'mm', 0}; From e1a0f21783ca7fd50bc0045f4d094875ec23ba50 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 19 Jun 2026 00:52:03 -0700 Subject: [PATCH 18/55] Bugfix: Handle region list --- toolbox/process/functions/process_fastgraph.m | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 3ae4a1f574..3f84b83aea 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -55,7 +55,7 @@ sProcess.options.colorscheme.Type = 'radio_linelabel'; sProcess.options.colorscheme.Value = 'Region'; % Select regions to include -sProcess.options.region.Comment = [{'Prefrontal', 'Frontal', 'Central', 'Parietal', 'Temporal', 'Occipital', 'Limbic'}, {'Select region(s) to include:'}]; +sProcess.options.region.Comment = [{'Prefrontal (PF)', 'Frontal (F)', 'Central (C)', 'Parietal (P)', 'Temporal (T)', 'Occipital (O)', 'Limbic (L)'}, {'Select region(s) to include:'}]; sProcess.options.region.Type = 'list_horizontal'; sProcess.options.region.Value = ''; % Add separator @@ -140,7 +140,7 @@ OPTIONS = GetOptions(sProcess); % Early exit if no region is selected - if ~any(OPTIONS.Region) + if isempty(OPTIONS.Region) bst_report('Error', sProcess, [], 'No region selected. Select at least one region to run the analysis.'); return; end @@ -627,8 +627,7 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL isKeep = ismember({atlas.Scouts.Label}, OPTIONS.AtlasScoutLabels); else % Region-based filtering - allRegionCodes = {'PF','F','C','P','T','O','L'}; - selectedRegions = allRegionCodes(OPTIONS.Region); + selectedRegions = regexprep(OPTIONS.Region, '^.*\((.*?)\).*$', '$1'); % Remove the leading character from Brainstorm scout region code scoutRegions = cellfun(@(x) x(2:end), {atlas.Scouts.Region}, 'UniformOutput', false); isKeep = ismember(scoutRegions, selectedRegions); From 78deddf9cb370dc715f5d14a5e9298a6b1b05148 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Fri, 19 Jun 2026 00:58:23 -0700 Subject: [PATCH 19/55] Use `FastGraph` capitalization in comments --- toolbox/process/functions/process_fastgraph.m | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 3f84b83aea..b75235a425 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -1,5 +1,5 @@ function varargout = process_fastgraph( varargin ) -% PROCESS_FASTGRAPH: Plot fastgraph for one or more SEEG recordings. +% PROCESS_FASTGRAPH: Plot FastGraph for one or more SEEG recordings. % For each stimulation pair, channels are split by hemisphere, sorted % by a user-selected metric, filtered by atlas region or scout label, and % plotted as stacked area plots @@ -35,7 +35,7 @@ %% ===== GET DESCRIPTION ===== function sProcess = GetDescription() %#ok % Describe the process and its UI options -sProcess.Comment = 'Plot Fastgraphs'; +sProcess.Comment = 'Plot FastGraphs'; sProcess.Category = 'Custom'; sProcess.SubGroup = 'FAST graph'; sProcess.Index = 1303; @@ -45,11 +45,11 @@ sProcess.OutputTypes = {'data'}; sProcess.nInputs = 1; sProcess.nMinFiles = 1; -% Scouts to use for plotting Fastgraph +% Scouts to use for plotting FastGraph sProcess.options.scouts.Comment = ''; sProcess.options.scouts.Type = 'scout'; sProcess.options.scouts.Value = {}; -% Color Fastgraph by region or by label +% Color FastGraph by region or by label sProcess.options.colorscheme.Comment = {'Region', 'Label', 'FastGraph color:'; ... 'Region', 'Label', ''}; sProcess.options.colorscheme.Type = 'radio_linelabel'; @@ -103,7 +103,7 @@ %% ===== GET OPTIONS ===== function OPTIONS = GetOptions(sProcess) OPTIONS = struct(); - % Atlas and scouts to use for plotting Fastgraph + % Atlas and scouts to use for plotting FastGraph ScoutsList = sProcess.options.scouts.Value; OPTIONS.Atlas = ScoutsList{1,1}; OPTIONS.AtlasScoutLabels = ScoutsList{1,2}; @@ -152,10 +152,10 @@ iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); % Get the midpoint location of each stimulation pair from channel stimLocs = GetStimLocs(sInputs, ChannelMat); - % Sort fastgraphs by stimulation-site location for LAPRAP style display + % Sort FastGraphs by stimulation-site location for LAPRAP style display sSortedFastgraphLocIdxs = SortLAPRAP(stimLocs); - % Load SEEG recordings after applying Fastgraph sorting + % Load SEEG recordings after applying FastGraph sorting [seegData, excludedContacts] = GetSeegData(sInputs, sSortedFastgraphLocIdxs, stimLocs, ChannelMat, OPTIONS); % Split SEEG contacts into left and right hemisphere groups sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat); @@ -169,11 +169,11 @@ chanNamesSeeg = chanTableWithAtlas(2:end, 1); atlasScoutLabelsSeeg = chanTableWithAtlas(2:end, cols); - % Create figure for Fastgraph + % Create figure for FastGraph figure; % Maximize figure set(gcf, 'Position', get(0,'Screensize')); - % Shared y-axis limits across Fastgraph subplots + % Shared y-axis limits across FastGraph subplots commonAxisLimits = []; % Reserve one extra subplot for the legend nSubplots = length(sInputs)+1; @@ -185,8 +185,8 @@ gap = [0.075 0.0175]; horzMargin = 0.03; vertMargin = 0.015; - % Generate one fastgraph per selected input - bst_progress('start', 'Process', 'Plotting Fastgraphs...', 0, 100); + % Generate one FastGraph per selected input + bst_progress('start', 'Process', 'Plotting FastGraphs...', 0, 100); for iSubplot = 1:nSubplots-1 % Show progress progressPrc = round(100 .* iSubplot ./ (nSubplots-1)); @@ -206,13 +206,13 @@ sSubplotDataSorted = ApplyDataSorting(subplotData, seegData, OPTIONS); % Create the subplot with custom spacing subtightplot(nRows, nCols, iSubplot, gap, horzMargin, vertMargin); - % Plot the Fastgraph for the current stimulation pair + % Plot the FastGraph for the current stimulation pair [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInputs, stimLocs, iSubplot, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); % Tighten axes to the plotted data and store the current axis handle axis tight axisLimits = axis; axSubplots(iSubplot) = gca; - % Update the shared y-axis limits so all Fastgraph subplots can + % Update the shared y-axis limits so all FastGraph subplots can % use the same vertical range for visual comparison if iSubplot == 1 commonAxisLimits = axisLimits; @@ -230,7 +230,7 @@ % Add the stimulation pair and atlas scout label as the subplot title AddFastgraphTitle(sInputs, sSortedFastgraphLocIdxs.All(iSubplot), chanNamesSeeg, atlasScoutLabelsSeeg); end - % Apply the shared y-axis limits to all Fastgraph subplots + % Apply the shared y-axis limits to all FastGraph subplots for iSubplot = 1:nSubplots-1 axSubplots(iSubplot).YLim = commonAxisLimits(3:4); end @@ -443,7 +443,7 @@ end %% ===== PLOT FASTGRAPH ===== -% Create one fastgraph subplot. +% Create one FastGraph subplot. % Left-hemisphere SEEG channels are plotted as positive stacked areas % Right-hemisphere SEEG channels are plotted as negative stacked areas function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInputs, stimLocs, iSubplot, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) @@ -467,7 +467,7 @@ plotWindowIdx = OPTIONS.PlotWindow(1) + 101 : OPTIONS.PlotWindow(2) + 101; timeMs = seegData{iSubplot}.Time(plotWindowIdx) * 1000; - fprintf('\n===== Fastgraph %d/%d: Stimulation site "%s" =====\n', iSubplot, numel(sInputs), sInputs(iSubplot).Comment) + fprintf('\n===== FastGraph %d/%d: Stimulation site "%s" =====\n', iSubplot, numel(sInputs), sInputs(iSubplot).Comment) % Extract SEEG data once for this subplot Fout = seegData{iSubplot}.F(iSeeg, :); @@ -682,7 +682,7 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL end %% ===== PLOT LEGEND ===== -% Shows the legend for the Fastgraph plots as in the paper +% Shows the legend for the FastGraph plots as in the paper function PlotLegend(axSubplotLegend, brainImg, xRange, yRange, xLabel, yLabel) % === Prepare the plot area === % Set the visible x- and y-axis limits From d7ac4858f87080b1b896f7ff58f21344c6294f50 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 24 Jun 2026 16:42:18 -0700 Subject: [PATCH 20/55] Tutorial: Update to `process_cutstim` --- toolbox/script/tutorial_fastgraph.m | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index 698170a413..e412493201 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -149,10 +149,12 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'baseline', []); % Process: Remove SPES artifacts -sFilesStimStartRmSpes = bst_process('CallProcess', 'process_remove_spes_artifacts', sFilesStimStart, [], ... - 'stimevent', 'STIM', ... - 'timeart', 0.005, ... % in ms - 'timespline', 0.003); % in ms +sFilesStimStartRmSpes = bst_process('CallProcess', 'process_cutstim', sFilesStimStart, [], ... + 'eventname', 'STIM', ... + 'timewindow', [0, 0.005], ... % in ms + 'sensortypes', 'SEEG', ... + 'method', 'spline', ... + 'splinebuffer', 0.003); % in ms % Process: Remove drift EMD sFilesStimStartEmd = bst_process('CallProcess', 'process_remove_drift_emd', sFilesStimStartRmSpes, [], ... From 09436065c920d28f6bceb451650bfc3efadff28c Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Wed, 24 Jun 2026 16:42:32 -0700 Subject: [PATCH 21/55] Tutorial: Update to `process_detrend_emd` --- toolbox/script/tutorial_fastgraph.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index e412493201..3aeb468877 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -157,8 +157,9 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'splinebuffer', 0.003); % in ms % Process: Remove drift EMD -sFilesStimStartEmd = bst_process('CallProcess', 'process_remove_drift_emd', sFilesStimStartRmSpes, [], ... - 'cutoff', 2); % in Hz +sFilesStimStartEmd = bst_process('CallProcess', 'process_detrend_emd', sFilesStimStartRmSpes, [], ... + 'sensortypes', 'SEEG', ... + 'emdcutoff', 2); % in Hz % Process: Load the STIM events sFilesStim = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartEmd, [], ... From 8e4fe5f4902ad441b995aaca7ef9351b5ca349a4 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 16 Jul 2026 11:54:12 -0700 Subject: [PATCH 22/55] Update header comments --- toolbox/io/export_channel_atlas.m | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/toolbox/io/export_channel_atlas.m b/toolbox/io/export_channel_atlas.m index 8dac6eb169..d2a759f89f 100644 --- a/toolbox/io/export_channel_atlas.m +++ b/toolbox/io/export_channel_atlas.m @@ -1,8 +1,8 @@ function [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, Modality, TsvFile, Radius, isProba, isInteractive) % EXPORT_CHANNEL_ATLAS: Compute anatomical labels for SEEG/ECOG contacts from volume and surface parcellations % -% USAGE: TsvFile = export_channel_atlas(ChannelFile, Modality='ECOG+SEEG', TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1) -% TsvFile = export_channel_atlas(ChannelFile, iChannels, TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1) +% USAGE: [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, Modality='ECOG+SEEG', TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1) +% [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, iChannels, TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1) % % INPUT: % - ChannelFile : Path to Brainstorm channel file to be processed @@ -14,6 +14,13 @@ % : If 0, use all available Coodinates, Parcellations (anat) and Atlases (surface), % and do not display output table % - iChannels : Limit export to a subset of channel indices +% OUTPUT: +% - TsvFile : Output text file (tab-separated values). Empty when no file was selected or requested. +% - ChanTable : Cell array containing the complete output table. The same information is written to TsvFile +% when an output file is requested. +% - The first column always contains channel names +% - The remaining columns contain available coordinates, anatomical labels, +% and optional probabilities corresponding to each channel % % REFERENCES: % - MERCIER M, 2021: From 145e39630bd512ef637dbe35b817586121e4436f Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 16 Jul 2026 14:04:18 -0700 Subject: [PATCH 23/55] Validate all input files share the same channel file --- toolbox/process/functions/process_fastgraph.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index b75235a425..12922b6bf2 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -136,6 +136,12 @@ % Initialize output OutputFiles = {}; + % Check that all input files use the same channel file + ChannelFiles = {sInputs.ChannelFile}; + if length(unique(ChannelFiles)) > 1 + bst_report('Error', sProcess, sInputs, 'All input files must use the same channel file.'); + return; + end % Get options OPTIONS = GetOptions(sProcess); From 59ab0317576c2e8ffff01397885ac5f7a4ce90e1 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 16 Jul 2026 14:11:47 -0700 Subject: [PATCH 24/55] Validate that all comments have valid bipolar channel names --- toolbox/process/functions/process_fastgraph.m | 36 ++++++++++++++++--- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 12922b6bf2..90601cf232 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -142,6 +142,37 @@ bst_report('Error', sProcess, sInputs, 'All input files must use the same channel file.'); return; end + + % ===== Check that every comment contains a bipolar channel name ===== + % Extract bipolar channel pairs from all comments + bipolarPattern = '([A-Za-z]+''?\d+)\s*-\s*([A-Za-z]+''?\d+)'; + bipolarChannels = regexp({sInputs.Comment}, bipolarPattern, 'tokens', 'once'); + % Check that every comment contains a bipolar pair + isBipolar = ~cellfun(@isempty, bipolarChannels); + if ~all(isBipolar) + iInvalid = find(~isBipolar); + bst_report('Error', sProcess, sInputs(iInvalid), ... + sprintf('Could not find a bipolar channel name in the file comment: "%s".\n', ... + sInputs(iInvalid).Comment)); + return; + end + + % ===== Check whether all channel names in comment are valid ===== + % Load the channel file + ChannelMat = in_bst_channel(ChannelFiles{1}); + channelNames = {ChannelMat.Channel.Name}; + % Flatten all extracted pairs + allBipolarChannels = [bipolarChannels{:}]; + % Check whether all extracted channel names exist + isChannelFound = ismember(allBipolarChannels, channelNames); + if ~all(isChannelFound) + missingChannels = unique(allBipolarChannels(~isChannelFound), 'stable'); + bst_report('Error', sProcess, sInputs, ... + sprintf('The following channels were not found in the channel file: %s.', ... + strjoin(missingChannels, ', '))); + return; + end + % Get options OPTIONS = GetOptions(sProcess); @@ -151,9 +182,6 @@ return; end - % Load the channel file - ChannelFile = file_fullpath(sInputs(1).ChannelFile); - ChannelMat = load(ChannelFile); % Get indices of SEEG channels iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); % Get the midpoint location of each stimulation pair from channel @@ -166,7 +194,7 @@ % Split SEEG contacts into left and right hemisphere groups sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat); % Compute anatomical labels for the contacts from volume/surface parcellations - [~, chanTableWithAtlas] = export_channel_atlas(ChannelFile, 'SEEG', [], 10, 0, 0); + [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 10, 0, 0); % Locate atlas related columns from channel table above hit = cellfun(@(x) ischar(x) && (~isempty(strfind(OPTIONS.Atlas, x)) || ~isempty(strfind(x, OPTIONS.Atlas))), chanTableWithAtlas(1,:)); % Columns whose header matches the atlas name From 3056d8d97023eb5e6f8f8e23f32d62d8e4433ab9 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 16 Jul 2026 18:21:21 -0700 Subject: [PATCH 25/55] `5mm` radius for bipolar contacts --- toolbox/process/functions/process_fastgraph.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 90601cf232..8c90424a9a 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -194,7 +194,7 @@ % Split SEEG contacts into left and right hemisphere groups sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat); % Compute anatomical labels for the contacts from volume/surface parcellations - [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 10, 0, 0); + [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 5, 0, 0); % Locate atlas related columns from channel table above hit = cellfun(@(x) ischar(x) && (~isempty(strfind(OPTIONS.Atlas, x)) || ~isempty(strfind(x, OPTIONS.Atlas))), chanTableWithAtlas(1,:)); % Columns whose header matches the atlas name From 39bf0dcbe3066bbcc4f15bdd0cfee687eb99f605 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Thu, 16 Jul 2026 18:37:08 -0700 Subject: [PATCH 26/55] FG subplot grid: Keep it simple --- toolbox/process/functions/process_fastgraph.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 8c90424a9a..a7baecad65 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -213,8 +213,8 @@ nSubplots = length(sInputs)+1; % Define the plot parameters % Subplot grid dimensions - nRows = floor(sqrt(nSubplots/1.5)); - nCols = ceil(nSubplots/floor(sqrt(nSubplots/1.5))); + nCols = ceil(sqrt(nSubplots)); + nRows = ceil(nSubplots / nCols); % Subplot spacing and margins gap = [0.075 0.0175]; horzMargin = 0.03; From 3ed08281bfc722e24d161b146b98e44cc615e815 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Tue, 21 Jul 2026 11:11:36 -0700 Subject: [PATCH 27/55] Clean --- toolbox/process/functions/process_fastgraph.m | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index a7baecad65..55888ec850 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -152,8 +152,7 @@ if ~all(isBipolar) iInvalid = find(~isBipolar); bst_report('Error', sProcess, sInputs(iInvalid), ... - sprintf('Could not find a bipolar channel name in the file comment: "%s".\n', ... - sInputs(iInvalid).Comment)); + sprintf('Could not find a bipolar channel name in the file comment: "%s".\n', sInputs(iInvalid).Comment)); return; end @@ -168,8 +167,7 @@ if ~all(isChannelFound) missingChannels = unique(allBipolarChannels(~isChannelFound), 'stable'); bst_report('Error', sProcess, sInputs, ... - sprintf('The following channels were not found in the channel file: %s.', ... - strjoin(missingChannels, ', '))); + sprintf('The following channels were not found in the channel file: %s.', strjoin(missingChannels, ', '))); return; end From c7a7dfe8463a946b8329e7708c8ff71860d1f9b6 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Tue, 21 Jul 2026 19:59:51 -0700 Subject: [PATCH 28/55] Refine FastGraph axes and legend rendering --- toolbox/process/functions/process_fastgraph.m | 166 ++++++++---------- 1 file changed, 73 insertions(+), 93 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 55888ec850..41a2076474 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -262,12 +262,19 @@ % Add the stimulation pair and atlas scout label as the subplot title AddFastgraphTitle(sInputs, sSortedFastgraphLocIdxs.All(iSubplot), chanNamesSeeg, atlasScoutLabelsSeeg); end - % Apply the shared y-axis limits to all FastGraph subplots - for iSubplot = 1:nSubplots-1 - axSubplots(iSubplot).YLim = commonAxisLimits(3:4); - end + % Format all FastGraph axes + axFastGraphs = axSubplots(1:nSubplots-1); + % Set axis labels for all FastGraph subplots + set([axFastGraphs.XLabel], 'String', 'Time (ms)'); + set([axFastGraphs.YLabel], 'String', 'Voltage (mV)'); + % Apply common axes properties + set(axFastGraphs, ... + 'YLim', commonAxisLimits(3:4), ... + 'XAxisLocation', 'bottom'); + % Add zero-reference lines and hemisphere labels + DecorateFastgraphAxes(axFastGraphs); % Link subplot axes so that zooming stays synchronized - linkaxes(axSubplots) + linkaxes(axFastGraphs) set(gcf,'units','normalized','outerposition',[0 0 1 1]) zoom on @@ -279,7 +286,7 @@ subtightplot(nRows, nCols, iSubplot+1, gap, horzMargin, vertMargin); % Plot the reference panel with the cortex snapshot and axis labels axSubplots(iSubplot+1) = gca; - PlotLegend(axSubplots(iSubplot+1), imgCortex, round(axSubplots(1).XLim), [0 1], 'Time (ms)', 'Voltage (mV)'); + PlotLegend(axSubplots(iSubplot+1), imgCortex, round(axSubplots(1).XLim), axSubplots(1).YLim); % Close progress bst_progress('stop'); @@ -714,109 +721,82 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL end %% ===== PLOT LEGEND ===== -% Shows the legend for the FastGraph plots as in the paper -function PlotLegend(axSubplotLegend, brainImg, xRange, yRange, xLabel, yLabel) - % === Prepare the plot area === - % Set the visible x- and y-axis limits - set(axSubplotLegend, 'XLim', xRange, 'YLim', yRange); - % Add x-axis label - axSubplotLegend.XLabel.String = xLabel; - % Move x-axis label closer to the axis (slightly upward) - axSubplotLegend.XLabel.Position = [mean(axSubplotLegend.XLim), axSubplotLegend.YLim(1) - 0.01, 0]; - % Add y-axis label - axSubplotLegend.YLabel.String = yLabel; - % Move the y-axis label closer to the axis (slightly right) - axSubplotLegend.YLabel.Position = [axSubplotLegend.XLim(1) - 5, mean(axSubplotLegend.YLim), 0]; - % Show ticks only at the minimum and maximum values of each axis - axSubplotLegend.XTick = [xRange(1), xRange(2)]; - axSubplotLegend.YTick = [yRange(1), yRange(2)]; - - % === Create overlay axes for the brain atlas image === - axImg = axes('Parent', ancestor(axSubplotLegend, 'figure'), ... +% Show the reference cortex image using the same axes layout as the plots +function PlotLegend(axLegend, brainImg, xLim, yLim) + % Configure legend axes + set(axLegend, ... + 'XLim', xLim, ... + 'YLim', yLim, ... + 'XAxisLocation', 'bottom'); + axLegend.XLabel.String = 'Time (ms)'; + axLegend.YLabel.String = 'Voltage (mV)'; + % Create overlay axes for the cortex image + hFig = ancestor(axLegend, 'figure'); + axImg = axes( ... + 'Parent', hFig, ... 'Units', 'pixels', ... 'Color', 'none'); - % Display the brain image inside the overlay axes - hImg = imshow(brainImg, 'Parent', axImg); - % Hide the overlay axes so only the image is visible + imshow(brainImg, 'Parent', axImg); axis(axImg, 'off'); - % Keep the original axes limits fixed so the image does not alter them - axis(axSubplotLegend, 'manual'); - % Initial placement - UpdateLegendImage(axSubplotLegend, axImg, brainImg); - % Update placement whenever the figure is resized/moved - hFig = ancestor(axSubplotLegend, 'figure'); - hFig.SizeChangedFcn = @(~,~)UpdateLegendImage(axSubplotLegend, axImg, brainImg); - - % Add left/right hemisphere labels with pixel-based spacing - AddLegendHemisphereLabels(axSubplotLegend, xRange, yRange); + axis(axLegend, 'manual'); + % Position the image initially and after resizing + UpdateLegendImage(axLegend, axImg, brainImg); + hFig.SizeChangedFcn = @(~,~) UpdateLegendImage(axLegend, axImg, brainImg); end -%% ===== ADD 'L/R' HEMISPHERE LABELS IN THE LEGEND ===== -% Add 'L/R' hemisphere labels to the legend axes -function AddLegendHemisphereLabels(axSubplotLegend, xRange, yRange) - % Position labels near the right side of the legend axes - xSpan = diff(xRange); - ySpan = diff(yRange); - xLR = xRange(2) - 0.08 * xSpan; - - % Get axes height in pixels - oldUnits = axSubplotLegend.Units; - axSubplotLegend.Units = 'pixels'; - axPos = axSubplotLegend.Position; - axSubplotLegend.Units = oldUnits; - - % Convert a fixed pixel gap into data units - pixelsPerDataY = axPos(4) / ySpan; - gapPx = max(14, axSubplotLegend.FontSize + 4); - gapData = gapPx / pixelsPerDataY; - - % Place labels above and below the x-axis - yAxisLevel = yRange(1); - yL = yAxisLevel + gapData; - yR = yAxisLevel - gapData; - - % Draw the labels - text(axSubplotLegend, xLR, yL, 'L', ... - 'FontSize', 8, ... - 'FontWeight', 'bold', ... - 'HorizontalAlignment', 'right', ... - 'VerticalAlignment', 'middle', ... - 'Clipping', 'off', ... - 'Margin', 1); +%% ===== DECORATE FASTGRAPH AXES ===== +% Add the zero-reference line and L/R hemisphere labels +function DecorateFastgraphAxes(axFastgraphs) + for ax = axFastgraphs + line(ax, ax.XLim, [0 0], ... + 'Color', [0 0 0], ... + 'LineWidth', 0.5, ... + 'HandleVisibility', 'off'); + AddHemisphereLabels(ax); + end +end - text(axSubplotLegend, xLR, yR, 'R', ... +%% ===== ADD HEMISPHERE LABELS ===== +% Add L/R labels immediately above and below the y = 0 reference line +function AddHemisphereLabels(ax) + % Common label properties + labelProperties = { ... + 'Parent', ax, ... + 'Units', 'normalized', ... 'FontSize', 8, ... 'FontWeight', 'bold', ... 'HorizontalAlignment', 'right', ... 'VerticalAlignment', 'middle', ... - 'Clipping', 'off', ... - 'Margin', 1); + 'Clipping', 'off'}; + % Find the normalized vertical position of y = 0 + yZeroNormalized = -ax.YLim(1) / diff(ax.YLim); + % Position labels close to the right edge + xPosition = 0.96; + % Normalized vertical spacing around the zero line + labelGap = 0.035; + % Place labels above and below the bottom x-axis + text(xPosition, yZeroNormalized + labelGap, 'L', labelProperties{:}); + text(xPosition, yZeroNormalized - labelGap, 'R', labelProperties{:}); end %% ===== UPDATE LEGEND IMAGE ===== % Update the overlay image position so it stays centered inside the % legend subplot when the figure is resized or moved across screens function UpdateLegendImage(axSubplotLegend, axImg, brainImg) - % Get original image size in pixels - imgH = size(brainImg, 1); - imgW = size(brainImg, 2); - % Read the legend subplot position in pixel units - oldUnits = axSubplotLegend.Units; - axSubplotLegend.Units = 'pixels'; - % Get the axes position in pixel units: [left, bottom, width, height] - pos = axSubplotLegend.Position; - axSubplotLegend.Units = oldUnits; - % Available subplot width and height in pixels - boxW = pos(3); - boxH = pos(4); - % Scale the image to fit inside the subplot while preserving aspect ratio - scale = min(boxW / imgW, boxH / imgH) * 0.75; + % Axes position in pixels + axPos = getpixelposition(axSubplotLegend); + % Original image dimensions + imgSize = size(brainImg); + imgH = imgSize(1); + imgW = imgSize(2); + % Scale while preserving image aspect ratio + scale = 0.75 * min(axPos(3) / imgW, axPos(4) / imgH); newW = imgW * scale; newH = imgH * scale; - % Center the image inside the legend subplot - xLeft = pos(1) + (boxW - newW) / 2; - yBottom = pos(2) + (boxH - newH) / 2; - % Update the overlay axes position in pixel coordinates - axImg.Units = 'pixels'; - axImg.Position = [xLeft, yBottom, newW, newH]; + % Center inside the legend subplot + xLeft = axPos(1) + (axPos(3) - newW) / 2; + yBottom = axPos(2) + (axPos(4) - newH) / 2; + set(axImg, ... + 'Units', 'pixels', ... + 'Position', [xLeft, yBottom, newW, newH]); end \ No newline at end of file From 7171dd55e34a3764fef67c39df1b00045c2a467a Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Tue, 21 Jul 2026 20:10:49 -0700 Subject: [PATCH 29/55] Simplify cortex snapshot capture --- toolbox/process/functions/process_fastgraph.m | 24 +++++++------------ 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 41a2076474..2f31aa9f95 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -697,25 +697,17 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL bst_figures('SetBackgroundColor', hFigSurf, [1 1 1]); % Select atlas panel_scout('SetCurrentAtlas', iAtlas); - % Set options - switch(OPTIONS.ColorScheme) - case 'Region' - panel_scout('SetScoutsOptions', 0, 0, 1, 'select', 0, 1, 0, 1); - case 'Label' - panel_scout('SetScoutsOptions', 0, 0, 1, 'select', 0, 1, 0, 0); - end + % Color scouts by region or individual label + isRegionColor = strcmp(OPTIONS.ColorScheme, 'Region'); + panel_scout('SetScoutsOptions', 0, 0, 1, 'select', 0, 1, 0, isRegionColor); % Show only selected scouts panel_scout('SetSelectedScouts', iSelectedScouts); - % Capture image + % Set background color + bst_figures('SetBackgroundColor', hFigSurf, [1 1 1]); + % Capture and crop the cortex image img = out_figure_image(hFigSurf); - % Crop background - bgColor = img(1,1,:); - mask = (img(:,:,1) == bgColor(1)) & ... - (img(:,:,2) == bgColor(2)) & ... - (img(:,:,3) == bgColor(3)); - goodRows = any(~mask, 2); - goodCols = any(~mask, 1); - imgCortex = img(goodRows, goodCols, :); + isBackground = all(img == 255, 3); + imgCortex = img(any(~isBackground, 2), any(~isBackground, 1), :); % Close figure close(hFigSurf); end From 2eb5129598930f453319e710658cc876ef1c76c2 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Tue, 21 Jul 2026 23:14:23 -0700 Subject: [PATCH 30/55] Add `Plot FastGraphs` to context menu --- toolbox/tree/tree_callbacks.m | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/toolbox/tree/tree_callbacks.m b/toolbox/tree/tree_callbacks.m index bfda27116a..91ed053460 100644 --- a/toolbox/tree/tree_callbacks.m +++ b/toolbox/tree/tree_callbacks.m @@ -1681,6 +1681,11 @@ fcnPopupProjectSources(0); fcnPopupScoutTimeSeries(jPopup); end + + % Plot FastGraphs + if ~isempty(AllMod) && ismember('SEEG', AllMod) && strcmpi(DataType, 'recordings') + gui_component('MenuItem', jPopup, [], 'Plot FastGraphs', IconLoader.ICON_TS_DISPLAY, [], @(h,ev)panel_process_select('ShowPanelForFile', GetAllFilenames(bstNodes, 'data'), 'process_fastgraph')); + end %% ===== POPUP: STAT/DATA ===== case 'pdata' From ca0cdaa16cd16af01f158cb249494158ccaeef48 Mon Sep 17 00:00:00 2001 From: Chinmay Chinara Date: Tue, 21 Jul 2026 23:47:43 -0700 Subject: [PATCH 31/55] Clean --- toolbox/process/functions/process_fastgraph.m | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 2f31aa9f95..0140e3f32d 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -268,9 +268,7 @@ set([axFastGraphs.XLabel], 'String', 'Time (ms)'); set([axFastGraphs.YLabel], 'String', 'Voltage (mV)'); % Apply common axes properties - set(axFastGraphs, ... - 'YLim', commonAxisLimits(3:4), ... - 'XAxisLocation', 'bottom'); + set(axFastGraphs, 'YLim', commonAxisLimits(3:4), 'XAxisLocation', 'bottom'); % Add zero-reference lines and hemisphere labels DecorateFastgraphAxes(axFastGraphs); % Link subplot axes so that zooming stays synchronized @@ -724,10 +722,7 @@ function PlotLegend(axLegend, brainImg, xLim, yLim) axLegend.YLabel.String = 'Voltage (mV)'; % Create overlay axes for the cortex image hFig = ancestor(axLegend, 'figure'); - axImg = axes( ... - 'Parent', hFig, ... - 'Units', 'pixels', ... - 'Color', 'none'); + axImg = axes('Parent', hFig, 'Units', 'pixels', 'Color', 'none'); imshow(brainImg, 'Parent', axImg); axis(axImg, 'off'); axis(axLegend, 'manual'); @@ -788,7 +783,5 @@ function UpdateLegendImage(axSubplotLegend, axImg, brainImg) % Center inside the legend subplot xLeft = axPos(1) + (axPos(3) - newW) / 2; yBottom = axPos(2) + (axPos(4) - newH) / 2; - set(axImg, ... - 'Units', 'pixels', ... - 'Position', [xLeft, yBottom, newW, newH]); + set(axImg, 'Units', 'pixels', 'Position', [xLeft, yBottom, newW, newH]); end \ No newline at end of file From 787ff2891106210fa66a5a2866b5bd42703ccf15 Mon Sep 17 00:00:00 2001 From: rcassani Date: Wed, 22 Jul 2026 17:13:39 -0400 Subject: [PATCH 32/55] Update GUI in `process_fastgraph.m` --- toolbox/process/functions/process_fastgraph.m | 79 ++++++++++--------- 1 file changed, 40 insertions(+), 39 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 0140e3f32d..7af0710723 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -37,7 +37,7 @@ % Describe the process and its UI options sProcess.Comment = 'Plot FastGraphs'; sProcess.Category = 'Custom'; -sProcess.SubGroup = 'FAST graph'; +sProcess.SubGroup = 'FastGraph'; sProcess.Index = 1303; sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/FastGraph'; % Definition of the input accepted by this process @@ -49,50 +49,50 @@ sProcess.options.scouts.Comment = ''; sProcess.options.scouts.Type = 'scout'; sProcess.options.scouts.Value = {}; -% Color FastGraph by region or by label -sProcess.options.colorscheme.Comment = {'Region', 'Label', 'FastGraph color:'; ... - 'Region', 'Label', ''}; -sProcess.options.colorscheme.Type = 'radio_linelabel'; -sProcess.options.colorscheme.Value = 'Region'; +% Color FastGraph by Region or by Scout +sProcess.options.colorscheme.Comment = {'Region', 'Scout', 'Color scheme:  '; ... + 'region', 'scout', ''}; +sProcess.options.colorscheme.Type = 'radio_linelabel'; +sProcess.options.colorscheme.Value = 'region'; +sProcess.options.colorscheme.Controller = struct('region', 'region'); % Select regions to include -sProcess.options.region.Comment = [{'Prefrontal (PF)', 'Frontal (F)', 'Central (C)', 'Parietal (P)', 'Temporal (T)', 'Occipital (O)', 'Limbic (L)'}, {'Select region(s) to include:'}]; +regionsStr = {'Prefrontal (PF)', 'Frontal (F)', 'Central (C)', 'Parietal (P)', 'Temporal (T)', 'Occipital (O)', 'Limbic (L)'}; +sProcess.options.region.Comment = [regionsStr, {'Select regions to include:'}]; sProcess.options.region.Type = 'list_horizontal'; -sProcess.options.region.Value = ''; -% Add separator -sProcess.options.separator1.Type = 'separator'; +sProcess.options.region.Value = regionsStr; +sProcess.options.region.Class = 'region'; % Method for sorting the data -sProcess.options.label5.Comment = 'Select method to sort the data:'; -sProcess.options.label5.Type = 'label'; -sProcess.options.sortmethod.Comment = {'Root Mean Square', 'Max Absolute'; 'Root Mean Square', 'Max Absolute'}; +sProcess.options.label1.Comment = 'Method for sorting data:'; +sProcess.options.label1.Type = 'label'; +sProcess.options.sortmethod.Comment = {'Root Mean Square', 'Max Absolute'; 'rms', 'maxabs'}; sProcess.options.sortmethod.Type = 'radio_label'; -sProcess.options.sortmethod.Value = 'Root Mean Square'; +sProcess.options.sortmethod.Value = 'rms'; % Sort window -sProcess.options.label6.Comment = ['Choose range to sort over:' ... - '' ... - 'Early latency:    0-60 ms
' ... - 'Middle latency: 60-250 ms
' ... - 'Late latency:     250-600 ms']; -sProcess.options.label6.Type = 'label'; -sProcess.options.sortwindow.Comment = 'Sort range: '; +sProcess.options.sortwindow.Comment = 'Time window to sort data: '; sProcess.options.sortwindow.Type = 'timewindow'; sProcess.options.sortwindow.Value = []; -% Add separator -sProcess.options.separator2.Type = 'separator'; +sProcess.options.label2.Comment = ['' ... + 'E.g., Early latency (0-60 ms), ' ... + 'Middle latency (60-250 ms), or ' ... + 'Late latency (250-600 ms)']; +sProcess.options.label2.Type = 'label'; +% Exclude contacts within a certain distance from the stimulation sites +sProcess.options.excluderadius.Comment = 'Exclusion zone radius:
'; +sProcess.options.excluderadius.Type = 'value'; +sProcess.options.excluderadius.Value = {20,'mm', 0}; +sProcess.options.label3.Comment = ['' ... + 'Exclude analysis of contacts within ' ... + 'this distance from the stimulation site']; +sProcess.options.label3.Type = 'label'; +sProcess.options.separator1.Type = 'separator'; % Plot window -sProcess.options.plotwindow.Comment = 'Plot range: '; +sProcess.options.plotwindow.Comment = 'Plot time range: '; sProcess.options.plotwindow.Type = 'timewindow'; sProcess.options.plotwindow.Value = []; % Edge transparency of plot sProcess.options.edgealpha.Comment = 'Edge transparency of plot: '; sProcess.options.edgealpha.Type = 'value'; sProcess.options.edgealpha.Value = {0.05,' ', 2}; -% Exclude contacts within a certain distance from the stimulation sites -sProcess.options.label7.Comment = ['' ... - 'Exclude analysis of contacts within this distance from the stimulation site']; -sProcess.options.label7.Type = 'label'; -sProcess.options.excluderadius.Comment = 'Exclusion zone radius: '; -sProcess.options.excluderadius.Type = 'value'; -sProcess.options.excluderadius.Value = {20,'mm', 0}; end %% ===== FORMAT COMMENT ===== @@ -453,8 +453,8 @@ sortWindowIdx = OPTIONS.SortWindow(1):OPTIONS.SortWindow(2); end % Sort channels within each hemisphere using the selected metric - switch OPTIONS.SortMethod - case 'Root Mean Square' + switch lower(OPTIONS.SortMethod) + case 'rms' if ~isempty(subplotData.leftData) leftDataRms = sqrt(sum(subplotData.leftData(:,sortWindowIdx).^2, 2)); leftDataRms(isnan(leftDataRms)) = -Inf; @@ -465,7 +465,7 @@ rightDataRms(isnan(rightDataRms)) = -Inf; [sSorted.Vals.Right, sSorted.Idxs.Right] = sort(rightDataRms,'ascend'); end - case 'Max Absolute' + case 'maxabs' if ~isempty(subplotData.leftData) leftDataMax = max(abs(subplotData.leftData(:,sortWindowIdx)),[],2); leftDataMax(isnan(leftDataMax)) = -Inf; @@ -613,10 +613,11 @@ % Matching scout found: assign region name region.Name = atlas.Scouts(iScout).Region(2:end); % Assign color based on the selected color scheme - if strcmp(OPTIONS.ColorScheme, 'Region') - region.Color = panel_scout('GetRegionColor', atlas.Scouts(iScout).Region); - else - region.Color = atlas.Scouts(iScout).Color; + switch lower(OPTIONS.ColorScheme) + case 'region' + region.Color = panel_scout('GetRegionColor', atlas.Scouts(iScout).Region); + case 'scout' + region.Color = atlas.Scouts(iScout).Color; end return; end @@ -696,7 +697,7 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL % Select atlas panel_scout('SetCurrentAtlas', iAtlas); % Color scouts by region or individual label - isRegionColor = strcmp(OPTIONS.ColorScheme, 'Region'); + isRegionColor = strcmpi(OPTIONS.ColorScheme, 'region'); panel_scout('SetScoutsOptions', 0, 0, 1, 'select', 0, 1, 0, isRegionColor); % Show only selected scouts panel_scout('SetSelectedScouts', iSelectedScouts); From 90ec8a7380808899d500760bf4b95d8194f166c6 Mon Sep 17 00:00:00 2001 From: rcassani Date: Wed, 22 Jul 2026 17:15:34 -0400 Subject: [PATCH 33/55] Update `tutorial_fastgraph.m` - Cleanup - Overwrite files when removing SPES artifcts and drift - Better organize imported files - Update call to `process_fastgraph` --- toolbox/script/tutorial_fastgraph.m | 67 +++++++++++++++-------------- 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/toolbox/script/tutorial_fastgraph.m b/toolbox/script/tutorial_fastgraph.m index 3aeb468877..394e022215 100644 --- a/toolbox/script/tutorial_fastgraph.m +++ b/toolbox/script/tutorial_fastgraph.m @@ -29,6 +29,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) % Authors: Chinmay Chinara, 2026 % John C. Mosher, 2026 + %% ===== PARSE INPUTS ===== % Output folder for reports if (nargin < 2) || isempty(reports_dir) || ~isfolder(reports_dir) @@ -41,19 +42,21 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) % Subject name SubjectName = 'Subject01'; + %% ===== FILES TO IMPORT ===== % Build the path of the files to import tutorial_dir = bst_fullfile(tutorial_dir, 'tutorial_fastgraph'); -MriFilePre = bst_fullfile(tutorial_dir, 'anatomy', 'pre_T1.nii.gz'); -MriCat12Path = fullfile(tutorial_dir, 'anatomy', 'cat12'); +MriFilePre = bst_fullfile(tutorial_dir, 'anatomy', 'pre_T1.nii.gz'); +MriCat12Path = bst_fullfile(tutorial_dir, 'anatomy', 'cat12'); BaselineFile = bst_fullfile(tutorial_dir, 'recordings', 'Baseline.edf'); ElecPosFile = bst_fullfile(tutorial_dir, 'recordings', 'Subject01_electrodes_mm.tsv'); % Check if the folder contains the required files -if ~file_exist(BaselineFile) +if ~file_exist(MriFilePre) || ~file_exist(BaselineFile) || ~file_exist(ElecPosFile) error(['The folder ' tutorial_dir ' does not contain the folder from the file tutorial_fastgraph.zip.']); end isMriSegmented = file_exist(bst_fullfile(MriCat12Path, 'Subject01.nii')); + %% ===== CREATE PROTOCOL ===== % The protocol name has to be a valid folder name (no spaces, no weird characters...) ProtocolName = 'TutorialFastgraph'; @@ -68,6 +71,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) % Start a new report bst_report('Start'); + %% ===== IMPORT MRI AND CT VOLUMES ===== if ~isMriSegmented % Process: Import MRI @@ -103,6 +107,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) % Reference MRI DbMriFilePre = sSubject.Anatomy(sSubject.iAnatomy).FileName; + %% ===== CREATE SEEG CONTACT IMPLANTATION ===== iStudyImplantation = db_add_condition(SubjectName, 'Implantation'); % Import locations and convert to subject coordinate system (SCS) @@ -112,6 +117,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) bst_report('Snapshot', hFigMri3d, ImplantationChannelFile, 'SEEG electrodes in 3D MRI slices'); close(hFigMri3d); + %% ===== ACCESS THE RECORDINGS ===== % Process: Create link to raw file sFileRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ... @@ -126,7 +132,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'vox2ras', 0); % Do not use the voxel=>subject transformation, already in SCS % Process: Customize SPES -bst_process('CallProcess', 'process_evt_detect_spes', sFileRaw, [], ... +sFileRaw = bst_process('CallProcess', 'process_evt_detect_spes', sFileRaw, [], ... 'stimstartlabel', 'SB', ... 'stimstoplabel', 'SE', ... 'stimchan', 'DC10', ... @@ -138,7 +144,7 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) % Process: Load the Stim Start blocks sFilesStimStart = bst_process('CallProcess', 'process_import_data_event', sFileRaw, [], ... 'subjectname', SubjectName, ... - 'condition', '', ... + 'condition', 'Epochs SB', ... 'eventname', 'SB', ... 'epochtime', [-2 32], ... % in s 'createcond', 0, ... @@ -149,22 +155,24 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'baseline', []); % Process: Remove SPES artifacts -sFilesStimStartRmSpes = bst_process('CallProcess', 'process_cutstim', sFilesStimStart, [], ... - 'eventname', 'STIM', ... - 'timewindow', [0, 0.005], ... % in ms - 'sensortypes', 'SEEG', ... - 'method', 'spline', ... - 'splinebuffer', 0.003); % in ms +sFilesStimStart = bst_process('CallProcess', 'process_cutstim', sFilesStimStart, [], ... + 'eventname', 'STIM', ... + 'timewindow', [0, 0.005], ... % in ms + 'sensortypes', 'SEEG', ... + 'method', 'spline', ... + 'splinebuffer', 0.003, ... % in ms + 'overwrite', 1); % Process: Remove drift EMD -sFilesStimStartEmd = bst_process('CallProcess', 'process_detrend_emd', sFilesStimStartRmSpes, [], ... - 'sensortypes', 'SEEG', ... - 'emdcutoff', 2); % in Hz +sFilesStimStart = bst_process('CallProcess', 'process_detrend_emd', sFilesStimStart, [], ... + 'sensortypes', 'SEEG', ... + 'emdcutoff', 2, ... ; % in Hz + 'overwrite', 1); % Process: Load the STIM events -sFilesStim = bst_process('CallProcess', 'process_import_data_event', sFilesStimStartEmd, [], ... +sFilesStim = bst_process('CallProcess', 'process_import_data_event', sFilesStimStart, [], ... 'subjectname', SubjectName, ... - 'condition', '', ... + 'condition', 'Epochs STIM', ... 'eventname', 'STIM', ... 'timewindow', [-2 32], ... % in s 'epochtime', [-0.100 0.900], ... % in ms @@ -182,22 +190,17 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) 'weighted', 0, ... 'keepevents', 0); -% Process: Plot Fastgraph +% Process: Plot FastGraphs bst_process('CallProcess', 'process_fastgraph', sFilesAvg, [], ... - 'scouts', {'Desikan-Killiany', {}}, ... - 'colorscheme', 'Region', ... % Color figures by region - 'regionprefrontal', 1, ... - 'regionfrontal', 1, ... - 'regioncentral', 1, ... - 'regionparietal', 1, ... - 'regiontemporal', 1, ... - 'regionoccipital', 1, ... - 'regionlimbic', 1, ... - 'sortmethod', 1, ... % "Root Mean Square" to sort data - 'sortwindow', [0.060, 0.250], ... % Range (middle latency) to sort the data (in ms) - 'plotwindow', [-0.100, 0.900], ... % Plot window (in ms) - 'edgealpha', 0.05, ... % Edge transparency of plot - 'excluderadius', 20); % Exclusion zone radius + 'scouts', {'Desikan-Killiany', {}}, ... + 'colorscheme', 'region', ... % Region + 'region', {'Prefrontal (PF)', 'Frontal (F)', 'Central (C)', 'Parietal (P)', 'Temporal (T)', 'Occipital (O)', 'Limbic (L)'}, ... + 'sortmethod', 'rms', ... % Root Mean Square + 'sortwindow', [0.060, 0.250], ... % Range (middle latency) to sort the data (in ms) + 'excluderadius', 20, ... % Exclusion zone radius (in mm) + 'plotwindow', [-0.100, 0.900], ... % Plot window (in ms) + 'edgealpha', 0.05); % Plot, edge transparency + %% ===== SAVE AND DISPLAY REPORT ===== ReportFile = bst_report('Save', []); @@ -207,4 +210,4 @@ function tutorial_fastgraph(tutorial_dir, reports_dir) bst_report('Open', ReportFile); end -disp([10 'DEMO> Fastgraph tutorial completed' 10]); \ No newline at end of file +disp([10 'BST> FastGraph tutorial completed' 10]); \ No newline at end of file From 162d8a333e9f15c8b6920716d0f6195fcf7a9510 Mon Sep 17 00:00:00 2001 From: rcassani Date: Mon, 27 Jul 2026 11:37:31 -0400 Subject: [PATCH 34/55] GUI: Improve labels --- toolbox/process/functions/process_fastgraph.m | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 7af0710723..3aa91f6ed2 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -72,17 +72,18 @@ sProcess.options.sortwindow.Type = 'timewindow'; sProcess.options.sortwindow.Value = []; sProcess.options.label2.Comment = ['' ... - 'E.g., Early latency (0-60 ms), ' ... - 'Middle latency (60-250 ms), or ' ... - 'Late latency (250-600 ms)']; + 'Examples:
'... + 'Early latency: 0-60 ms,
' ... + 'Middle latency: 60-250 ms, or
' ... + 'Late latency: 250-600 ms']; sProcess.options.label2.Type = 'label'; % Exclude contacts within a certain distance from the stimulation sites sProcess.options.excluderadius.Comment = 'Exclusion zone radius:
'; sProcess.options.excluderadius.Type = 'value'; sProcess.options.excluderadius.Value = {20,'mm', 0}; sProcess.options.label3.Comment = ['' ... - 'Exclude analysis of contacts within ' ... - 'this distance from the stimulation site']; + 'Exclude contacts within this distance ' ... + 'from the stimulation site']; sProcess.options.label3.Type = 'label'; sProcess.options.separator1.Type = 'separator'; % Plot window From 6a5f9599cba2f0b4da2779678b18d04839ae7536 Mon Sep 17 00:00:00 2001 From: rcassani Date: Mon, 27 Jul 2026 21:23:48 -0400 Subject: [PATCH 35/55] Bugfix: Left groups END with an apostrophe --- toolbox/process/functions/process_fastgraph.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 3aa91f6ed2..6dbb6e2ef2 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -425,8 +425,8 @@ % Use channel group names to assign hemisphere sContactGroupLocIdxs.Left = zeros(1, length(iSeeg)); for i = 1:length(iSeeg) - % Left groups start with an apostrophe - sContactGroupLocIdxs.Left(i) = strcmp(ChannelMat.Channel(iSeeg(i)).Group(1), ''''); + % Left groups end with an apostrophe + sContactGroupLocIdxs.Left(i) = strcmp(ChannelMat.Channel(iSeeg(i)).Group(end), ''''); end % Remaining contacts belong to the right hemisphere sContactGroupLocIdxs.Right = ~sContactGroupLocIdxs.Left; From 471cb406e69a4543cff20937cc38f293fc6f638e Mon Sep 17 00:00:00 2001 From: rcassani Date: Mon, 27 Jul 2026 21:24:19 -0400 Subject: [PATCH 36/55] Same output for `GroupSeegContacts` with or without locations --- toolbox/process/functions/process_fastgraph.m | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 6dbb6e2ef2..79cf65e0bc 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -416,7 +416,7 @@ % Split SEEG contacts into left and right hemisphere groups function sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat) % Initialize output structure - sContactGroupLocIdxs = struct(); + sContactGroupLocIdxs = struct('Left', [], 'Right', [], 'All', []); % Get index of SEEG channel type iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); % Check if valid location data is available @@ -428,14 +428,14 @@ % Left groups end with an apostrophe sContactGroupLocIdxs.Left(i) = strcmp(ChannelMat.Channel(iSeeg(i)).Group(end), ''''); end + sContactGroupLocIdxs.Left = find(sContactGroupLocIdxs.Left); % Remaining contacts belong to the right hemisphere - sContactGroupLocIdxs.Right = ~sContactGroupLocIdxs.Left; + sContactGroupLocIdxs.Right = find(~sContactGroupLocIdxs.Left); + % Combined contacts + sContactGroupLocIdxs.All = [sContactGroupLocIdxs.Left, sContactGroupLocIdxs.Right]; else % Store SEEG contact coordinates - contactLocs = zeros(length(iSeeg), 3); - for i = 1:length(iSeeg) - contactLocs(i, :) = ChannelMat.Channel(iSeeg(i)).Loc'; - end + contactLocs = cat(2, [ChannelMat.Channel(iSeeg).Loc])'; % Use coordinates to split contacts by hemisphere sContactGroupLocIdxs = SortLAPRAP(contactLocs); end From 077d56ffbb18e22c4f673438bcd7ceb3f66051da Mon Sep 17 00:00:00 2001 From: rcassani Date: Mon, 27 Jul 2026 21:26:45 -0400 Subject: [PATCH 37/55] Clean up --- toolbox/process/functions/process_fastgraph.m | 157 +++++++----------- 1 file changed, 61 insertions(+), 96 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 79cf65e0bc..0cafb4a06a 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -206,88 +206,61 @@ figure; % Maximize figure set(gcf, 'Position', get(0,'Screensize')); - % Shared y-axis limits across FastGraph subplots - commonAxisLimits = []; % Reserve one extra subplot for the legend - nSubplots = length(sInputs)+1; - % Define the plot parameters + nFastGraphs = length(sInputs); + hFastGraphAxes = gobjects(nFastGraphs, 0); % Subplot grid dimensions - nCols = ceil(sqrt(nSubplots)); - nRows = ceil(nSubplots / nCols); + nCols = ceil(sqrt(nFastGraphs+1)); + nRows = ceil((nFastGraphs+1) / nCols); % Subplot spacing and margins gap = [0.075 0.0175]; horzMargin = 0.03; vertMargin = 0.015; % Generate one FastGraph per selected input bst_progress('start', 'Process', 'Plotting FastGraphs...', 0, 100); - for iSubplot = 1:nSubplots-1 + for iFastGraph = 1:nFastGraphs % Show progress - progressPrc = round(100 .* iSubplot ./ (nSubplots-1)); + progressPrc = round(100 .* iFastGraph ./ nFastGraphs); bst_progress('set', progressPrc); % Data to be plotted for the current subplot subplotData = struct(); - Fout = seegData{iSubplot}.F(iSeeg, :); - % Keep only left-hemisphere channels if present - if any(sContactGroupLocIdxs.Left) - subplotData.leftData = Fout(sContactGroupLocIdxs.Left,:); - end - % Keep only left-hemisphere channels if present - if any(sContactGroupLocIdxs.Right) - subplotData.rightData = Fout(sContactGroupLocIdxs.Right,:); - end + Fout = seegData{iFastGraph}.F(iSeeg, :); + % Separate L and R hemisphere data + subplotData.leftData = Fout(sContactGroupLocIdxs.Left,:); + subplotData.rightData = Fout(sContactGroupLocIdxs.Right,:); % Sort channels within each hemisphere using the selected metric and time window sSubplotDataSorted = ApplyDataSorting(subplotData, seegData, OPTIONS); % Create the subplot with custom spacing - subtightplot(nRows, nCols, iSubplot, gap, horzMargin, vertMargin); + hFastGraphAxes(iFastGraph) = subtightplot(nRows, nCols, iFastGraph, gap, horzMargin, vertMargin); % Plot the FastGraph for the current stimulation pair - [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInputs, stimLocs, iSubplot, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); - % Tighten axes to the plotted data and store the current axis handle - axis tight - axisLimits = axis; - axSubplots(iSubplot) = gca; - % Update the shared y-axis limits so all FastGraph subplots can - % use the same vertical range for visual comparison - if iSubplot == 1 - commonAxisLimits = axisLimits; - else - commonAxisLimits(3) = min(commonAxisLimits(3), axisLimits(3)); - commonAxisLimits(4) = max(commonAxisLimits(4), axisLimits(4)); - end + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInputs, stimLocs, iFastGraph, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); % Apply edge transparency to the subplot - if exist('hLeftAreaPLot','var') - set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); - end - if exist('hRightAreaPLot','var') - set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); - end + set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); + set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); % Add the stimulation pair and atlas scout label as the subplot title - AddFastgraphTitle(sInputs, sSortedFastgraphLocIdxs.All(iSubplot), chanNamesSeeg, atlasScoutLabelsSeeg); - end - % Format all FastGraph axes - axFastGraphs = axSubplots(1:nSubplots-1); - % Set axis labels for all FastGraph subplots - set([axFastGraphs.XLabel], 'String', 'Time (ms)'); - set([axFastGraphs.YLabel], 'String', 'Voltage (mV)'); - % Apply common axes properties - set(axFastGraphs, 'YLim', commonAxisLimits(3:4), 'XAxisLocation', 'bottom'); - % Add zero-reference lines and hemisphere labels - DecorateFastgraphAxes(axFastGraphs); - % Link subplot axes so that zooming stays synchronized - linkaxes(axFastGraphs) - set(gcf,'units','normalized','outerposition',[0 0 1 1]) - zoom on - - % === Use the final subplot to display legend === + AddFastgraphTitle(sInputs, sSortedFastgraphLocIdxs.All(iFastGraph), chanNamesSeeg, atlasScoutLabelsSeeg); + end + + % === Common feature on FastGraph plots === + % Axes style + axis(hFastGraphAxes, 'tight'); + % Share axes, sets the same XY Limits + linkaxes(hFastGraphAxes, 'xy'); + % Set axis labels + xlabel(hFastGraphAxes, 'Time (ms)'); + ylabel(hFastGraphAxes, 'Voltage (mV)'); + % Line and label to distinguish hemispheres + SetHemisphereLabels(hFastGraphAxes); + + % === Plot brain legend === bst_progress('text', 'Plotting legend...'); % Generate a cortex snapshot with atlas scout for display imgCortex = GenerateCortexSnapshot(sInputs, OPTIONS); % Create the legend subplot with the same spacing settings - subtightplot(nRows, nCols, iSubplot+1, gap, horzMargin, vertMargin); + axBrain = subtightplot(nRows, nCols, nRows*nCols, gap, horzMargin, vertMargin); % Plot the reference panel with the cortex snapshot and axis labels - axSubplots(iSubplot+1) = gca; - PlotLegend(axSubplots(iSubplot+1), imgCortex, round(axSubplots(1).XLim), axSubplots(1).YLim); - - % Close progress + PlotLegend(axBrain, imgCortex, round(hFastGraphAxes(1).XLim), hFastGraphAxes(1).YLim); + % Close progress bst_progress('stop'); end @@ -341,7 +314,7 @@ % Contacts exactly on the midline (y == 0) are assigned to the left hemisphere. function sSortedLocIdxs = SortLAPRAP(contactLocs) % Initialize output structure - sSortedLocIdxs = struct(); + sSortedLocIdxs = struct('Left', [], 'Right', [], 'All', []); % Original row index of each location contactIdxs = (1:size(contactLocs, 1))'; % Append original row indices so duplicate coordinates keep input order @@ -354,14 +327,17 @@ rightContactLocs = contactLocsWithIdx(isRightHemisphere, :); % Sort left and right hemisphere contacts by x-coordinate in descending order (-xCoordColumn). % Use original index (idxColumn) as a secondary key so repeated locations remain grouped - % and keep their original input order + % and keep their original input order. Store sorted original indices for each hemisphere xCoordColumn = 1; - idxColumn = 4; - leftContactLocs = sortrows(leftContactLocs, [xCoordColumn, idxColumn], {'descend' 'ascend'}); - rightContactLocs = sortrows(rightContactLocs, [xCoordColumn, idxColumn], {'descend' 'ascend'}); - % Store sorted original indices for each hemisphere - sSortedLocIdxs.Left = leftContactLocs(:, 4)'; - sSortedLocIdxs.Right = rightContactLocs(:, 4)'; + idxColumn = 4; + if ~isempty(leftContactLocs) + leftContactLocs = sortrows(leftContactLocs, [xCoordColumn, idxColumn], {'descend' 'ascend'}); + sSortedLocIdxs.Left = leftContactLocs(:, 4)'; + end + if ~isempty(rightContactLocs) + rightContactLocs = sortrows(rightContactLocs, [xCoordColumn, idxColumn], {'descend' 'ascend'}); + sSortedLocIdxs.Right = rightContactLocs(:, 4)'; + end % Combined sorted indices sSortedLocIdxs.All = [sSortedLocIdxs.Left, sSortedLocIdxs.Right]; end @@ -735,38 +711,27 @@ function PlotLegend(axLegend, brainImg, xLim, yLim) %% ===== DECORATE FASTGRAPH AXES ===== % Add the zero-reference line and L/R hemisphere labels -function DecorateFastgraphAxes(axFastgraphs) - for ax = axFastgraphs - line(ax, ax.XLim, [0 0], ... - 'Color', [0 0 0], ... - 'LineWidth', 0.5, ... - 'HandleVisibility', 'off'); - AddHemisphereLabels(ax); +function SetHemisphereLabels(hFastGraphAxes) + % Positions for elements, assumes all hFastGraphAxes have the same XY Limits + XLim = hFastGraphAxes(1).XLim; + YLim = hFastGraphAxes(1).YLim; + xPositionText = XLim(1) + (0.95 * diff(XLim)); + yPositionText = min(abs(YLim)) * 0.15; + % Add 0 mV line and labels for L and R hemispheres + for hAxes = hFastGraphAxes + line(hAxes, XLim, [0 0], ... + 'Color', [0 0 0], ... + 'LineWidth', 0.5); + textProperties = { ... + 'Parent', hAxes, ... + 'FontSize', 14, ... + 'FontWeight', 'bold'}; + % Place text L/R above and below 0 mV line + text(xPositionText, yPositionText, 'L', textProperties{:}); + text(xPositionText, -yPositionText, 'R', textProperties{:}); end end -%% ===== ADD HEMISPHERE LABELS ===== -% Add L/R labels immediately above and below the y = 0 reference line -function AddHemisphereLabels(ax) - % Common label properties - labelProperties = { ... - 'Parent', ax, ... - 'Units', 'normalized', ... - 'FontSize', 8, ... - 'FontWeight', 'bold', ... - 'HorizontalAlignment', 'right', ... - 'VerticalAlignment', 'middle', ... - 'Clipping', 'off'}; - % Find the normalized vertical position of y = 0 - yZeroNormalized = -ax.YLim(1) / diff(ax.YLim); - % Position labels close to the right edge - xPosition = 0.96; - % Normalized vertical spacing around the zero line - labelGap = 0.035; - % Place labels above and below the bottom x-axis - text(xPosition, yZeroNormalized + labelGap, 'L', labelProperties{:}); - text(xPosition, yZeroNormalized - labelGap, 'R', labelProperties{:}); -end %% ===== UPDATE LEGEND IMAGE ===== % Update the overlay image position so it stays centered inside the From 05d472356b19810f58be41dc20fae402b7038c73 Mon Sep 17 00:00:00 2001 From: rcassani Date: Mon, 27 Jul 2026 21:27:05 -0400 Subject: [PATCH 38/55] Show figure when complete --- toolbox/process/functions/process_fastgraph.m | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 0cafb4a06a..2c0840a047 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -203,7 +203,8 @@ atlasScoutLabelsSeeg = chanTableWithAtlas(2:end, cols); % Create figure for FastGraph - figure; + hFig = figure; + hFig.Visible = 'off'; % Maximize figure set(gcf, 'Position', get(0,'Screensize')); % Reserve one extra subplot for the legend @@ -262,6 +263,9 @@ PlotLegend(axBrain, imgCortex, round(hFastGraphAxes(1).XLim), hFastGraphAxes(1).YLim); % Close progress bst_progress('stop'); + + % Show figure + hFig.Visible = 'on'; end %% ===== GET STIMULATION SITE CONTACT LOCATION ===== From 342c1d9e1a17fc98c36cda8ec069b6bbcde1506d Mon Sep 17 00:00:00 2001 From: rcassani Date: Tue, 28 Jul 2026 09:03:10 -0400 Subject: [PATCH 39/55] `GetOptions`, same order as in GUI process options --- toolbox/process/functions/process_fastgraph.m | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 2c0840a047..abd83cb6e9 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -105,31 +105,30 @@ function OPTIONS = GetOptions(sProcess) OPTIONS = struct(); % Atlas and scouts to use for plotting FastGraph - ScoutsList = sProcess.options.scouts.Value; - OPTIONS.Atlas = ScoutsList{1,1}; - OPTIONS.AtlasScoutLabels = ScoutsList{1,2}; + OPTIONS.Atlas = sProcess.options.scouts.Value{1,1}; + OPTIONS.AtlasScoutLabels = sProcess.options.scouts.Value{1,2}; % Color figure by region or by label OPTIONS.ColorScheme = sProcess.options.colorscheme.Value; % Select regions to include OPTIONS.Region = sProcess.options.region.Value; % Method for sorting the data OPTIONS.SortMethod = sProcess.options.sortmethod.Value; - % Sort window + % Time window for sorting the data [s] if isfield(sProcess.options, 'sortwindow') && isfield(sProcess.options.sortwindow, 'Value') && iscell(sProcess.options.sortwindow.Value) && ~isempty(sProcess.options.sortwindow.Value) OPTIONS.SortWindow = round((sProcess.options.sortwindow.Value{1} * 1000)) + 101; else OPTIONS.SortWindow = []; end - % Plot window + % Exclude contacts within a certain distance of stimulation sites + OPTIONS.ExcludeRadius = sProcess.options.excluderadius.Value{1}; + % Time window for plotting [s] if isfield(sProcess.options, 'plotwindow') && isfield(sProcess.options.plotwindow, 'Value') && iscell(sProcess.options.plotwindow.Value) && ~isempty(sProcess.options.plotwindow.Value) OPTIONS.PlotWindow = round((sProcess.options.plotwindow.Value{1} * 1000)); else OPTIONS.PlotWindow = []; end - % Edge transparency of plot + % Edge transparency for plotting OPTIONS.EdgeAlpha = sProcess.options.edgealpha.Value{1}; - % Exclude contacts within a certain distance of stimulation sites - OPTIONS.ExcludeRadius = sProcess.options.excluderadius.Value{1}; end %% ===== RUN ===== From 8ba8ad445c10f1c54fe7ef3c3bab3116e012145f Mon Sep 17 00:00:00 2001 From: rcassani Date: Tue, 28 Jul 2026 09:17:40 -0400 Subject: [PATCH 40/55] Bugfix: Do not use hardcoded values to compute --- toolbox/process/functions/process_fastgraph.m | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index abd83cb6e9..db5c57b0f7 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -115,7 +115,7 @@ OPTIONS.SortMethod = sProcess.options.sortmethod.Value; % Time window for sorting the data [s] if isfield(sProcess.options, 'sortwindow') && isfield(sProcess.options.sortwindow, 'Value') && iscell(sProcess.options.sortwindow.Value) && ~isempty(sProcess.options.sortwindow.Value) - OPTIONS.SortWindow = round((sProcess.options.sortwindow.Value{1} * 1000)) + 101; + OPTIONS.SortWindow = sProcess.options.sortwindow.Value{1}; else OPTIONS.SortWindow = []; end @@ -123,7 +123,7 @@ OPTIONS.ExcludeRadius = sProcess.options.excluderadius.Value{1}; % Time window for plotting [s] if isfield(sProcess.options, 'plotwindow') && isfield(sProcess.options.plotwindow, 'Value') && iscell(sProcess.options.plotwindow.Value) && ~isempty(sProcess.options.plotwindow.Value) - OPTIONS.PlotWindow = round((sProcess.options.plotwindow.Value{1} * 1000)); + OPTIONS.PlotWindow = sProcess.options.plotwindow.Value{1}; else OPTIONS.PlotWindow = []; end @@ -430,7 +430,7 @@ if isempty(OPTIONS.SortWindow) sortWindowIdx = 1:size(seegData{1}.F,2); else - sortWindowIdx = OPTIONS.SortWindow(1):OPTIONS.SortWindow(2); + sortWindowIdx = bst_closest(OPTIONS.SortWindow, seegData{1}.Time); end % Sort channels within each hemisphere using the selected metric switch lower(OPTIONS.SortMethod) @@ -481,7 +481,11 @@ % Check whether stimulation locations are available hasStimLocs = any(stimLocs(:)); % Select the time samples to display - plotWindowIdx = OPTIONS.PlotWindow(1) + 101 : OPTIONS.PlotWindow(2) + 101; + if isempty(OPTIONS.PlotWindow) + plotWindowIdx = 1:size(seegData{1}.F,2); + else + plotWindowIdx = bst_closest(OPTIONS.PlotWindow, seegData{1}.Time); + end timeMs = seegData{iSubplot}.Time(plotWindowIdx) * 1000; fprintf('\n===== FastGraph %d/%d: Stimulation site "%s" =====\n', iSubplot, numel(sInputs), sInputs(iSubplot).Comment) From 4f964e635bc5f62f97ea473199e9c7d0c44116da Mon Sep 17 00:00:00 2001 From: rcassani Date: Tue, 28 Jul 2026 09:18:33 -0400 Subject: [PATCH 41/55] Bugfix: Check for regions only for `region` color scheme --- toolbox/process/functions/process_fastgraph.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index db5c57b0f7..af40e25587 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -174,8 +174,8 @@ % Get options OPTIONS = GetOptions(sProcess); - % Early exit if no region is selected - if isempty(OPTIONS.Region) + % Check regions for 'region' color scheme + if strcmpi(OPTIONS.ColorScheme, 'region') && isempty(OPTIONS.Region) bst_report('Error', sProcess, [], 'No region selected. Select at least one region to run the analysis.'); return; end From 2f16fca3ecb1f5ddcf8f94b8a5a96e3a2a51250a Mon Sep 17 00:00:00 2001 From: rcassani Date: Tue, 28 Jul 2026 09:21:31 -0400 Subject: [PATCH 42/55] GUI: Input time windows in `[ms]` --- toolbox/process/functions/process_fastgraph.m | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index af40e25587..e90155c7a2 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -70,7 +70,7 @@ % Sort window sProcess.options.sortwindow.Comment = 'Time window to sort data: '; sProcess.options.sortwindow.Type = 'timewindow'; -sProcess.options.sortwindow.Value = []; +sProcess.options.sortwindow.Value = {[], 'ms', []}; sProcess.options.label2.Comment = ['' ... 'Examples:
'... 'Early latency: 0-60 ms,
' ... @@ -89,18 +89,20 @@ % Plot window sProcess.options.plotwindow.Comment = 'Plot time range: '; sProcess.options.plotwindow.Type = 'timewindow'; -sProcess.options.plotwindow.Value = []; +sProcess.options.plotwindow.Value = {[], 'ms', []}; % Edge transparency of plot sProcess.options.edgealpha.Comment = 'Edge transparency of plot: '; sProcess.options.edgealpha.Type = 'value'; sProcess.options.edgealpha.Value = {0.05,' ', 2}; end + %% ===== FORMAT COMMENT ===== function Comment = FormatComment(sProcess) %#ok Comment = sProcess.Comment; end + %% ===== GET OPTIONS ===== function OPTIONS = GetOptions(sProcess) OPTIONS = struct(); @@ -131,12 +133,13 @@ OPTIONS.EdgeAlpha = sProcess.options.edgealpha.Value{1}; end + %% ===== RUN ===== function OutputFiles = Run(sProcess, sInputs) %#ok % Initialize output OutputFiles = {}; - % Check that all input files use the same channel file + % ===== Check that all input files use the same channel file ===== ChannelFiles = {sInputs.ChannelFile}; if length(unique(ChannelFiles)) > 1 bst_report('Error', sProcess, sInputs, 'All input files must use the same channel file.'); @@ -156,7 +159,7 @@ return; end - % ===== Check whether all channel names in comment are valid ===== + % ===== Check channels for bipolar channels are valid channel names ===== % Load the channel file ChannelMat = in_bst_channel(ChannelFiles{1}); channelNames = {ChannelMat.Channel.Name}; From 0189ccbbf7cfb8ff78cb7766d765c08786c35a78 Mon Sep 17 00:00:00 2001 From: rcassani Date: Tue, 28 Jul 2026 12:37:26 -0400 Subject: [PATCH 43/55] Export channel atlas: Allow filtering which atlases to use --- toolbox/io/export_channel_atlas.m | 42 ++++++++++++++++--- toolbox/process/functions/process_fastgraph.m | 2 +- 2 files changed, 37 insertions(+), 7 deletions(-) diff --git a/toolbox/io/export_channel_atlas.m b/toolbox/io/export_channel_atlas.m index d2a759f89f..bbdc50ae68 100644 --- a/toolbox/io/export_channel_atlas.m +++ b/toolbox/io/export_channel_atlas.m @@ -1,8 +1,8 @@ -function [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, Modality, TsvFile, Radius, isProba, isInteractive) +function [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, Modality, TsvFile, Radius, isProba, isInteractive, AtlasFilter) % EXPORT_CHANNEL_ATLAS: Compute anatomical labels for SEEG/ECOG contacts from volume and surface parcellations % -% USAGE: [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, Modality='ECOG+SEEG', TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1) -% [TsvFile, ChanTable] = export_channel_atlas(ChannelFile, iChannels, TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1) +% USAGE: TsvFile = export_channel_atlas(ChannelFile, Modality='ECOG+SEEG', TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1, AtlasFilter='') +% TsvFile = export_channel_atlas(ChannelFile, iChannels, TsvFile=[ask], Radius=[ask], isProba=[ask], isInteractive=1, AtlasFilter='') % % INPUT: % - ChannelFile : Path to Brainstorm channel file to be processed @@ -10,10 +10,12 @@ % - iChannels : Array of integers, export only the selected channel indices % - TsvFile : Output text file (tab-separated values) % - Radius : Size in millimeters of the neighborhood to consider around each contact +% - isProba : If 1, for each volume atlas, add a column indicating the spatial probability (100 * nVoxelWithLabel / nVoxelsInSphere) % - IsInteractive : If 1, display the output table at the end of the process -% : If 0, use all available Coodinates, Parcellations (anat) and Atlases (surface), -% and do not display output table -% - iChannels : Limit export to a subset of channel indices +% : If 0, use all available Coodinates, and the Parcellations (anat) and Atlases (surface) filtered by 'AtlasFilter' +% and do not display output table +% - AtlasFilter : Used with IsInteractive=0. If absent or empty, use all the Parcellations (anat) and Atlases (surface). +% Otherwise use the string in AtlasFilter to indicate the Parcellations and Atlases to return. % OUTPUT: % - TsvFile : Output text file (tab-separated values). Empty when no file was selected or requested. % - ChanTable : Cell array containing the complete output table. The same information is written to TsvFile @@ -54,6 +56,9 @@ % ===== PASRSE INPUTS ===== +if (nargin < 7) || isempty(AtlasFilter) + AtlasFilter = []; +end if (nargin < 6) || isempty(isInteractive) isInteractive = 1; end @@ -247,6 +252,31 @@ end +% ===== FILTER ATLASES ===== +if ~isInteractive && ~isempty(AtlasFilter) + % === Volume === + iColVol = find(~cellfun(@(c)isempty(strfind(c, tagVol)), Columns(:,2))); + % Try match, if nothing found, try regular expression + iVolValid = find(strcmpi(Columns(iColVol,1), AtlasFilter)); + if isempty(iVolValid) + iVolValid = find(~cellfun(@isempty, regexp(Columns(iColVol,1), AtlasFilter))); + end + % Remove non-matching volume atlases + Columns(setdiff(iColVol, iColVol(iVolValid)), :) = []; + % Update iColSurf in case volumes were removed + iColSurf = find(~cellfun(@(c)isempty(strfind(c, tagSurf)), Columns(:,2))); + + % === Surface === + % Try match, if nothing found, try regular expression + iSurfValid = find(strcmpi(SurfAtlases, AtlasFilter)); + if isempty(iSurfValid) + iSurfValid = find(~cellfun(@isempty, regexp(SurfAtlases, AtlasFilter))); + end + % Keep only matching surface atlases + SurfAtlases = SurfAtlases(iSurfValid); +end + + % ===== SPHERE PROBE ===== % Ask sphere radius to users if isempty(Radius) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index e90155c7a2..f8b5f724e6 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -195,7 +195,7 @@ % Split SEEG contacts into left and right hemisphere groups sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat); % Compute anatomical labels for the contacts from volume/surface parcellations - [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 5, 0, 0); + [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 5, 0, 0, OPTIONS.Atlas); % Locate atlas related columns from channel table above hit = cellfun(@(x) ischar(x) && (~isempty(strfind(OPTIONS.Atlas, x)) || ~isempty(strfind(x, OPTIONS.Atlas))), chanTableWithAtlas(1,:)); % Columns whose header matches the atlas name From 7813043fd05333d84f62a1c5496003a4bc72662b Mon Sep 17 00:00:00 2001 From: rcassani Date: Tue, 28 Jul 2026 12:51:42 -0400 Subject: [PATCH 44/55] Bugfix: Do not use hardcoded values to compute window samples (Part2) --- toolbox/process/functions/process_fastgraph.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index f8b5f724e6..ea74b3fd8b 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -434,6 +434,7 @@ sortWindowIdx = 1:size(seegData{1}.F,2); else sortWindowIdx = bst_closest(OPTIONS.SortWindow, seegData{1}.Time); + sortWindowIdx = [sortWindowIdx(1):sortWindowIdx(2)]; end % Sort channels within each hemisphere using the selected metric switch lower(OPTIONS.SortMethod) @@ -488,6 +489,7 @@ plotWindowIdx = 1:size(seegData{1}.F,2); else plotWindowIdx = bst_closest(OPTIONS.PlotWindow, seegData{1}.Time); + plotWindowIdx = [plotWindowIdx(1):plotWindowIdx(2)]; end timeMs = seegData{iSubplot}.Time(plotWindowIdx) * 1000; From 87aa6ae8031016a542f55f9325421ddf4a7d7dbc Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 10:38:48 -0400 Subject: [PATCH 45/55] Use Group name to sort SEEG contacts IFF for SEEG contacts w/o Loc --- toolbox/process/functions/process_fastgraph.m | 38 ++++++++----------- 1 file changed, 16 insertions(+), 22 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index ea74b3fd8b..f879c382b2 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -396,31 +396,25 @@ %% ===== SPLIT CONTACTS TO LEFT/RIGHT HEMISPHERE ===== % Split SEEG contacts into left and right hemisphere groups -function sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat) - % Initialize output structure - sContactGroupLocIdxs = struct('Left', [], 'Right', [], 'All', []); +function sContactGroupLocIdxs = GroupSeegContacts(ChannelMat) % Get index of SEEG channel type - iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); - % Check if valid location data is available - noLocations = isempty(stimLocs) || ~any(stimLocs(:)); - if noLocations - % Use channel group names to assign hemisphere - sContactGroupLocIdxs.Left = zeros(1, length(iSeeg)); - for i = 1:length(iSeeg) - % Left groups end with an apostrophe - sContactGroupLocIdxs.Left(i) = strcmp(ChannelMat.Channel(iSeeg(i)).Group(end), ''''); + iSeegs = channel_find(ChannelMat.Channel, 'SEEG'); + % For each SEEG Concact, if no valid location, add temporary Loc based on the Group Name + for ix = 1 : iSeegs + iSeeg = iSeegs(ix); + if isempty(ChannelMat.Channel(iSeeg).Loc) || all(ChannelMat.Channel(iSeeg).Loc == 0) || any(isnan(ChannelMat.Channel(iSeeg).Loc)) + % SEEG groups in Left hemisphere end with an apostrophe + if strcmp(ChannelMat.Channel(iSeeg).Group(end), '''') + ChannelMat.Channel(iSeeg).Loc = [ 1; 0; 0]; % Left hemisphere in SCS + else + ChannelMat.Channel(iSeeg).Loc = [-1; 0; 0]; % Right hemisphere in SCS + end end - sContactGroupLocIdxs.Left = find(sContactGroupLocIdxs.Left); - % Remaining contacts belong to the right hemisphere - sContactGroupLocIdxs.Right = find(~sContactGroupLocIdxs.Left); - % Combined contacts - sContactGroupLocIdxs.All = [sContactGroupLocIdxs.Left, sContactGroupLocIdxs.Right]; - else - % Store SEEG contact coordinates - contactLocs = cat(2, [ChannelMat.Channel(iSeeg).Loc])'; - % Use coordinates to split contacts by hemisphere - sContactGroupLocIdxs = SortLAPRAP(contactLocs); end + % Store SEEG contact coordinates + contactLocs = cat(2, [ChannelMat.Channel(iSeegs).Loc])'; + % Use coordinates to split contacts by hemisphere + sContactGroupLocIdxs = SortLAPRAP(contactLocs); end %% ===== WITHIN-HEMISPHERE DATA SORTING ===== From 3c1ba570109fa138fbfd1142d8c798ec3fab24dd Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 10:39:05 -0400 Subject: [PATCH 46/55] Clean up --- toolbox/process/functions/process_fastgraph.m | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index f879c382b2..d5c7ac115d 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -271,9 +271,9 @@ end %% ===== GET STIMULATION SITE CONTACT LOCATION ===== -% Get the midpoint location of each stimulation pair from channel +% Get the midpoint location of each stimulation pair in Comment function stimLocs = GetStimLocs(sInputs, ChannelMat) - % Preallocate one [x y z] midpoint per stimulation pair + % Preallocate one [x y z] SCS midpoint per stimulation pair stimLocs = zeros(numel(sInputs), 3); % Get channel names once for lookup chanNames = {ChannelMat.Channel.Name}; @@ -285,12 +285,9 @@ if numel(parts) ~= 2 continue; end - % Clean extracted comment - contact1 = parts{1}; - contact2 = parts{2}; % Get the contact names - contact1Parts = strsplit(contact1); - contact2Parts = strsplit(contact2); + contact1Parts = strsplit(parts{1}, ' '); + contact2Parts = strsplit(parts{2}, ' '); contact1 = contact1Parts{end}; contact2 = contact2Parts{1}; % Find the channel indices From aba0a53cb64b8dc94bd671f094c1ab603f825757 Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 10:39:17 -0400 Subject: [PATCH 47/55] Reorganize --- toolbox/process/functions/process_fastgraph.m | 62 +++++++++++-------- 1 file changed, 35 insertions(+), 27 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index d5c7ac115d..1db9a201a2 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -138,7 +138,15 @@ function OutputFiles = Run(sProcess, sInputs) %#ok % Initialize output OutputFiles = {}; + % Get options + OPTIONS = GetOptions(sProcess); + % ===== Check regions for 'region' color scheme ===== + if strcmpi(OPTIONS.ColorScheme, 'region') && isempty(OPTIONS.Region) + bst_report('Error', sProcess, [], 'No region selected. Select at least one region to run the analysis.'); + return; + end + % ===== Check that all input files use the same channel file ===== ChannelFiles = {sInputs.ChannelFile}; if length(unique(ChannelFiles)) > 1 @@ -174,27 +182,17 @@ return; end - % Get options - OPTIONS = GetOptions(sProcess); - - % Check regions for 'region' color scheme - if strcmpi(OPTIONS.ColorScheme, 'region') && isempty(OPTIONS.Region) - bst_report('Error', sProcess, [], 'No region selected. Select at least one region to run the analysis.'); - return; - end - - % Get indices of SEEG channels - iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); - % Get the midpoint location of each stimulation pair from channel - stimLocs = GetStimLocs(sInputs, ChannelMat); - % Sort FastGraphs by stimulation-site location for LAPRAP style display - sSortedFastgraphLocIdxs = SortLAPRAP(stimLocs); - - % Load SEEG recordings after applying FastGraph sorting - [seegData, excludedContacts] = GetSeegData(sInputs, sSortedFastgraphLocIdxs, stimLocs, ChannelMat, OPTIONS); - % Split SEEG contacts into left and right hemisphere groups - sContactGroupLocIdxs = GroupSeegContacts(stimLocs, ChannelMat); - % Compute anatomical labels for the contacts from volume/surface parcellations + % === Sort SEEG contacts and Stimulus by location + % Sort ONLY SEEG Contacts using Location or Group name into Left | Right groups + sContactLocIdxs = GroupSeegContacts(ChannelMat); % iSeeg(sCon) + % Get stimulus location from each Input, based on the Comment + stimScsLocs = GetStimLocs(sInputs, ChannelMat); + % Sort stimulus using Location into Left | Right groups, Anterior->Posterior within group + sStimLocIdxs = SortLAPRAP(stimScsLocs); + % Sort Inputs by their Stimulus location + % sInputs = sInputs(sStimLocIdxs); + + % === Anatomical labels for SEEG contacts [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 5, 0, 0, OPTIONS.Atlas); % Locate atlas related columns from channel table above hit = cellfun(@(x) ischar(x) && (~isempty(strfind(OPTIONS.Atlas, x)) || ~isempty(strfind(x, OPTIONS.Atlas))), chanTableWithAtlas(1,:)); @@ -204,12 +202,18 @@ chanNamesSeeg = chanTableWithAtlas(2:end, 1); atlasScoutLabelsSeeg = chanTableWithAtlas(2:end, cols); - % Create figure for FastGraph + % Get indices of SEEG channels + iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); + % Load ONLY SEEG recordings after applying FastGraph sorting + [seegData, excludedContacts] = GetSeegData(sInputs, sStimLocIdxs, stimScsLocs, ChannelMat, OPTIONS); + + + % ===== Create figure for FastGraph ===== hFig = figure; hFig.Visible = 'off'; % Maximize figure set(gcf, 'Position', get(0,'Screensize')); - % Reserve one extra subplot for the legend + % Reserve one extra subplot for the legend (brain surface with Scouts) nFastGraphs = length(sInputs); hFastGraphAxes = gobjects(nFastGraphs, 0); % Subplot grid dimensions @@ -221,27 +225,31 @@ vertMargin = 0.015; % Generate one FastGraph per selected input bst_progress('start', 'Process', 'Plotting FastGraphs...', 0, 100); + + % ===== Get data and Plot each FastGraph and ===== for iFastGraph = 1:nFastGraphs % Show progress progressPrc = round(100 .* iFastGraph ./ nFastGraphs); bst_progress('set', progressPrc); + + % Data to be plotted for the current subplot subplotData = struct(); Fout = seegData{iFastGraph}.F(iSeeg, :); % Separate L and R hemisphere data - subplotData.leftData = Fout(sContactGroupLocIdxs.Left,:); - subplotData.rightData = Fout(sContactGroupLocIdxs.Right,:); + subplotData.leftData = Fout(sContactLocIdxs.Left,:); + subplotData.rightData = Fout(sContactLocIdxs.Right,:); % Sort channels within each hemisphere using the selected metric and time window sSubplotDataSorted = ApplyDataSorting(subplotData, seegData, OPTIONS); % Create the subplot with custom spacing hFastGraphAxes(iFastGraph) = subtightplot(nRows, nCols, iFastGraph, gap, horzMargin, vertMargin); % Plot the FastGraph for the current stimulation pair - [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInputs, stimLocs, iFastGraph, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInputs, stimScsLocs, iFastGraph, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); % Apply edge transparency to the subplot set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); % Add the stimulation pair and atlas scout label as the subplot title - AddFastgraphTitle(sInputs, sSortedFastgraphLocIdxs.All(iFastGraph), chanNamesSeeg, atlasScoutLabelsSeeg); + AddFastgraphTitle(sInputs, sStimLocIdxs.All(iFastGraph), chanNamesSeeg, atlasScoutLabelsSeeg); end % === Common feature on FastGraph plots === From 8ce75c3119a77d9190f7f5103b1ce958044b5d15 Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 12:32:26 -0400 Subject: [PATCH 48/55] Refactor - Avoid loading data from all Inputs at the same time - Load only SEEG data, ignore other channels --- toolbox/process/functions/process_fastgraph.m | 169 ++++++++---------- 1 file changed, 75 insertions(+), 94 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 1db9a201a2..7952c0e97a 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -184,13 +184,13 @@ % === Sort SEEG contacts and Stimulus by location % Sort ONLY SEEG Contacts using Location or Group name into Left | Right groups - sContactLocIdxs = GroupSeegContacts(ChannelMat); % iSeeg(sCon) - % Get stimulus location from each Input, based on the Comment - stimScsLocs = GetStimLocs(sInputs, ChannelMat); - % Sort stimulus using Location into Left | Right groups, Anterior->Posterior within group - sStimLocIdxs = SortLAPRAP(stimScsLocs); - % Sort Inputs by their Stimulus location - % sInputs = sInputs(sStimLocIdxs); + sContactLocIdxs = SortSeegContacts(ChannelMat); + % Sort stimulus location from each Input into Left | Right groups, Anterior->Posterior within group + stimLocs = GetStimLocs(sInputs, ChannelMat); + sStimLocIdxs = SortLAPRAP(stimLocs); + % Sort Inputs and StimLocs by their Stimulus location: I.e. SubPlot order + sInputs = sInputs(sStimLocIdxs.All); + stimLocs = stimLocs(sStimLocIdxs.All, :); % === Anatomical labels for SEEG contacts [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 5, 0, 0, OPTIONS.Atlas); @@ -202,12 +202,6 @@ chanNamesSeeg = chanTableWithAtlas(2:end, 1); atlasScoutLabelsSeeg = chanTableWithAtlas(2:end, cols); - % Get indices of SEEG channels - iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); - % Load ONLY SEEG recordings after applying FastGraph sorting - [seegData, excludedContacts] = GetSeegData(sInputs, sStimLocIdxs, stimScsLocs, ChannelMat, OPTIONS); - - % ===== Create figure for FastGraph ===== hFig = figure; hFig.Visible = 'off'; @@ -228,28 +222,32 @@ % ===== Get data and Plot each FastGraph and ===== for iFastGraph = 1:nFastGraphs - % Show progress - progressPrc = round(100 .* iFastGraph ./ nFastGraphs); - bst_progress('set', progressPrc); + sInput = sInputs(iFastGraph); + stimLoc = stimLocs(iFastGraph, :); + % Show progress + bst_progress('set', round(100 .* (iFastGraph-1) ./ nFastGraphs)); + fprintf('\n===== FastGraph %d/%d: Stimulation site "%s" =====\n', iFastGraph, nFastGraphs, sInput.Comment); + % Load ONLY SEEG recordings + [seegData, excludedContacts] = GetSeegData(sInput, stimLoc, ChannelMat, OPTIONS); % Data to be plotted for the current subplot subplotData = struct(); - Fout = seegData{iFastGraph}.F(iSeeg, :); % Separate L and R hemisphere data - subplotData.leftData = Fout(sContactLocIdxs.Left,:); - subplotData.rightData = Fout(sContactLocIdxs.Right,:); + subplotData.leftData = seegData.F(sContactLocIdxs.Left,:); + subplotData.rightData = seegData.F(sContactLocIdxs.Right,:); % Sort channels within each hemisphere using the selected metric and time window sSubplotDataSorted = ApplyDataSorting(subplotData, seegData, OPTIONS); + % Create the subplot with custom spacing hFastGraphAxes(iFastGraph) = subtightplot(nRows, nCols, iFastGraph, gap, horzMargin, vertMargin); % Plot the FastGraph for the current stimulation pair - [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInputs, stimScsLocs, iFastGraph, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); % Apply edge transparency to the subplot set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); % Add the stimulation pair and atlas scout label as the subplot title - AddFastgraphTitle(sInputs, sStimLocIdxs.All(iFastGraph), chanNamesSeeg, atlasScoutLabelsSeeg); + AddFastgraphTitle(sInput, chanNamesSeeg, atlasScoutLabelsSeeg); end % === Common feature on FastGraph plots === @@ -264,6 +262,7 @@ SetHemisphereLabels(hFastGraphAxes); % === Plot brain legend === + bst_progress('set', 100); bst_progress('text', 'Plotting legend...'); % Generate a cortex snapshot with atlas scout for display imgCortex = GenerateCortexSnapshot(sInputs, OPTIONS); @@ -356,52 +355,42 @@ %% ===== LOAD AND FILTER SEEG DATA ===== % Load each selected SEEG block and optionally exclude contacts based on % distance from the stimulation site -function [seegData, excludedContacts] = GetSeegData(sInputs, sSortedFastgraphLocIdxs, stimLocs, ChannelMat, OPTIONS) - % Intialize output - seegData = cell(numel(sInputs), 1); - excludedContacts = cell(numel(sInputs), 1); +function [seegData, excludedContacts] = GetSeegData(sInput, stimLoc, ChannelMat, OPTIONS) % Get index of SEEG channel types iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); - for k = 1:numel(sInputs) - % Load current file - data = load(file_fullpath(sInputs(sSortedFastgraphLocIdxs.All(k)).FileName)); - % Mark bad channels as NaN - data.F(data.ChannelFlag<0, :) = NaN; - if ~isempty(stimLocs) - % Current stimulation center - stimCenter = stimLocs(sSortedFastgraphLocIdxs.All(k), :); - % Compute distance from stimulation site to each SEEG contact (mm) - contactDist = zeros(1, numel(ChannelMat.Channel)); - for j = iSeeg - contactDist(j) = norm(stimCenter - ChannelMat.Channel(j).Loc', 2) * 1000; - end - % Exclude stimulation contacts themselves - iStimContacts = (contactDist > 0) & (contactDist <= 2); - % Exclude contacts within user-provided distance from the stimulation sites - % iExcluded = contactDist > OPTIONS.ExcludeRadius; - iExcluded = (contactDist > 2) & (contactDist <= OPTIONS.ExcludeRadius); - % Keep only valid SEEG contacts - isSeeg = strcmp('SEEG',{ChannelMat.Channel.Type}); - validContacts = isSeeg & ~iExcluded & ~iStimContacts; - excludedContacts{k} = ~validContacts; - % Report excluded contacts - fprintf('Contacts excluded for being within the %d mm exclusion zone "%s":\n', OPTIONS.ExcludeRadius, sInputs(k).Comment); - fprintf('%s %s ', ChannelMat.Channel(iStimContacts).Name, ChannelMat.Channel(iExcluded).Name); - fprintf('\n\n'); - % Remove excluded channels - data.F(excludedContacts{k}, :) = NaN; - else - % If no stimulation locations are available, keep only SEEG channels - excludedContacts{k} = ~iSeeg; - end - % Store SEEG data for the current block - seegData{k} = data; + % Load current file + allData = load(file_fullpath(sInput.FileName)); + % Set data from bad channels to NaN + allData.F(allData.ChannelFlag<0, :) = NaN; + % Keep only data from SEEG + seegData.F = allData.F(iSeeg, :); + seegData.Time = allData.Time; + if all(stimLoc ~= 0) + % Compute distance [mm] from stimulation site to each SEEG contact (mm) + contactLocs = cat(2, [ChannelMat.Channel(iSeeg).Loc])'; + contactDist = sqrt(sum((contactLocs - repmat(stimLoc, length(iSeeg), 1)).^2, 2)) * 1000; + % Exclude stimulation contacts themselves + iStimContacts = (contactDist >= 0) & (contactDist <= 2); + % Exclude contacts within user-provided distance from the stimulation sites + iExcluded = (contactDist > 2) & (contactDist <= OPTIONS.ExcludeRadius); + % Keep only valid SEEG contacts + validContacts = ~iExcluded & ~iStimContacts; + excludedContacts = ~validContacts; + % Report excluded contacts + fprintf('Contacts excluded for being within the %d mm exclusion zone "%s":\n', OPTIONS.ExcludeRadius, sInput.Comment); + fprintf('%s %s ', ChannelMat.Channel(iSeeg(iStimContacts)).Name, ChannelMat.Channel(iSeeg(iExcluded)).Name); + fprintf('\n\n'); + % Set data from excluded channels to NaN + seegData.F(excludedContacts, :) = NaN; + else + % If no stimulation locations are available, keep only SEEG channels + excludedContacts = ~iSeeg; end end -%% ===== SPLIT CONTACTS TO LEFT/RIGHT HEMISPHERE ===== -% Split SEEG contacts into left and right hemisphere groups -function sContactGroupLocIdxs = GroupSeegContacts(ChannelMat) +%% ===== SORT CONTACTS INTO LEFT/RIGHT HEMISPHERE ===== +% Sort SEEG contacts into left and right hemisphere groups +function sContactGroupLocIdxs = SortSeegContacts(ChannelMat) % Get index of SEEG channel type iSeegs = channel_find(ChannelMat.Channel, 'SEEG'); % For each SEEG Concact, if no valid location, add temporary Loc based on the Group Name @@ -430,9 +419,9 @@ sSorted = struct(); % Get sample indices used for sorting if isempty(OPTIONS.SortWindow) - sortWindowIdx = 1:size(seegData{1}.F,2); + sortWindowIdx = 1:size(seegData.F,2); else - sortWindowIdx = bst_closest(OPTIONS.SortWindow, seegData{1}.Time); + sortWindowIdx = bst_closest(OPTIONS.SortWindow, seegData.Time); sortWindowIdx = [sortWindowIdx(1):sortWindowIdx(2)]; end % Sort channels within each hemisphere using the selected metric @@ -466,36 +455,29 @@ % Create one FastGraph subplot. % Left-hemisphere SEEG channels are plotted as positive stacked areas % Right-hemisphere SEEG channels are plotted as negative stacked areas -function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInputs, stimLocs, iSubplot, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) +function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) % Initialize output handles hLeftAreaPlot = []; hRightAreaPlot = []; % Get cortex to be used for region/color lookup - sSubject = bst_get('Subject', sInputs(1).SubjectName); + sSubject = bst_get('Subject', sInput.SubjectName); CortexFile = sSubject.Surface(sSubject.iCortex).FileName; sCortex = bst_memory('LoadSurface', CortexFile); % Resolve selected scouts selectedScoutLabels = ResolveScoutSelection(sCortex, OPTIONS); - % Get indices of all SEEG channels - iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); - % Match channel names against atlas table names - [~, iChanLocs] = ismember({ChannelMat.Channel.Name}, chanNamesSeeg); - % Check whether stimulation locations are available - hasStimLocs = any(stimLocs(:)); + % Check whether stimulation location is available + hasStimLocs = any(stimLoc); + % Select the time samples to display if isempty(OPTIONS.PlotWindow) - plotWindowIdx = 1:size(seegData{1}.F,2); + plotWindowIdx = 1:size(seegData.F,2); else - plotWindowIdx = bst_closest(OPTIONS.PlotWindow, seegData{1}.Time); + plotWindowIdx = bst_closest(OPTIONS.PlotWindow, seegData.Time); plotWindowIdx = [plotWindowIdx(1):plotWindowIdx(2)]; end - timeMs = seegData{iSubplot}.Time(plotWindowIdx) * 1000; + timeMs = seegData.Time(plotWindowIdx) * 1000; - fprintf('\n===== FastGraph %d/%d: Stimulation site "%s" =====\n', iSubplot, numel(sInputs), sInputs(iSubplot).Comment) - % Extract SEEG data once for this subplot - Fout = seegData{iSubplot}.F(iSeeg, :); - % Loop over left and right hemispheres for iSide = 1:2 if iSide == 1 @@ -520,29 +502,28 @@ % Reorder SEEG channels for the current hemisphere contactIdxs = groupLocIdxs(sortedIdxs); - plotLocs = iSeeg(contactIdxs); - hemiData = abs(Fout(contactIdxs, :)); + hemiData = abs(seegData.F(contactIdxs, :)); % Get atlas scout labels for these channels - channelScoutLabels = cell(1, numel(plotLocs)); - for i = 1:numel(plotLocs) - channelScoutLabels{i} = atlasScoutLabelsSeeg{iChanLocs(plotLocs(i))}; + channelScoutLabels = cell(1, numel(contactIdxs)); + for i = 1:numel(contactIdxs) + channelScoutLabels{i} = atlasScoutLabelsSeeg{contactIdxs(i)}; end % Filter channels using resolved scout selection if hasStimLocs toPlot = ismember(channelScoutLabels, selectedScoutLabels); else - toPlot = true(1, numel(plotLocs)); + toPlot = true(1, numel(contactIdxs)); end % Keep track of number of channel before filtering - nChannelsBeforeFilter = numel(plotLocs); + nChannelsBeforeFilter = numel(contactIdxs); % Keep only channels that pass the filters - plotLocs = plotLocs(toPlot); + contactIdxs = contactIdxs(toPlot); hemiData = hemiData(toPlot, :); channelScoutLabels = channelScoutLabels(toPlot); % Skip plotting if no channels remain after atlas/scout filtering fprintf('\n%s contacts and atlas scout labels:\n', sideName); - if isempty(plotLocs) + if isempty(contactIdxs) if nChannelsBeforeFilter > 0 fprintf('Nothing to plot. All contacts were filtered out by the selected atlas/scout regions.\n'); else @@ -556,10 +537,10 @@ % Print labels and assign colors isAllContactsExcluded = 1; - for i = 1:numel(plotLocs) + for i = 1:numel(contactIdxs) atlasScoutLabelSeeg = channelScoutLabels{i}; - if ~excludedContacts{iSubplot}(plotLocs(i)) - fprintf('%s - %s\n', ChannelMat.Channel(plotLocs(i)).Name, atlasScoutLabelSeeg); + if ~excludedContacts(contactIdxs(i)) + fprintf('%s - %s\n', chanNamesSeeg{contactIdxs(i)}, atlasScoutLabelSeeg); isAllContactsExcluded = 0; end region = GetRegionFromScouts(sCortex, atlasScoutLabelSeeg, OPTIONS); @@ -615,9 +596,9 @@ %% ===== FASTGRAPH TITLE ===== % Build the title shown above each subplot using the stimulation pair and % the atlas label associated with the first contact -function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutLabelsSeeg) +function AddFastgraphTitle(sInput, chanNamesSeeg, atlasScoutLabelsSeeg) % Split the comment into the two parts - parts = strsplit(sInputs(iSortedFastgraph).Comment, '-'); + parts = strsplit(sInput.Comment, '-'); % Clean extracted comment contact1 = strtrim(parts{1}); % Get the contact names @@ -630,7 +611,7 @@ function AddFastgraphTitle(sInputs, iSortedFastgraph, chanNamesSeeg, atlasScoutL else contact1AtlasScoutLabel = '?'; end - title(sprintf('%s\n%s', sInputs(iSortedFastgraph).Comment, contact1AtlasScoutLabel),'fontsize', 8); + title(sprintf('%s\n%s', sInput.Comment, contact1AtlasScoutLabel),'fontsize', 8); end %% ===== RESOLVE SELECTED SCOUTS ===== From bffced1442cc046b4ae69127756ae52a5fd28826 Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 13:44:03 -0400 Subject: [PATCH 49/55] Improve performance on linking axes --- toolbox/process/functions/process_fastgraph.m | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 7952c0e97a..1686aa0f5e 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -217,10 +217,9 @@ gap = [0.075 0.0175]; horzMargin = 0.03; vertMargin = 0.015; - % Generate one FastGraph per selected input - bst_progress('start', 'Process', 'Plotting FastGraphs...', 0, 100); % ===== Get data and Plot each FastGraph and ===== + bst_progress('start', 'Process', 'Plotting FastGraphs...', 0, 100); for iFastGraph = 1:nFastGraphs sInput = sInputs(iFastGraph); stimLoc = stimLocs(iFastGraph, :); @@ -251,9 +250,11 @@ end % === Common feature on FastGraph plots === - % Axes style - axis(hFastGraphAxes, 'tight'); - % Share axes, sets the same XY Limits + % Set limits for axes before linking for improved performance + xlim(hFastGraphAxes, [seegData.Time(1), seegData.Time(end)]*1000); + ylims = ylim(hFastGraphAxes); + ylims = cat(1,ylims{:}); + ylim(hFastGraphAxes, [min(ylims(:,1)), max(ylims(:,2))]); linkaxes(hFastGraphAxes, 'xy'); % Set axis labels xlabel(hFastGraphAxes, 'Time (ms)'); From 39ceb533367aad43834e2e1e73bf626675a140f5 Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 14:18:56 -0400 Subject: [PATCH 50/55] Improve messages in the command window --- toolbox/process/functions/process_fastgraph.m | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 1686aa0f5e..72f3eb3409 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -226,7 +226,7 @@ % Show progress bst_progress('set', round(100 .* (iFastGraph-1) ./ nFastGraphs)); - fprintf('\n===== FastGraph %d/%d: Stimulation site "%s" =====\n', iFastGraph, nFastGraphs, sInput.Comment); + fprintf('\n===== FastGraph %d/%d: Stimulation file "%s" =====\n', iFastGraph, nFastGraphs, sInput.Comment); % Load ONLY SEEG recordings [seegData, excludedContacts] = GetSeegData(sInput, stimLoc, ChannelMat, OPTIONS); @@ -378,9 +378,12 @@ validContacts = ~iExcluded & ~iStimContacts; excludedContacts = ~validContacts; % Report excluded contacts - fprintf('Contacts excluded for being within the %d mm exclusion zone "%s":\n', OPTIONS.ExcludeRadius, sInput.Comment); - fprintf('%s %s ', ChannelMat.Channel(iSeeg(iStimContacts)).Name, ChannelMat.Channel(iSeeg(iExcluded)).Name); - fprintf('\n\n'); + fprintf('Contacts excluded for being at the stimulation location ( <= 2 mm):\n'); + fprintf('%s ', ChannelMat.Channel(iSeeg(iStimContacts)).Name); + fprintf('\n'); + fprintf('Contacts excluded for being within the %d mm exclusion zone:\n', OPTIONS.ExcludeRadius); + fprintf('%s ', ChannelMat.Channel(iSeeg(iExcluded)).Name); + fprintf('\n'); % Set data from excluded channels to NaN seegData.F(excludedContacts, :) = NaN; else @@ -537,11 +540,12 @@ hAreaPlot = area(timeMs, signFactor * hemiData(:, plotWindowIdx)'); % Print labels and assign colors + strMaxLen = max(cellfun(@length, chanNamesSeeg)); isAllContactsExcluded = 1; for i = 1:numel(contactIdxs) atlasScoutLabelSeeg = channelScoutLabels{i}; if ~excludedContacts(contactIdxs(i)) - fprintf('%s - %s\n', chanNamesSeeg{contactIdxs(i)}, atlasScoutLabelSeeg); + fprintf('%-*s - %s\n', strMaxLen, chanNamesSeeg{contactIdxs(i)}, atlasScoutLabelSeeg); isAllContactsExcluded = 0; end region = GetRegionFromScouts(sCortex, atlasScoutLabelSeeg, OPTIONS); @@ -561,6 +565,7 @@ hold off; end end + fprintf('\n'); end %% ===== ATLAS REGION FROM SCOUTS ===== From d6adfbe90574ef0a78fa7028c88dde4bf3968d6f Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 14:55:59 -0400 Subject: [PATCH 51/55] Remove unused input argument --- toolbox/process/functions/process_fastgraph.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 72f3eb3409..0bc21b2273 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -241,7 +241,7 @@ % Create the subplot with custom spacing hFastGraphAxes(iFastGraph) = subtightplot(nRows, nCols, iFastGraph, gap, horzMargin, vertMargin); % Plot the FastGraph for the current stimulation pair - [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactLocIdxs, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); % Apply edge transparency to the subplot set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); @@ -459,7 +459,7 @@ % Create one FastGraph subplot. % Left-hemisphere SEEG channels are plotted as positive stacked areas % Right-hemisphere SEEG channels are plotted as negative stacked areas -function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, ChannelMat, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) +function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) % Initialize output handles hLeftAreaPlot = []; hRightAreaPlot = []; From 21a0d1053fdeeadf63bc0f79e74a92e577a1afa5 Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 6 Aug 2026 15:26:14 -0400 Subject: [PATCH 52/55] Simplify indexing and functions inputs --- toolbox/process/functions/process_fastgraph.m | 79 +++++++++---------- 1 file changed, 38 insertions(+), 41 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 0bc21b2273..0b6da25d73 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -229,19 +229,17 @@ fprintf('\n===== FastGraph %d/%d: Stimulation file "%s" =====\n', iFastGraph, nFastGraphs, sInput.Comment); % Load ONLY SEEG recordings - [seegData, excludedContacts] = GetSeegData(sInput, stimLoc, ChannelMat, OPTIONS); - % Data to be plotted for the current subplot - subplotData = struct(); - % Separate L and R hemisphere data - subplotData.leftData = seegData.F(sContactLocIdxs.Left,:); - subplotData.rightData = seegData.F(sContactLocIdxs.Right,:); - % Sort channels within each hemisphere using the selected metric and time window - sSubplotDataSorted = ApplyDataSorting(subplotData, seegData, OPTIONS); + seegData = GetSeegData(sInput, stimLoc, ChannelMat, OPTIONS); + % Add indices for L and R hemisphere data + seegData.LeftIx = sContactLocIdxs.Left; + seegData.RightIx = sContactLocIdxs.Right; + % Sort L and R indices using the selected metric and time window + seegData = SortHemiIndices(seegData, OPTIONS); % Create the subplot with custom spacing hFastGraphAxes(iFastGraph) = subtightplot(nRows, nCols, iFastGraph, gap, horzMargin, vertMargin); % Plot the FastGraph for the current stimulation pair - [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactLocIdxs, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, seegData, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); % Apply edge transparency to the subplot set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); @@ -356,7 +354,7 @@ %% ===== LOAD AND FILTER SEEG DATA ===== % Load each selected SEEG block and optionally exclude contacts based on % distance from the stimulation site -function [seegData, excludedContacts] = GetSeegData(sInput, stimLoc, ChannelMat, OPTIONS) +function seegData = GetSeegData(sInput, stimLoc, ChannelMat, OPTIONS) % Get index of SEEG channel types iSeeg = channel_find(ChannelMat.Channel, 'SEEG'); % Load current file @@ -376,7 +374,7 @@ iExcluded = (contactDist > 2) & (contactDist <= OPTIONS.ExcludeRadius); % Keep only valid SEEG contacts validContacts = ~iExcluded & ~iStimContacts; - excludedContacts = ~validContacts; + seegData.excludedContacts = ~validContacts; % Report excluded contacts fprintf('Contacts excluded for being at the stimulation location ( <= 2 mm):\n'); fprintf('%s ', ChannelMat.Channel(iSeeg(iStimContacts)).Name); @@ -385,10 +383,10 @@ fprintf('%s ', ChannelMat.Channel(iSeeg(iExcluded)).Name); fprintf('\n'); % Set data from excluded channels to NaN - seegData.F(excludedContacts, :) = NaN; + seegData.F(seegData.excludedContacts, :) = NaN; else % If no stimulation locations are available, keep only SEEG channels - excludedContacts = ~iSeeg; + seegData.excludedContacts = ~iSeeg; end end @@ -415,12 +413,10 @@ sContactGroupLocIdxs = SortLAPRAP(contactLocs); end -%% ===== WITHIN-HEMISPHERE DATA SORTING ===== -% Sort left and right hemisphere channel data within a selected time +%% ===== WITHIN-HEMISPHERE SORTING OF INDICES ===== +% Sort left and right indices for SEEG data within a selected time % window using either RMS amplitude or maximum absolute amplitude -function sSorted = ApplyDataSorting(subplotData, seegData, OPTIONS) - % Initialize output structure - sSorted = struct(); +function seegData = SortHemiIndices(seegData, OPTIONS) % Get sample indices used for sorting if isempty(OPTIONS.SortWindow) sortWindowIdx = 1:size(seegData.F,2); @@ -431,26 +427,30 @@ % Sort channels within each hemisphere using the selected metric switch lower(OPTIONS.SortMethod) case 'rms' - if ~isempty(subplotData.leftData) - leftDataRms = sqrt(sum(subplotData.leftData(:,sortWindowIdx).^2, 2)); + if ~isempty(seegData.LeftIx) + leftDataRms = sqrt(sum(seegData.F(seegData.LeftIx, sortWindowIdx).^2, 2)); leftDataRms(isnan(leftDataRms)) = -Inf; - [sSorted.Vals.Left, sSorted.Idxs.Left] = sort(leftDataRms,'ascend'); + [~, iSort] = sort(leftDataRms(:), 'ascend'); + seegData.LeftIx = seegData.LeftIx(iSort); end - if ~isempty(subplotData.rightData) - rightDataRms = sqrt(sum(subplotData.rightData(:,sortWindowIdx).^2, 2)); + if ~isempty(seegData.RightIx) + rightDataRms = sqrt(sum(seegData.F(seegData.RightIx, sortWindowIdx).^2, 2)); rightDataRms(isnan(rightDataRms)) = -Inf; - [sSorted.Vals.Right, sSorted.Idxs.Right] = sort(rightDataRms,'ascend'); + [~, iSort] = sort(rightDataRms(:), 'ascend'); + seegData.RightIx = seegData.RightIx(iSort); end case 'maxabs' - if ~isempty(subplotData.leftData) - leftDataMax = max(abs(subplotData.leftData(:,sortWindowIdx)),[],2); + if ~isempty(seegData.LeftIx) + leftDataMax = max(abs(seegData.F(seegData.LeftIx, sortWindowIdx)), [], 2); leftDataMax(isnan(leftDataMax)) = -Inf; - [sSorted.Vals.Left, sSorted.Idxs.Left] = sort(leftDataMax,1,'ascend'); + [~, iSort] = sort(leftDataMax(:), 'ascend'); + seegData.LeftIx = seegData.LeftIx(iSort); end - if ~isempty(subplotData.rightData) - rightDataMax = max(abs(subplotData.rightData(:,sortWindowIdx)),[],2); + if ~isempty(seegData.RightIx) + rightDataMax = max(abs(seegData.F(seegData.RightIx, sortWindowIdx)), [], 2); rightDataMax(isnan(rightDataMax)) = -Inf; - [sSorted.Vals.Right, sSorted.Idxs.Right] = sort(rightDataMax,1,'ascend'); + [~, iSort] = sort(rightDataMax(:), 'ascend'); + seegData.RightIx = seegData.RightIx(iSort); end end end @@ -459,7 +459,7 @@ % Create one FastGraph subplot. % Left-hemisphere SEEG channels are plotted as positive stacked areas % Right-hemisphere SEEG channels are plotted as negative stacked areas -function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, subplotData, sSubplotDataSorted, seegData, excludedContacts, sContactGroupLocIdxs, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) +function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, seegData, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) % Initialize output handles hLeftAreaPlot = []; hRightAreaPlot = []; @@ -486,27 +486,24 @@ for iSide = 1:2 if iSide == 1 % Left hemisphere settings - if isempty(subplotData.leftData) + if isempty(seegData.LeftIx) continue; end sideName = 'Left'; - groupLocIdxs = sContactGroupLocIdxs.Left; - sortedIdxs = sSubplotDataSorted.Idxs.Left; + contactIdxs = seegData.LeftIx; signFactor = 1; else % Right hemisphere settings - if isempty(subplotData.rightData) + if isempty(seegData.RightIx) continue; end sideName = 'Right'; - groupLocIdxs = sContactGroupLocIdxs.Right; - sortedIdxs = sSubplotDataSorted.Idxs.Right; + contactIdxs = seegData.RightIx; signFactor = -1; end - % Reorder SEEG channels for the current hemisphere - contactIdxs = groupLocIdxs(sortedIdxs); - hemiData = abs(seegData.F(contactIdxs, :)); + % Get SEEG channels for the current hemisphere + hemiData = abs(seegData.F(contactIdxs, :)); % Get atlas scout labels for these channels channelScoutLabels = cell(1, numel(contactIdxs)); for i = 1:numel(contactIdxs) @@ -544,7 +541,7 @@ isAllContactsExcluded = 1; for i = 1:numel(contactIdxs) atlasScoutLabelSeeg = channelScoutLabels{i}; - if ~excludedContacts(contactIdxs(i)) + if ~seegData.excludedContacts(contactIdxs(i)) fprintf('%-*s - %s\n', strMaxLen, chanNamesSeeg{contactIdxs(i)}, atlasScoutLabelSeeg); isAllContactsExcluded = 0; end From 5b02a5a97ce5ca27f5530ebb7c148b5e24866b95 Mon Sep 17 00:00:00 2001 From: rcassani Date: Wed, 12 Aug 2026 11:47:22 -0400 Subject: [PATCH 53/55] Process: Add `anatparcel` process option --- toolbox/process/panel_process_select.m | 91 ++++++++++++++++++++++++-- 1 file changed, 87 insertions(+), 4 deletions(-) diff --git a/toolbox/process/panel_process_select.m b/toolbox/process/panel_process_select.m index a7eea51f58..79fb33306a 100644 --- a/toolbox/process/panel_process_select.m +++ b/toolbox/process/panel_process_select.m @@ -1299,10 +1299,10 @@ function UpdateProcessOptions() jCombo.setEnabled(isListEnable); % Set current atlas - AtlasSelection_Callback(iProcess, optNames{iOpt}, AtlasList, jCombo, jList, []); + AtlasSelection_Callback(iProcess, optNames{iOpt}, AtlasList, 0, jCombo, jList, []); drawnow; % Set callbacks - java_setcb(jCombo, 'ItemStateChangedCallback', @(h,ev)AtlasSelection_Callback(iProcess, optNames{iOpt}, AtlasList, jCombo, jList, ev)); + java_setcb(jCombo, 'ItemStateChangedCallback', @(h,ev)AtlasSelection_Callback(iProcess, optNames{iOpt}, AtlasList, 0, jCombo, jList, ev)); java_setcb(jList, 'ValueChangedCallback', @(h,ev)ScoutSelection_Callback(iProcess, optNames{iOpt}, AtlasList, jCombo, jList, jCheck, ev)); if ~isempty(jCheck) java_setcb(jCheck, 'ActionPerformedCallback', @(h,ev)ScoutSelection_Callback(iProcess, optNames{iOpt}, AtlasList, jCombo, jList, jCheck, [])); @@ -1629,6 +1629,48 @@ function UpdateProcessOptions() % Set preferred size for the container prefPanelSize = java_scaled('dimension', prefPanelSize(1), prefPanelSize(2)); + case 'anatparcel' + % Get available and selected anatomical parcellations + [AnatAtlasList, iAnatAtlasList] = GetAnatAtlasList(sProcess, optNames{iOpt}); + if isempty(AnatAtlasList) + gui_component('label', jPanelOpt, [], 'Error: No anatomical atlases available.'); + else + % Create list + jList = java_create('javax.swing.JList'); + jList.setLayoutOrientation(jList.HORIZONTAL_WRAP); + jList.setVisibleRowCount(-1); + jList.setCellRenderer(BstStringListRenderer(fontSize)); + % Comment + gui_component('label', jPanelOpt, [], 'Select anatomical parcellations:'); + % Horizontal glue + gui_component('label', jPanelOpt, 'hfill', ' ', [],[],[],[]); + % Atlas selection box + jCombo = gui_component('combobox', jPanelOpt, 'right', [], {AnatAtlasList(:,1)}, [], []); + % Try to re-use previously defined atlas + iDefault = []; + if ~isempty(option.Value) && iscell(option.Value) && (size(option.Value,2) >= 2) && ischar(option.Value{1,1}) + iPrev = find(strcmpi(option.Value{1,1}, AnatAtlasList(:,1))); + if ~isempty(iPrev) + iDefault = iPrev; + end + end + if isempty(iDefault) + iDefault = iAnatAtlasList; + end + % Set current atlas + jCombo.setSelectedIndex(iDefault - 1); + AtlasSelection_Callback(iProcess, optNames{iOpt}, AnatAtlasList, 0, jCombo, jList, []); + drawnow; + % Set callbacks + java_setcb(jCombo, 'ItemStateChangedCallback', @(h,ev)AtlasSelection_Callback(iProcess, optNames{iOpt}, AnatAtlasList, 1, jCombo, jList, ev)); + java_setcb(jList, 'ValueChangedCallback', @(h,ev)ScoutSelection_Callback(iProcess, optNames{iOpt}, AnatAtlasList, jCombo, jList, [], ev)); + % Create scroll panel + jScroll = javax.swing.JScrollPane(jList); + jPanelOpt.add('br hfill vfill', jScroll); + % Set preferred size for the container + prefPanelSize = java_scaled('dimension', 250,180); + end + end jPanelOpt.setPreferredSize(prefPanelSize); end @@ -2239,8 +2281,49 @@ function Cluster_ValueChangedCallback(iProcess, optName, jList, jCheck, ev) end + %% ===== OPTIONS: GET ANATOMY ATLAS LIST ===== + function [AnatAtlasList, iAnatAtlasList] = GetAnatAtlasList(sProcess, optName) + import org.brainstorm.list.*; + % Initialize returned list + AnatAtlasList = {}; + iAnatAtlasList = []; + % Get the current file + if isfield(sProcess.options.(optName), 'InputTypesB') && ~isempty(sFiles2) + curFile = sFiles2(1); + else + curFile = sFiles(1); + end + if isempty(curFile) + return; + end + if isempty(curFile.SubjectFile) + return + end + % Read the subject structure + sSubject = bst_get('Subject', curFile.SubjectFile); + iAnatAtlases = find(cellfun(@(c) ~isempty(strfind(c, '_volatlas')) || ~isempty(strfind(c, 'tissues')), {sSubject.Anatomy.FileName})); + if isempty(iAnatAtlases) + return + end + % Get the names of all parcels in each anatomical atlas + AnatAtlasList = cell(length(iAnatAtlases),2); + for ix = 1 : length(iAnatAtlases) + iAnatAtlas = iAnatAtlases(ix); + AnatAtlasList{ix,1} = sSubject.Anatomy(iAnatAtlas).Comment; + tmp = load(file_fullpath(sSubject.Anatomy(iAnatAtlas).FileName), 'Labels'); + if ~isempty(fields(tmp)) + AnatAtlasList{ix,2} = tmp.Labels(:,2)'; + else + AnatAtlasList{ix,2} = []; + end + end + % Selected anat atlas + iAnatAtlasList = 1; + end + + %% ===== OPTIONS: ATLAS SELECTION CALLBACK ===== - function AtlasSelection_Callback(iProcess, optName, AtlasList, jCombo, jList, ev) + function AtlasSelection_Callback(iProcess, optName, AtlasList, selectAllScouts, jCombo, jList, ev) import org.brainstorm.list.*; % Skip deselected event if ~isempty(ev) && (ev.getStateChange() ~= ev.SELECTED) @@ -2279,7 +2362,7 @@ function AtlasSelection_Callback(iProcess, optName, AtlasList, jCombo, jList, ev end end % If a previous scout selection was not found: select all the scouts - if isempty(iSelScouts) + if isempty(iSelScouts) || selectAllScouts iSelScouts = 1:length(ScoutNames); end % Select scouts in the list From 079c9de5baaf43cfb19267303c3d9303ae98148f Mon Sep 17 00:00:00 2001 From: rcassani Date: Thu, 13 Aug 2026 11:35:10 -0400 Subject: [PATCH 54/55] Store SEEG contact Name and Scout together. Remove unnecessary indexing --- toolbox/process/functions/process_fastgraph.m | 35 +++++++------------ 1 file changed, 12 insertions(+), 23 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 0b6da25d73..0c8d0475aa 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -199,8 +199,7 @@ % Columns whose header matches the atlas name cols = find(any(hit, 1)); % Extract SEEG channel names and their atlas scout labels - chanNamesSeeg = chanTableWithAtlas(2:end, 1); - atlasScoutLabelsSeeg = chanTableWithAtlas(2:end, cols); + chanSeegNameScout = [chanTableWithAtlas(2:end, 1), chanTableWithAtlas(2:end, cols)]; % ===== Create figure for FastGraph ===== hFig = figure; @@ -239,12 +238,12 @@ % Create the subplot with custom spacing hFastGraphAxes(iFastGraph) = subtightplot(nRows, nCols, iFastGraph, gap, horzMargin, vertMargin); % Plot the FastGraph for the current stimulation pair - [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, seegData, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS); + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, seegData, chanSeegNameScout, OPTIONS); % Apply edge transparency to the subplot set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); % Add the stimulation pair and atlas scout label as the subplot title - AddFastgraphTitle(sInput, chanNamesSeeg, atlasScoutLabelsSeeg); + AddFastgraphTitle(sInput, chanSeegNameScout); end % === Common feature on FastGraph plots === @@ -459,7 +458,7 @@ % Create one FastGraph subplot. % Left-hemisphere SEEG channels are plotted as positive stacked areas % Right-hemisphere SEEG channels are plotted as negative stacked areas -function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, seegData, chanNamesSeeg, atlasScoutLabelsSeeg, OPTIONS) +function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, seegData, chanSeegNameScout, OPTIONS) % Initialize output handles hLeftAreaPlot = []; hRightAreaPlot = []; @@ -502,16 +501,9 @@ signFactor = -1; end - % Get SEEG channels for the current hemisphere - hemiData = abs(seegData.F(contactIdxs, :)); - % Get atlas scout labels for these channels - channelScoutLabels = cell(1, numel(contactIdxs)); - for i = 1:numel(contactIdxs) - channelScoutLabels{i} = atlasScoutLabelsSeeg{contactIdxs(i)}; - end % Filter channels using resolved scout selection if hasStimLocs - toPlot = ismember(channelScoutLabels, selectedScoutLabels); + toPlot = ismember(chanSeegNameScout(contactIdxs, 2), selectedScoutLabels); else toPlot = true(1, numel(contactIdxs)); end @@ -519,8 +511,6 @@ nChannelsBeforeFilter = numel(contactIdxs); % Keep only channels that pass the filters contactIdxs = contactIdxs(toPlot); - hemiData = hemiData(toPlot, :); - channelScoutLabels = channelScoutLabels(toPlot); % Skip plotting if no channels remain after atlas/scout filtering fprintf('\n%s contacts and atlas scout labels:\n', sideName); @@ -534,18 +524,17 @@ end % Plot stacked area traces for the current hemisphere - hAreaPlot = area(timeMs, signFactor * hemiData(:, plotWindowIdx)'); + hAreaPlot = area(timeMs, signFactor * abs(seegData.F(contactIdxs, plotWindowIdx))'); % Print labels and assign colors - strMaxLen = max(cellfun(@length, chanNamesSeeg)); + strMaxLen = max(cellfun(@length, chanSeegNameScout(:,1))); isAllContactsExcluded = 1; for i = 1:numel(contactIdxs) - atlasScoutLabelSeeg = channelScoutLabels{i}; if ~seegData.excludedContacts(contactIdxs(i)) - fprintf('%-*s - %s\n', strMaxLen, chanNamesSeeg{contactIdxs(i)}, atlasScoutLabelSeeg); + fprintf('%-*s - %s\n', strMaxLen, chanSeegNameScout{contactIdxs(i), 1}, chanSeegNameScout{contactIdxs(i), 2}); isAllContactsExcluded = 0; end - region = GetRegionFromScouts(sCortex, atlasScoutLabelSeeg, OPTIONS); + region = GetRegionFromScouts(sCortex, chanSeegNameScout{contactIdxs(i), 2}, OPTIONS); hAreaPlot(i).FaceColor = region.Color; end if isAllContactsExcluded @@ -599,7 +588,7 @@ %% ===== FASTGRAPH TITLE ===== % Build the title shown above each subplot using the stimulation pair and % the atlas label associated with the first contact -function AddFastgraphTitle(sInput, chanNamesSeeg, atlasScoutLabelsSeeg) +function AddFastgraphTitle(sInput, chanSeegNameScout) % Split the comment into the two parts parts = strsplit(sInput.Comment, '-'); % Clean extracted comment @@ -608,9 +597,9 @@ function AddFastgraphTitle(sInput, chanNamesSeeg, atlasScoutLabelsSeeg) contact1Parts = strsplit(contact1); contact1 = contact1Parts{end}; % Look up atlas label for the first contact - iContact1 = find(strcmp(chanNamesSeeg, contact1), 1); + iContact1 = find(strcmp(chanSeegNameScout(:,1), contact1), 1); if ~isempty(iContact1) - contact1AtlasScoutLabel = atlasScoutLabelsSeeg{iContact1}; + contact1AtlasScoutLabel = chanSeegNameScout{iContact1, 2}; else contact1AtlasScoutLabel = '?'; end From cf102a2f227066a92ef2400f18ee83fc014dd7e4 Mon Sep 17 00:00:00 2001 From: rcassani Date: Fri, 14 Aug 2026 14:22:45 -0400 Subject: [PATCH 55/55] Refactor to use Anatomical Atlas (volatlas) instead of Surf Atlas --- toolbox/process/functions/process_fastgraph.m | 257 ++++++++++-------- 1 file changed, 149 insertions(+), 108 deletions(-) diff --git a/toolbox/process/functions/process_fastgraph.m b/toolbox/process/functions/process_fastgraph.m index 0c8d0475aa..6be439b805 100644 --- a/toolbox/process/functions/process_fastgraph.m +++ b/toolbox/process/functions/process_fastgraph.m @@ -1,7 +1,7 @@ function varargout = process_fastgraph( varargin ) % PROCESS_FASTGRAPH: Plot FastGraph for one or more SEEG recordings. % For each stimulation pair, channels are split by hemisphere, sorted -% by a user-selected metric, filtered by atlas region or scout label, and +% by a user-selected metric, filtered by anatomical parcels or regions, and % plotted as stacked area plots % % USAGE: @@ -45,18 +45,20 @@ sProcess.OutputTypes = {'data'}; sProcess.nInputs = 1; sProcess.nMinFiles = 1; -% Scouts to use for plotting FastGraph -sProcess.options.scouts.Comment = ''; -sProcess.options.scouts.Type = 'scout'; -sProcess.options.scouts.Value = {}; -% Color FastGraph by Region or by Scout -sProcess.options.colorscheme.Comment = {'Region', 'Scout', 'Color scheme:  '; ... - 'region', 'scout', ''}; +% Anatomical parcels to be use for plotting FastGraph +sProcess.options.parcels.Comment = ''; +sProcess.options.parcels.Type = 'anatparcel'; +sProcess.options.parcels.Value = {}; +% Color FastGraph by Parcel or Region +sProcess.options.colorscheme.Comment = {'Parcel', 'Region', 'Color scheme:  '; ... + 'parcel', 'region', ''}; sProcess.options.colorscheme.Type = 'radio_linelabel'; -sProcess.options.colorscheme.Value = 'region'; +sProcess.options.colorscheme.Value = 'parcel'; sProcess.options.colorscheme.Controller = struct('region', 'region'); % Select regions to include -regionsStr = {'Prefrontal (PF)', 'Frontal (F)', 'Central (C)', 'Parietal (P)', 'Temporal (T)', 'Occipital (O)', 'Limbic (L)'}; +regionsStr = {'Prefrontal(PF)', 'Frontal (F)', 'Central (C)', 'Parietal (P)', ... + 'Temporal (T)', 'Occipital (O)', 'Limbic (L)', 'White', ... + 'CSF', 'Other'}; sProcess.options.region.Comment = [regionsStr, {'Select regions to include:'}]; sProcess.options.region.Type = 'list_horizontal'; sProcess.options.region.Value = regionsStr; @@ -106,13 +108,14 @@ %% ===== GET OPTIONS ===== function OPTIONS = GetOptions(sProcess) OPTIONS = struct(); - % Atlas and scouts to use for plotting FastGraph - OPTIONS.Atlas = sProcess.options.scouts.Value{1,1}; - OPTIONS.AtlasScoutLabels = sProcess.options.scouts.Value{1,2}; + % Anatomy Atlas and Parcels to use for plotting FastGraph + OPTIONS.AnatAtlas = sProcess.options.parcels.Value{1,1}; + OPTIONS.AnatAtlasParcels = sProcess.options.parcels.Value{1,2}; % Color figure by region or by label OPTIONS.ColorScheme = sProcess.options.colorscheme.Value; % Select regions to include - OPTIONS.Region = sProcess.options.region.Value; + OPTIONS.AllRegions = sProcess.options.region.Comment(1:end-1); + OPTIONS.Regions = sProcess.options.region.Value; % Method for sorting the data OPTIONS.SortMethod = sProcess.options.sortmethod.Value; % Time window for sorting the data [s] @@ -142,7 +145,7 @@ OPTIONS = GetOptions(sProcess); % ===== Check regions for 'region' color scheme ===== - if strcmpi(OPTIONS.ColorScheme, 'region') && isempty(OPTIONS.Region) + if strcmpi(OPTIONS.ColorScheme, 'region') && isempty(OPTIONS.Regions) bst_report('Error', sProcess, [], 'No region selected. Select at least one region to run the analysis.'); return; end @@ -192,21 +195,72 @@ sInputs = sInputs(sStimLocIdxs.All); stimLocs = stimLocs(sStimLocIdxs.All, :); + % === Get sAnatAtlas information + sSubject = bst_get('Subject', sInputs(1).SubjectName); + iAnatAtlas = find(strcmp(OPTIONS.AnatAtlas, {sSubject.Anatomy.Comment})); + if isempty(iAnatAtlas) + errMsg = 'TODO: Anat Atlas was not found in Subject'; + return + end + sAnatAtlas = load(file_fullpath(sSubject.Anatomy(iAnatAtlas).FileName), 'Labels'); + % Colors in sAnataAtlas.Labels are in the 0-255 range, convert them to 0-1 range + sAnatAtlas.Labels(:, 3) = cellfun(@(x) x / 255, sAnatAtlas.Labels(:, 3), 'UniformOutput', false); + + % Handle Region color scheme + if strcmpi(OPTIONS.ColorScheme, 'region') + % 1. Get Region for each Parcel in sAnatAtlas. Cortical version of atlas used as reference + [sAnatAtlas, errMsg] = GetParcelRegion(sSubject, OPTIONS.AnatAtlas, sAnatAtlas); + if ~isempty(errMsg) + % Error + end + % 2. Update colors in sAnatAtlas by Parcel region + defaultScoutColors = panel_scout('GetScoutsColorTable'); + allRegionIds = regexprep(OPTIONS.AllRegions, '^.*\((.*?)\).*$', '$1'); + regionColorTable = allRegionIds'; + regionColorTable(1:7, 2) = num2cell(defaultScoutColors(1:7, :), 2); % PF, F, C, P, T, O and L + regionColorTable{ 8, 2} = [220, 220, 220] / 255; % White + regionColorTable{ 9, 2} = [ 44, 152, 254] / 255; % CSF + regionColorTable{ 10, 2} = [130, 130, 130] / 255; % Other + for iRegion = 1 : size(regionColorTable, 1) + iAnatAtlasLabel = ismember(sAnatAtlas.Labels(:,4), regionColorTable{iRegion, 1}); + if any(iAnatAtlasLabel) + [sAnatAtlas.Labels{iAnatAtlasLabel, 3}] = deal(regionColorTable{iRegion, 2}); + end + end + % 3. If no parcel was selected, select parcels from selected regions + if isempty(OPTIONS.AnatAtlasParcels) + regionSelIds = regexprep(OPTIONS.Regions, '^.*\((.*?)\).*$', '$1'); + isKeep = ismember(sAnatAtlas.Labels{:, 4}, regionSelIds); + OPTIONS.AnatAtlasParcels = sAnatAtlas.Labels{isKeep, 1}'; + end + end + % Add Parcel 'N/A': Name and Color + iNA = size(sAnatAtlas.Labels, 1) + 1; + sAnatAtlas.Labels{iNA,2} = 'N/A'; + sAnatAtlas.Labels{iNA,3} = [0,0,0]; + % === Anatomical labels for SEEG contacts - [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 5, 0, 0, OPTIONS.Atlas); - % Locate atlas related columns from channel table above - hit = cellfun(@(x) ischar(x) && (~isempty(strfind(OPTIONS.Atlas, x)) || ~isempty(strfind(x, OPTIONS.Atlas))), chanTableWithAtlas(1,:)); - % Columns whose header matches the atlas name - cols = find(any(hit, 1)); - % Extract SEEG channel names and their atlas scout labels - chanSeegNameScout = [chanTableWithAtlas(2:end, 1), chanTableWithAtlas(2:end, cols)]; + [~, chanTableWithAtlas] = export_channel_atlas(ChannelFiles{1}, 'SEEG', [], 5, 0, 0, OPTIONS.AnatAtlas); + % Locate anatomy atlas in columns from channel table. Full match to get Anatomical parcellations (volatlas) only + iCol = find(strcmp(OPTIONS.AnatAtlas, chanTableWithAtlas(1,:))); + if length(iCol) > 1 + bst_report('Error', sProcess, sInputs, 'Two or more anatomical atlases have the same name, you should rename them to be unique'); + end + % Number of SEEG channels + nSeegChan = size(chanTableWithAtlas,1) - 1; + % SEEG channels: Names, Anatomical Parcel, and Color + seegLocInfo = repmat(struct('Name', '', 'Parcel', '', 'Color', []), nSeegChan, 1); + [seegLocInfo.Name] = deal(chanTableWithAtlas{2:end, 1}); + [seegLocInfo.Parcel] = deal(chanTableWithAtlas{2:end, iCol}); + [~, iAnatAtlasLabel] = ismember({seegLocInfo.Parcel}, sAnatAtlas.Labels(:,2)); + [seegLocInfo.Color] = deal(sAnatAtlas.Labels{iAnatAtlasLabel,3}); % ===== Create figure for FastGraph ===== hFig = figure; hFig.Visible = 'off'; % Maximize figure set(gcf, 'Position', get(0,'Screensize')); - % Reserve one extra subplot for the legend (brain surface with Scouts) + % Reserve one extra subplot for the legend (brain figure) nFastGraphs = length(sInputs); hFastGraphAxes = gobjects(nFastGraphs, 0); % Subplot grid dimensions @@ -238,12 +292,12 @@ % Create the subplot with custom spacing hFastGraphAxes(iFastGraph) = subtightplot(nRows, nCols, iFastGraph, gap, horzMargin, vertMargin); % Plot the FastGraph for the current stimulation pair - [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, seegData, chanSeegNameScout, OPTIONS); + [hLeftAreaPLot, hRightAreaPLot] = PlotFastgraph(sInput, stimLoc, seegData, seegLocInfo, OPTIONS); % Apply edge transparency to the subplot set(hLeftAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); set(hRightAreaPLot,'edgealpha', OPTIONS.EdgeAlpha); - % Add the stimulation pair and atlas scout label as the subplot title - AddFastgraphTitle(sInput, chanSeegNameScout); + % Add the stimulation pair and atlas parcels label as the subplot title + AddFastgraphTitle(sInput, seegLocInfo); end % === Common feature on FastGraph plots === @@ -262,12 +316,12 @@ % === Plot brain legend === bst_progress('set', 100); bst_progress('text', 'Plotting legend...'); - % Generate a cortex snapshot with atlas scout for display - imgCortex = GenerateCortexSnapshot(sInputs, OPTIONS); + % Generate a brain snapshot for display + % imgCortex = GenerateCortexSnapshot(sInputs, OPTIONS); % Create the legend subplot with the same spacing settings axBrain = subtightplot(nRows, nCols, nRows*nCols, gap, horzMargin, vertMargin); % Plot the reference panel with the cortex snapshot and axis labels - PlotLegend(axBrain, imgCortex, round(hFastGraphAxes(1).XLim), hFastGraphAxes(1).YLim); + % PlotLegend(axBrain, imgCortex, round(hFastGraphAxes(1).XLim), hFastGraphAxes(1).YLim); % Close progress bst_progress('stop'); @@ -458,17 +512,12 @@ % Create one FastGraph subplot. % Left-hemisphere SEEG channels are plotted as positive stacked areas % Right-hemisphere SEEG channels are plotted as negative stacked areas -function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, seegData, chanSeegNameScout, OPTIONS) +function [hLeftAreaPlot, hRightAreaPlot] = PlotFastgraph(sInput, stimLoc, seegData, seegLocInfo, OPTIONS) % Initialize output handles hLeftAreaPlot = []; hRightAreaPlot = []; - - % Get cortex to be used for region/color lookup - sSubject = bst_get('Subject', sInput.SubjectName); - CortexFile = sSubject.Surface(sSubject.iCortex).FileName; - sCortex = bst_memory('LoadSurface', CortexFile); - % Resolve selected scouts - selectedScoutLabels = ResolveScoutSelection(sCortex, OPTIONS); + % Selected parcels + selectedParcels = OPTIONS.AnatAtlasParcels; % Check whether stimulation location is available hasStimLocs = any(stimLoc); @@ -501,9 +550,9 @@ signFactor = -1; end - % Filter channels using resolved scout selection + % Filter channels using resolved parcel selection if hasStimLocs - toPlot = ismember(chanSeegNameScout(contactIdxs, 2), selectedScoutLabels); + toPlot = ismember({seegLocInfo(contactIdxs).Parcel}, selectedParcels); else toPlot = true(1, numel(contactIdxs)); end @@ -512,11 +561,11 @@ % Keep only channels that pass the filters contactIdxs = contactIdxs(toPlot); - % Skip plotting if no channels remain after atlas/scout filtering - fprintf('\n%s contacts and atlas scout labels:\n', sideName); + % Skip plotting if no channels remain after atlas/parcel filtering + fprintf('\n%s contacts and anatomical atlas parcels labels:\n', sideName); if isempty(contactIdxs) if nChannelsBeforeFilter > 0 - fprintf('Nothing to plot. All contacts were filtered out by the selected atlas/scout regions.\n'); + fprintf('Nothing to plot. All contacts were filtered out by the selected atlas/parcels regions.\n'); else fprintf('Nothing to plot. No contacts are available for this hemisphere.\n'); end @@ -527,15 +576,14 @@ hAreaPlot = area(timeMs, signFactor * abs(seegData.F(contactIdxs, plotWindowIdx))'); % Print labels and assign colors - strMaxLen = max(cellfun(@length, chanSeegNameScout(:,1))); + strMaxLen = max(cellfun(@length, {seegLocInfo.Name})); isAllContactsExcluded = 1; for i = 1:numel(contactIdxs) if ~seegData.excludedContacts(contactIdxs(i)) - fprintf('%-*s - %s\n', strMaxLen, chanSeegNameScout{contactIdxs(i), 1}, chanSeegNameScout{contactIdxs(i), 2}); + fprintf('%-*s - %s\n', strMaxLen, seegLocInfo(contactIdxs(i)).Name, seegLocInfo(contactIdxs(i)).Parcel); isAllContactsExcluded = 0; end - region = GetRegionFromScouts(sCortex, chanSeegNameScout{contactIdxs(i), 2}, OPTIONS); - hAreaPlot(i).FaceColor = region.Color; + hAreaPlot(i).FaceColor = seegLocInfo(contactIdxs(i)).Color; end if isAllContactsExcluded fprintf('Nothing plotted. All contacts lie within the stimulation-site exclusion zone.\n'); @@ -554,41 +602,10 @@ fprintf('\n'); end -%% ===== ATLAS REGION FROM SCOUTS ===== -% Map an atlas scout label to a Brainstorm region code and plot color -function region = GetRegionFromScouts(sCortex, inputAtlasScoutLabel, OPTIONS) - % Default output if no matching scout is found - region.Name = '?'; - region.Color = [0.5 0.5 0.5]; - % Find the atlas selected by the user - iAtlas = find(strcmpi({sCortex.Atlas.Name}, OPTIONS.Atlas), 1); - if isempty(iAtlas) - return; - end - % Get the selected atlas - atlas = sCortex.Atlas(iAtlas); - % Match the input atlas scout label against atlas scouts - for iScout = 1:numel(atlas.Scouts) - atlasScoutLabel = atlas.Scouts(iScout).Label(1:end-2); - if ~isempty(strfind(lower(inputAtlasScoutLabel), lower(atlasScoutLabel))) - % Matching scout found: assign region name - region.Name = atlas.Scouts(iScout).Region(2:end); - % Assign color based on the selected color scheme - switch lower(OPTIONS.ColorScheme) - case 'region' - region.Color = panel_scout('GetRegionColor', atlas.Scouts(iScout).Region); - case 'scout' - region.Color = atlas.Scouts(iScout).Color; - end - return; - end - end -end - %% ===== FASTGRAPH TITLE ===== % Build the title shown above each subplot using the stimulation pair and % the atlas label associated with the first contact -function AddFastgraphTitle(sInput, chanSeegNameScout) +function AddFastgraphTitle(sInput, seegLocInfo) % Split the comment into the two parts parts = strsplit(sInput.Comment, '-'); % Clean extracted comment @@ -597,43 +614,67 @@ function AddFastgraphTitle(sInput, chanSeegNameScout) contact1Parts = strsplit(contact1); contact1 = contact1Parts{end}; % Look up atlas label for the first contact - iContact1 = find(strcmp(chanSeegNameScout(:,1), contact1), 1); + iContact1 = find(strcmp({seegLocInfo.Name}, contact1), 1); if ~isempty(iContact1) - contact1AtlasScoutLabel = chanSeegNameScout{iContact1, 2}; + contact1AtlasParcelLabel = seegLocInfo(iContact1).Parcel; else - contact1AtlasScoutLabel = '?'; + contact1AtlasParcelLabel = '?'; end - title(sprintf('%s\n%s', sInput.Comment, contact1AtlasScoutLabel),'fontsize', 8); + title(sprintf('%s\n%s', sInput.Comment, contact1AtlasParcelLabel),'fontsize', 8); end -%% ===== RESOLVE SELECTED SCOUTS ===== -% Resolve which atlas scouts should be used based on either: -% 1) explicit scout labels selected by the user, or -% 2) selected anatomical regions from the checkboxes -function [selectedScoutLabels, iSelectedScouts, iAtlas] = ResolveScoutSelection(sCortex, OPTIONS) - % Default outputs - selectedScoutLabels = {}; - iSelectedScouts = []; - iAtlas = []; - % Find selected atlas - iAtlas = find(strcmpi({sCortex.Atlas.Name}, OPTIONS.Atlas), 1); - if isempty(iAtlas) - return; + +%% ===== GET REGION FOR PARCEL ===== +function [sAnatAtlas, errMsg] = GetParcelRegion(sSubject, AnatAtlasName, sAnatAtlas) + errMsg = []; + % 1. Search for cortical version of anatomical atlas to retrieve region labels + sSurf = load(file_fullpath(sSubject.Surface(sSubject.iCortex).FileName), 'Atlas'); + iAnatAtlas = find(strcmp(AnatAtlasName, {sSurf.Atlas.Name})); + if isempty(iAnatAtlas) + errMsg = 'Anatomical atlas does not have cortical version. Needed for regions. Try Parcel colors.'; + return + elseif length(iAnatAtlas) > 1 + errMsg = 'Two or more surface atlases have the same name, you should rename them to be unique'; + return end - atlas = sCortex.Atlas(iAtlas); - if ~isempty(OPTIONS.AtlasScoutLabels) - % Explicit scout-label filtering - isKeep = ismember({atlas.Scouts.Label}, OPTIONS.AtlasScoutLabels); - else - % Region-based filtering - selectedRegions = regexprep(OPTIONS.Region, '^.*\((.*?)\).*$', '$1'); - % Remove the leading character from Brainstorm scout region code - scoutRegions = cellfun(@(x) x(2:end), {atlas.Scouts.Region}, 'UniformOutput', false); - isKeep = ismember(scoutRegions, selectedRegions); + sSurfAtlas = sSurf.Atlas(iAnatAtlas); + % Normalize labels from anatomical atlas and surface atlas + anatAtlasLabels = lower(strrep(sAnatAtlas.Labels(:, 2), ' ', '')); + surfAtlasLabels = lower(strrep({sSurfAtlas.Scouts.Label}, ' ', '')); + + % 2. Obtain region for each parcel in sAnatAtlas + for iAnatAtlasLabel = 1 : length(anatAtlasLabels) + anatAtlasLabel = anatAtlasLabels{iAnatAtlasLabel}; + iScoutFound = find(strcmpi(anatAtlasLabel, surfAtlasLabels)); + % Match + if ~isempty(iScoutFound) + % Region without hemisphere indicator + sAnatAtlas.Labels{iAnatAtlasLabel, 4} = sSurfAtlas.Scouts(iScoutFound(1)).Region(2:end); + continue + end + % Try some common fixes for the label + anatAtlasLabel2 = anatAtlasLabel; + anatAtlasLabel2 = strrep(anatAtlasLabel2, 'antcing', 'anteriorcingulate'); + anatAtlasLabel2 = strrep(anatAtlasLabel2, 'midfront', 'middlefrontal'); + iScoutFound = find(strcmpi(anatAtlasLabel2, surfAtlasLabels)); + if ~isempty(iScoutFound) + % Region without hemisphere indicator + sAnatAtlas.Labels{iAnatAtlasLabel, 4} = sSurfAtlas.Scouts(iScoutFound(1)).Region(2:end); + continue + end + % White matter + if ~isempty(regexp(anatAtlasLabel, '^white[l|r]?$', 'once')) + sAnatAtlas.Labels{iAnatAtlasLabel, 4} = 'White'; + continue + end + % CSF + if strcmpi(anatAtlasLabel, 'csf') + sAnatAtlas.Labels{iAnatAtlasLabel, 4} = 'CSF'; + continue + end + % Other + sAnatAtlas.Labels{iAnatAtlasLabel, 4} = 'Other'; end - % Return selected scout indices and labels - iSelectedScouts = find(isKeep); - selectedScoutLabels = {atlas.Scouts(iSelectedScouts).Label}; end %% ===== GENERATE IMAGE FOR LEGEND =====