Added atlas option for voxels

fabriciocravo · fabriciocravo · commit cf44ebdc7e28 · 2025-09-23T17:02:13.000-04:00
diff --git a/.DS_Store b/.DS_Store
diff --git a/atlas_storage/shen_2mm_268_parcellation__in_subnetworks.nii b/atlas_storage/shen_2mm_268_parcellation__in_subnetworks.nii
diff --git a/atlas_storage/test_hcp_fc_atlas.mat b/atlas_storage/test_hcp_fc_atlas.mat
diff --git a/config_files/setparams.m b/config_files/setparams.m
@@ -29,9 +29,9 @@
 
 % Datasets - Commented for easy use
 % Params.data_dir = './data/s_abcd_fc_rosenblatt.mat';
-Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
-% Params.data_dir = './data/s_hcp_act_noble_1.mat';
-Params.output = 'test_tfce_comparions';
+% Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
+Params.data_dir = './data/s_hcp_act_noble_1.mat';
+Params.output = 'hcp_activation';
 
 % Save specifications - if NaN output becomes dataset file name
 Params.save_directory = './power_calculator_results/';
@@ -56,7 +56,7 @@
 Params.other_scripts_dir='./NBS_benchmarking/support_scripts/';
 
 %%% Resampling parameters %%%
-Params.parallel = true; % run stuff sequentially or in parallel
+Params.parallel = false; % run stuff sequentially or in parallel
 Params.n_workers = 5; % num parallel workers for parfor, best if # workers = # cores
 Params.n_repetitions = 100;  % 500 recommended
 Params.batch_size = 25;
@@ -79,8 +79,8 @@
 Params.pthresh_second_level = 0.05;  % FWER or FDR rate 
 Params.tpr_dthresh = 0; % Threshold for true positives vs negatives
 Params.save_significance_thresh = 0.15;
-Params.all_cluster_stat_types = {'Exact_FC_TFCE', 'Fast_TFCE_dh1', 'Fast_TFCE_dh5', 'Fast_TFCE_dh10'};
 % Params.all_cluster_stat_types = {'Parametric', 'Size_cpp', 'Fast_TFCE_cpp', 'Constrained_cpp', 'Omnibus_cNBS'};
+Params.all_cluster_stat_types = {'Constrained_cpp', 'Omnibus_cNBS'};
 
 Params.all_submethods = {'FWER', 'FDR'};
 
diff --git a/create_flat_function.m b/create_flat_function.m
@@ -5,9 +5,8 @@
     parse(p, varargin{:});
     variable_type = p.Results.variable_type;
     
+    % Ok, it stayed like this because of legacy code, sorry
     switch variable_type
-        case 'node'
-            flat_function = @(x) x(:);
         otherwise
             flat_function = @(x) x(mask);
     end
diff --git a/file_handlers/create_edge_level_stats_in_file.m b/file_handlers/create_edge_level_stats_in_file.m
@@ -1,22 +1,22 @@
-function create_edge_level_stats_in_file(file_path, file_type, n_var, n_repetitions, edge_groups)
+function create_edge_level_stats_in_file(file_path, file_type, n_var, n_repetitions, edge_groups, n_networks)
 
 
     switch file_type
         
         case 'full_file'
 
             edge_level_stats = NaN(n_var, n_repetitions);
-            network_level_stats = NaN(length(unique(edge_groups)) - 1, 1);
+            network_level_stats = NaN(n_networks, 1);
 
             save(file_path, 'edge_level_stats', 'network_level_stats', '-append');
 
         case 'compact_file'
             
             edge_level_stats = zeros(n_var, 1);
-            network_level_stats = zeros(length(unique(edge_groups)) - 1, 1);
+            network_level_stats = zeros(n_networks, 1);
 
             edge_mean_squared_error = zeros(n_var, 1);
-            network_mean_squared_error = zeros(length(unique(edge_groups)) - 1, 1);
+            network_mean_squared_error = zeros(n_networks, 1);
 
             save(file_path, 'edge_level_stats', 'network_level_stats', ...
                 'edge_mean_squared_error', 'network_mean_squared_error', '-append');
diff --git a/file_handlers/extract_relevant_meta_data_from_RP.m b/file_handlers/extract_relevant_meta_data_from_RP.m
@@ -52,6 +52,7 @@
     
     % File related 
     meta_data.atlas_file = RP.atlas_file;
+    meta_data.n_networks = RP.n_networks;
     
     % Sub related data - total subs and this run subs
     meta_data.total_n_subs_1 = RP.n_subs_1;
diff --git a/power_calculator_tools/.DS_Store b/power_calculator_tools/.DS_Store
diff --git a/power_calculator_tools/atlas_data_set_map.m b/power_calculator_tools/atlas_data_set_map.m
@@ -49,7 +49,7 @@ case strcmp(data_set_choice, 'abcd_fc')
                 atlas_file = './atlas_storage/map268_subnetwork.mat';
 
             case strcmp(data_set_choice, 'hcp_act')
-                atlas_file = NaN;
+                atlas_file = './atlas_storage/shen_2mm_268_parcellation__in_subnetworks.nii';
 
             case startsWith(data_set_choice, 'test') && endsWith(data_set_choice, 'fc')
                 atlas_file = './atlas_storage/test_hcp_fc_atlas.mat';
diff --git a/power_calculator_tools/check_calculation_status.m b/power_calculator_tools/check_calculation_status.m
@@ -167,7 +167,7 @@
             meta_data = create_meta_data_file(file_path, meta_data, RP, ids_sampled, existing_repetitions);
 
             create_edge_level_stats_in_file(file_path, RP.subsample_file_type, ...
-                RP.n_var, RP.n_repetitions, RP.edge_groups)
+                RP.n_var, RP.n_repetitions, RP.edge_groups, RP.n_networks)
 
             fprintf('Initialized results file with repetition IDs: %s\n', file_path);
         end
diff --git a/power_calculator_tools/extract_atlas_related_parameters.m b/power_calculator_tools/extract_atlas_related_parameters.m
@@ -1,4 +1,4 @@
-function [n_node_nets, trilmask_net, edge_groups] = extract_atlas_related_parameters(RP, Y)
+function [n_node_nets, trilmask_net, edge_groups, n_networks] = extract_atlas_related_parameters(RP, Y)
 %% extract_atlas_related_parameters
 % **Description**
 % Extracts atlas-related parameters used for network-based statistics 
@@ -42,19 +42,52 @@
     n_node_nets = NaN;
     trilmask_net = NaN;
     edge_groups = [];
+    n_networks = 0;
 
     % Only apply atlas to the network-based stats
     if ~isnan(RP.atlas_file)
-       
-        template_net = summarize_matrix_by_atlas(Y(:, 1), RP.atlas_file, 'suppressimg', 1, 'mask', RP.mask);
-    
-        n_node_nets = size(template_net, 1); % square
-        % trilmask_net = tril(true(n_node_nets));
-    
-        edge_groups = load_atlas_edge_groups(RP.atlas_file);
-        % edge_groups = tril(edge_groups,-1);
-        edge_groups = triu(edge_groups, 1);
+
+        [~, ~, atlas_ext] = fileparts(RP.atlas_file);
+
+
+        switch atlas_ext
+
+            case '.nii'
+                [edge_groups, n_node_nets] = nii_case_atlas_extraction(RP.atlas_file);
+
+            case '.mat'
+                [edge_groups, n_node_nets] = mat_case_atlas_extraction(Y(:, 1),  RP.atlas_file, RP.mask);
+
+            otherwise 
+                error('Atlas extension not supported')
+
+        end
         
+        
+        unique_vals = unique(edge_groups);
+        n_networks = sum(unique_vals > 0);      
     end
 
-end
+
+end
+
+
+function [edge_groups, n_nodes_nets] = mat_case_atlas_extraction(Y_exp, atlas_file, mask)
+    
+    template_net = summarize_matrix_by_atlas(Y_exp, atlas_file, 'suppressimg', 1, 'mask', mask);
+    
+    n_nodes_nets = size(template_net, 1); % square
+    % trilmask_net = tril(true(n_node_nets));
+
+    edge_groups = load_atlas_edge_groups(atlas_file);
+    edge_groups = triu(edge_groups, 1);
+
+end
+
+function [edge_groups, n_nodes_nets] = nii_case_atlas_extraction(atlas_file)
+    
+    edge_groups = niftiread(atlas_file);
+    n_nodes_nets = max(edge_groups(:));
+
+end
+
diff --git a/power_calculator_tools/flatten_edge_groups.m b/power_calculator_tools/flatten_edge_groups.m
diff --git a/power_calculator_tools/set_up_nbs_parameters.m b/power_calculator_tools/set_up_nbs_parameters.m
@@ -199,6 +199,8 @@
         % If no edge groups are provided and it's not preaveraged, disable use_edge_groups
         nbs.STATS.use_edge_groups = 0;
     end
+    % Force it to pass, the check is not done at this level anymore
+    UI.edge_groups.ok = 1;
     
     % Number of nodes
     nbs.STATS.N = DIMS.nodes; 
diff --git a/run_benchmarking.m b/run_benchmarking.m
@@ -42,7 +42,8 @@ function run_benchmarking(RP, Y, X)
         RP.n_subs_subset = RP.list_of_nsubset{id_nsub_list};
         RP = set_n_subs_subset(RP);
 
-        [RP.node_nets, RP.trilmask_net, RP.edge_groups] = extract_atlas_related_parameters(RP, Y);
+        [RP.node_nets, RP.trilmask_net, RP.edge_groups, RP.n_networks] = ...
+            extract_atlas_related_parameters(RP, Y);
         
         % Prepare benchmarking setup
         [UI, RP] = setup_benchmarking(RP);
diff --git a/statistical_methods/Constrained_cpp.m b/statistical_methods/Constrained_cpp.m
@@ -24,12 +24,12 @@
             edge_stats = params.edge_stats;
             permuted_edge_stats = params.permuted_edge_data;
 
-            flat_edge_groups = flatten_edge_groups(STATS.edge_groups, STATS.mask);
+            flat_edge_groups = double(STATS.flatten_matrix(STATS.edge_groups));
             
             pvals = struct();
             [p_FWER, p_FDR] = constrained_pval_cpp(edge_stats, permuted_edge_stats, ...
                 flat_edge_groups, STATS.alpha);
-            
+
             if STATS.submethods.FWER
                 pvals.FWER = p_FWER;
             end
diff --git a/statistical_methods/Omnibus_cNBS.m b/statistical_methods/Omnibus_cNBS.m
@@ -56,6 +56,10 @@
     % Aggregate edge statistics by network
     for i = 1:length(edge_stats)
         network_index = flatted_edge_groups(i);
+        if network_index == 0
+            continue;
+        end
+
         return_vector(network_index) = return_vector(network_index) + edge_stats(i);
     end
 end
diff --git a/statistical_methods/mex_scripts/constrained_pval_cpp.cpp b/statistical_methods/mex_scripts/constrained_pval_cpp.cpp
@@ -119,6 +119,10 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
     int max_network_idx = 0;
     for (mwSize i = 0; i < num_edges; i++) {
         int idx = static_cast<int>(network_indices[i]);
+        
+        // Zero is not a network
+        if(idx == 0) continue;
+
         network_set.insert(idx);
         if (idx > max_network_idx) {
             max_network_idx = idx;
@@ -137,6 +141,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
     // Calculate network statistics for observed data
     for (mwSize e = 0; e < num_edges; e++) {
         int network_idx = static_cast<int>(network_indices[e]);
+        
+        if(network_idx == 0) continue;
+
         network_stats[network_idx] += edge_stats[e];
     }
     
@@ -154,6 +161,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
         // Calculate network statistics for this permutation
         for (mwSize e = 0; e < num_edges; e++) {
             int network_idx = static_cast<int>(network_indices[e]);
+            
+            if(network_idx == 0) continue;
+
             perm_network_stats[network_idx] += perm_data[e];
         }
         
