Skip to content

Commit 7fc9368

Browse files
fabriciocravofabriciocravo
committed
Added classical tfce cpp imp
1 parent 110689a commit 7fc9368

File tree

13 files changed

+631
-4
lines changed

13 files changed

+631
-4
lines changed

‎.DS_Store‎

0 Bytes
Binary file not shown.

‎TFCE_Benchmarkin/compare_full_tfce.m‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
% Test script to compare the two TFCE implementations
2+
% Creates a simple symmetric matrix with sequential values for easy debugging
3+
% Each edge will have a unique value incrementing by 0.1
4+
5+
% Create a test matrix (symmetric connectivity matrix)
6+
n = 268; % Full connectivity matrix for 268 ROIs (35,778 edges)
7+
img = zeros(n, n);
8+
9+
% Create random functional connectivity values between -1 and 1
10+
% Fill the upper triangular part of the matrix with random FC values
11+
for i = 1:n
12+
for j = (i+1):n
13+
img(i,j) = (rand() * 2) - 1; % Random value between -1 and 1
14+
end
15+
end
16+
17+
18+
% Make it symmetric
19+
img = img + img';
20+
21+
% Display matrix info
22+
fprintf('Matrix size: %dx%d\n', n, n);
23+
fprintf('Number of unique edges: %d\n', n*(n-1)/2);
24+
fprintf('Max edge value: %.1f\n', max(img(:)));
25+
26+
% Run the first implementation
27+
tfced1 = matlab_tfce_transform(img); % Explicitly set dh=0.10
28+
29+
% Run the second implementation
30+
% [tfced2, comps, comp_sizes] = matlab_tfce_transform(img); % Explicitly set dh=0.10
31+
tfced2 = apply_tfce(img);
32+
33+
% Compare results
34+
disp('Comparison of results:');
35+
disp(['Mean absolute difference: ', num2str(mean(abs(tfced1(:) - tfced2(:))))]);
36+
disp(['Maximum absolute difference: ', num2str(max(abs(tfced1(:) - tfced2(:))))]);
37+
disp(['Correlation between outputs: ', num2str(corr(tfced1(:), tfced2(:)))]);
38+
39+
% Visualize differences
40+
figure('Color', 'white', 'Position', [100, 100, 1200, 1000]);
41+
42+
% Set default font sizes
43+
set(0, 'DefaultAxesFontSize', 14);
44+
set(0, 'DefaultTextFontSize', 16);
45+
46+
% Note: If 'redblue' colormap is not available, you can use:
47+
% - colormap(redblue()) if you have a custom redblue function
48+
% - colormap(diverging_map()) for a custom diverging map
49+
% - colormap([linspace(0,1,128)' linspace(0,0,128)' linspace(1,0,128)'; ...
50+
% linspace(1,1,128)' linspace(0,1,128)' linspace(0,0,128)']) for custom red-white-blue
51+
% - colormap('coolwarm') or colormap('jet') as alternatives
52+
figure('Color', 'white', 'Position', [100, 100, 1200, 1000]);
53+
54+
% Set default font sizes
55+
set(0, 'DefaultAxesFontSize', 14);
56+
set(0, 'DefaultTextFontSize', 16);
57+
58+
subplot(2,2,1);
59+
imagesc(img);
60+
title('Input Matrix', 'FontSize', 18, 'FontWeight', 'bold');
61+
colormap(gca, redblue()); % or use 'coolwarm' if redblue is not available
62+
colorbar('FontSize', 12);
63+
set(gca, 'FontSize', 14);
64+
xlabel('ROI', 'FontSize', 14);
65+
ylabel('ROI', 'FontSize', 14);
66+
67+
subplot(2,2,2);
68+
imagesc(tfced1);
69+
title('Traditional TFCE Output', 'FontSize', 18, 'FontWeight', 'bold');
70+
colormap(gca, 'redblue');
71+
colorbar('FontSize', 12);
72+
set(gca, 'FontSize', 14);
73+
xlabel('ROI', 'FontSize', 14);
74+
ylabel('ROI', 'FontSize', 14);
75+
76+
subplot(2,2,3);
77+
imagesc(tfced2);
78+
title('Incremental Cluster TFCE Output', 'FontSize', 18, 'FontWeight', 'bold');
79+
colormap(gca, 'redblue');
80+
colorbar('FontSize', 12);
81+
set(gca, 'FontSize', 14);
82+
xlabel('ROI', 'FontSize', 14);
83+
ylabel('ROI', 'FontSize', 14);
84+
85+
subplot(2,2,4);
86+
imagesc(abs(tfced1 - tfced2));
87+
title('Absolute Difference', 'FontSize', 18, 'FontWeight', 'bold');
88+
colormap(gca, 'redblue');
89+
colorbar('FontSize', 12);
90+
set(gca, 'FontSize', 14);
91+
xlabel('ROI', 'FontSize', 14);
92+
ylabel('ROI', 'FontSize', 14);

‎TFCE_Benchmarkin/compare_with_simple_edge.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@
2525
%disp('Input Matrix:');
2626
%disp(img);
2727

28-
[tfced1, cc_size] = apply_tfce(img); % Explicitly set dh=0.10
28+
tfced1 = apply_tfce(img); % Explicitly set dh=0.10
2929

3030
% Run the second implementation
3131
% [tfced2, comps, comp_sizes] = matlab_tfce_transform(img); % Explicitly set dh=0.10

‎TFCE_Benchmarkin/redblue.m‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
function cmap = redblue(n)
2+
% REDBLUE Creates a red-white-blue diverging colormap
3+
% CMAP = REDBLUE(N) returns an N-by-3 matrix containing a diverging
4+
% colormap with red for negative values, white for zero, and blue for
5+
% positive values. If N is not specified, N = 256 is used.
6+
%
7+
% Example:
8+
% colormap(redblue)
9+
% colormap(redblue(128))
10+
%
11+
% This colormap is ideal for data that ranges from negative to positive
12+
% values, such as correlation matrices or functional connectivity.
13+
14+
if nargin < 1
15+
n = 256;
16+
end
17+
18+
if n < 2
19+
error('Colormap must have at least 2 colors');
20+
end
21+
22+
% Create the colormap
23+
if mod(n,2) == 0
24+
% Even number of colors
25+
n_half = n/2;
26+
27+
% Red to white (bottom half)
28+
red_to_white = [
29+
linspace(0.7, 1, n_half)' ... % R: dark red to white
30+
linspace(0, 1, n_half)' ... % G: dark red to white
31+
linspace(0, 1, n_half)' % B: dark red to white
32+
];
33+
34+
% White to blue (top half)
35+
white_to_blue = [
36+
linspace(1, 0, n_half)' ... % R: white to dark blue
37+
linspace(1, 0, n_half)' ... % G: white to dark blue
38+
linspace(1, 0.7, n_half)' % B: white to dark blue
39+
];
40+
41+
cmap = [red_to_white; white_to_blue];
42+
else
43+
% Odd number of colors - ensure white is exactly in the middle
44+
n_half = floor(n/2);
45+
46+
% Red to white (bottom half, not including white)
47+
red_to_white = [
48+
linspace(0.7, 1, n_half)' ... % R
49+
linspace(0, 1, n_half)' ... % G
50+
linspace(0, 1, n_half)' % B
51+
];
52+
53+
% White (middle)
54+
white = [1 1 1];
55+
56+
% White to blue (top half, not including white)
57+
white_to_blue = [
58+
linspace(1, 0, n_half)' ... % R
59+
linspace(1, 0, n_half)' ... % G
60+
linspace(1, 0.7, n_half)' % B
61+
];
62+
63+
cmap = [red_to_white; white; white_to_blue];
64+
end
65+
66+
% Ensure values are in [0,1] range
67+
cmap = max(0, min(1, cmap));
68+
end
69+

‎config_files/setparams.m‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@
3131
% Params.data_dir = './data/s_abcd_fc_rosenblatt.mat';
3232
Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
3333
% Params.data_dir = './data/s_hcp_act_noble_1.mat';
34-
Params.output = 'tfce_power_comp';
34+
Params.output = 'tfce_cpp_speed_cmp';
3535

3636
% Save specifications - if NaN output becomes dataset file name
3737
Params.save_directory = './power_calculator_results/';
@@ -81,7 +81,7 @@
8181
Params.save_significance_thresh = 0.15;
8282
% Params.all_cluster_stat_types = {'Parametric', 'Size_cpp', 'Fast_TFCE_cpp', 'Constrained_cpp', 'Omnibus_cNBS'};
8383
Params.all_cluster_stat_types = {'IC_TFCE_FC_cpp_dh1','IC_TFCE_FC_cpp_dh5','IC_TFCE_FC_cpp_dh10', ...
84-
'IC_TFCE_FC_cpp_dh25'};
84+
'IC_TFCE_FC_cpp_dh25', 'TFCE_cpp_dh1', 'TFCE_cpp_dh5', 'TFCE_cpp_dh10', 'TFCE_cpp_dh25'};
8585

8686
Params.all_submethods = {'FWER', 'FDR'};
8787

‎mex_binaries/compile_mex.m‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,8 @@ function compile_mex()
2525
'statistical_methods/mex_scripts/size_pval_cpp.cpp', ...
2626
'statistical_methods/mex_scripts/sparse_size_pval_cpp.cpp', ...
2727
'statistical_methods/mex_scripts/sparse_tfce_cpp.cpp', ...
28-
'statistical_methods/mex_scripts/exact_tfce_cpp.cpp'
28+
'statistical_methods/mex_scripts/exact_tfce_cpp.cpp', ...
29+
'/statistical_methods/mex_scripts/traditional_tfce_cpp.cpp'
2930
};
3031

3132
% Print header

‎statistical_methods/mex_scripts/.DS_Store‎

0 Bytes
Binary file not shown.

‎statistical_methods/mex_scripts/traditional_tfce_cpp.cpp‎

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
#include <vector>
2+
#include <cmath>
3+
#include <algorithm>
4+
#include <queue>
5+
#include <cstring>
6+
#include "mex.h"
7+
8+
// Find connected components using BFS and return cluster edge counts for each node
9+
void findConnectedComponents(const std::vector<std::vector<double>>& adjMatrix,
10+
double threshold,
11+
std::vector<int>& clusterSizePerNode) {
12+
int n = adjMatrix.size();
13+
std::vector<bool> visited(n, false);
14+
std::fill(clusterSizePerNode.begin(), clusterSizePerNode.end(), 0);
15+
16+
for (int start = 0; start < n; ++start) {
17+
if (visited[start]) continue;
18+
19+
// BFS to find component
20+
std::queue<int> q;
21+
std::vector<int> component;
22+
q.push(start);
23+
visited[start] = true;
24+
component.push_back(start);
25+
26+
// Count edges in this component
27+
int edgeCount = 0;
28+
29+
while (!q.empty()) {
30+
int node = q.front();
31+
q.pop();
32+
33+
for (int neighbor = 0; neighbor < n; ++neighbor) {
34+
if (adjMatrix[node][neighbor] >= threshold) {
35+
// Count this edge (only upper triangle to avoid double counting)
36+
if (node < neighbor) {
37+
edgeCount++;
38+
}
39+
40+
if (!visited[neighbor]) {
41+
visited[neighbor] = true;
42+
q.push(neighbor);
43+
component.push_back(neighbor);
44+
}
45+
}
46+
}
47+
}
48+
49+
// Assign cluster size (number of edges) to all nodes in component
50+
for (int node : component) {
51+
clusterSizePerNode[node] = edgeCount;
52+
}
53+
}
54+
}
55+
56+
// Main TFCE computation
57+
std::vector<std::vector<double>> computeTFCE(const std::vector<std::vector<double>>& img,
58+
double H, double E, double dh) {
59+
int n = img.size();
60+
std::vector<std::vector<double>> tfced(n, std::vector<double>(n, 0.0));
61+
62+
// Find maximum value in the matrix
63+
double maxVal = 0.0;
64+
for (int i = 0; i < n; ++i) {
65+
for (int j = i + 1; j < n; ++j) { // Only upper triangle
66+
maxVal = std::max(maxVal, img[i][j]);
67+
}
68+
}
69+
70+
// Generate thresholds - matching Fast_TFCE implementation exactly
71+
std::vector<double> thresholds;
72+
for (double t = 0.0; t <= maxVal + dh; t += dh) {
73+
thresholds.push_back(t);
74+
}
75+
76+
// Vector to store cluster sizes for each node at current threshold
77+
std::vector<int> clusterSizePerNode(n);
78+
79+
// For each threshold (skip index 0, matching Fast_TFCE)
80+
for (size_t h = 1; h < thresholds.size(); ++h) {
81+
double thresh = thresholds[h];
82+
// Find connected components at this threshold
83+
findConnectedComponents(img, thresh, clusterSizePerNode);
84+
85+
// Calculate TFCE contribution for each edge
86+
for (int i = 0; i < n; ++i) {
87+
for (int j = i + 1; j < n; ++j) {
88+
if (img[i][j] >= thresh && clusterSizePerNode[i] > 0) {
89+
// Both nodes should have the same cluster size
90+
double contribution = std::pow(clusterSizePerNode[i], E) *
91+
std::pow(thresh, H) * dh;
92+
93+
tfced[i][j] += contribution;
94+
tfced[j][i] += contribution; // Symmetric
95+
}
96+
}
97+
}
98+
}
99+
100+
return tfced;
101+
}
102+
103+
// MEX gateway function
104+
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
105+
// Check input arguments
106+
if (nrhs != 4) {
107+
mexErrMsgIdAndTxt("TFCE:invalidNumInputs",
108+
"Four inputs required: matrix, H, E, dh");
109+
}
110+
if (nlhs > 1) {
111+
mexErrMsgIdAndTxt("TFCE:invalidNumOutputs",
112+
"Too many output arguments");
113+
}
114+
115+
// Get input matrix
116+
if (!mxIsDouble(prhs[0])) {
117+
mexErrMsgIdAndTxt("TFCE:invalidInput",
118+
"Input matrix must be double");
119+
}
120+
121+
double *matrixPtr = mxGetPr(prhs[0]);
122+
int rows = mxGetM(prhs[0]);
123+
int cols = mxGetN(prhs[0]);
124+
125+
if (rows != cols) {
126+
mexErrMsgIdAndTxt("TFCE:invalidInput",
127+
"Input must be a square matrix");
128+
}
129+
130+
// Get parameters
131+
double H, E, dh;
132+
H = mxGetScalar(prhs[1]);
133+
E = mxGetScalar(prhs[2]);
134+
dh = mxGetScalar(prhs[3]);
135+
136+
// Convert MATLAB matrix to C++ vector (column-major to row-major)
137+
std::vector<std::vector<double>> img(rows, std::vector<double>(cols));
138+
for (int i = 0; i < rows; ++i) {
139+
for (int j = 0; j < cols; ++j) {
140+
img[i][j] = matrixPtr[j * rows + i];
141+
// Apply preprocessing similar to MATLAB version
142+
if (img[i][j] > 1000) {
143+
img[i][j] = 100;
144+
}
145+
// Set diagonal to zero
146+
if (i == j) {
147+
img[i][j] = 0;
148+
}
149+
}
150+
}
151+
152+
// Compute TFCE
153+
std::vector<std::vector<double>> tfced = computeTFCE(img, H, E, dh);
154+
155+
// Create output matrix
156+
plhs[0] = mxCreateDoubleMatrix(rows, cols, mxREAL);
157+
double *outputPtr = mxGetPr(plhs[0]);
158+
159+
// Convert back to MATLAB format (row-major to column-major)
160+
for (int i = 0; i < rows; ++i) {
161+
for (int j = 0; j < cols; ++j) {
162+
outputPtr[j * rows + i] = tfced[i][j];
163+
}
164+
}
165+
}

‎statistical_methods/tfce_dh_methods/.DS_Store‎

2 KB
Binary file not shown.

0 commit comments

Comments
 (0)