/
trk_compile_data.m
215 lines (188 loc) · 9.81 KB
/
trk_compile_data.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
function [track_means,starting_pts_out,nPts] = trk_compile_data(subsDir,subIDs,tract_info,outDir,starting_pts_in, saveTrk, saveASCII)
%TRK_COMPILE_DATA - Compiles along-tract data for subjects/hemispheres/tracts
%
% Syntax: [track_means,starting_pts_out,nPts] = trk_compile_data(subsDir,subIDs,tract_info,outDir,starting_pts_in,saveTrk,saveASCII)
%
% Inputs:
% subsDir - Path to subject directory [char]
% subIDs - List of subject ID folders in subsDir [nSubs x 1 cell or num array]
% tract_info - Dataset with tract names, default origins, and viewpoints
% outDir - Path to output directory. (Default: PWD)
% starting_pts_in - Dataset with tract origins. Useful for reusing past
% origins that were determined interactively
% saveTrk - Save the mean tract geometries to .trk files? [logical]
% (Default: 0)
% saveASCII - Save raw streamlines to ASCII files? [logical] (Default: 0)
%
% Outputs:
% track_means - Structure array with mean tract geometries and along-tract
% scalar mean estimates.
% starting_pts_out - Dataset with tract origins used for flipping. Possibly a
% combination of defaults stored in the tract_info dataset, and
% interactively chosen ones.
% nPts - The number of interpolation points used for each tract
% and subject [nTracts x nSubjects]
%
% Output files:
% starting_pts_out.txt - Same as above
% Tracking_QC_*.pdf - Quality control images for each tract. Allows you to
% quickly see if the flipping/interpolating appears successful.
% trk_data.txt - Data table with scalar means and standard deviations for
% each subject/tract/hemisphere (for stats in R).
% trk_props_long.txt - Data table with number of streamlines for each
% subject/tract/hemisphere (for stats in R).
% <subID>_<trkName>.txt ? (Optionally) Raw streamlines from the original
% tract group in a plain text ASCII format for easy plotting in R.
% http://github.com/johncolby/along-tract-stats/wiki/single-subject
% <subID>_<trkName>_mean.trk ? (Optionally) Mean tract geometry with attached
% along-tract cross-sectional mean scalar metric (e.g. FA) for display in
% TrackVis.
% http://www.colbyimaging.com/wiki/neuroimaging/along-tract-stats#trackvis_visualization
%
%
% Example:
% exDir = '/path/to/along-tract-stats/example';
% subIDs = {'subject1'};
% tract_info = dataset('file', fullfile(exDir, 'tract_info.txt'));
% starting_pts_in = dataset('file', fullfile(exDir, 'starting_pts_out.txt'));
% [track_means,starting_pts_out,nPts] = trk_compile_data(exDir,subIDs,tract_info,[],starting_pts_in,1,1);
%
% Other m-files required: read_avw, trk_read, trk_interp, trk_flip, trk_restruc,
% trk_add_sc, trk_mean_sc, trk_plot, trk_write_ascii, dataset, export
% Subfunctions: none
% MAT-files required: none
%
% See also: DATASET
% Author: John Colby (johncolby@ucla.edu)
% UCLA Developmental Cognitive Neuroimaging Group (Sowell Lab)
% Sept 2010
%% Check and format input arguments
if nargin < 7 || isempty(saveASCII), saveASCII = 0; end
if nargin < 6 || isempty(saveTrk), saveTrk = 0; end
if nargin < 5 || isempty(starting_pts_in)
starting_pts_in = dataset();
elseif ~iscell(starting_pts_in.Subject) % Reformat to cells if subIDs are all numeric
starting_pts_in.Subject = cellfun(@num2str, num2cell(starting_pts_in.Subject), 'UniformOutput', 0);
end
if nargin < 4 || isempty(outDir), outDir = pwd; end
if nargin < 3, error('Need more input arguments'), end
if isnumeric(subIDs), subIDs = cellstr(num2str(subIDs)); end
tract_info.startPt = cell2mat(cellfun(@eval, tract_info.startPt, 'UniformOutput',false));
tract_info.view = cell2mat(cellfun(@eval, tract_info.view, 'UniformOutput',false));
%% Setup output files and variables
% Save one file for overall track properties
fid1 = fopen(fullfile(outDir, 'trk_props_long.txt'), 'wt');
fprintf(fid1, 'ID\tHemisphere\tTract\tStreamlines');
% Save another file with the mean FA +/- SD at each point along the track
fid2 = fopen(fullfile(outDir, 'trk_data.txt'), 'wt');
fprintf(fid2, 'ID\tPoint\tHemisphere\tTract\tFA\tSD');
% Initialize variables
track_means = struct([]);
starting_pts_out = dataset();
nPts = zeros(length(tract_info),length(subIDs));
%% Main loop to extract along-tract properties
% Loop over tracks
for iTrk=1:length(tract_info)
ipage = 1;
newpage = 0;
fh = figure; hold on
% Loop over subjects
for i=1:length(subIDs)
pt_start = [];
if (i-1)/25 == ipage, newpage = 1; end
try
% Load scalar volume
% Note: Modify path according to your directory setup
subStr = subIDs{i};
trkName = sprintf('%s_%s', tract_info.Tract{iTrk}, tract_info.Hemisphere{iTrk});
volPath = fullfile(subsDir, subStr, 'dti_fa.nii.gz');
volume = read_avw(volPath);
% Load tract group
% Note: Modify path according to your directory setup
trkPath = fullfile(subsDir, subStr, strcat(trkName, '.trk'));
[header tracks] = trk_read(trkPath);
% Determine # of interpolation points
if isnan(tract_info.nPts(iTrk))
nPts(iTrk,i) = round(mean(trk_length(tracks))/header.voxel_size(1));
else
nPts(iTrk,i) = tract_info.nPts(iTrk);
end
% Interpolate streamlines
tracks_interp = trk_interp(tracks, nPts(iTrk,i), [], 1);
nPts(iTrk,i) = size(tracks_interp, 1);
% Determine 'pt_start' near tract origin. First look in
% 'starting_pts_in' if available, and then the 'tract_info' defaults
if ~isempty(starting_pts_in)
pt_start = double(starting_pts_in(strcmp(subStr, starting_pts_in.Subject) &...
strcmp(tract_info.Tract(iTrk), starting_pts_in.Tract) &...
strcmp(tract_info.Hemisphere(iTrk), starting_pts_in.Hemisphere),4:6));
end
if isempty(pt_start)
pt_start = tract_info.startPt(iTrk,:);
end
% Reorient streamlines according to 'pt_start'. Determine
% interactively if needed
[tracks_interp pt_start] = trk_flip(header, tracks_interp, pt_start, volume);
starting_pts_out = [starting_pts_out; dataset({subStr}, tract_info.Hemisphere(iTrk), tract_info.Tract(iTrk), pt_start(1), pt_start(2), pt_start(3),...
'VarNames',{'Subject','Hemisphere','Tract','PointX','PointY','PointZ'})];
tracks_interp_str = trk_restruc(tracks_interp);
% Extract scalar values from 'volume'
[header_sc tracks_sc] = trk_add_sc(header,tracks_interp_str,volume,'FA');
% Determine the mean scalar at each cross section along the tract group
[scalar_mean scalar_sd] = trk_mean_sc(header_sc,tracks_sc);
% Write outputs
fprintf(fid1, '\n%s\t%s\t%s\t%d', subStr, tract_info.Hemisphere{iTrk}, tract_info.Tract{iTrk}, header.n_count);
for iPt=1:nPts(iTrk,i)
fprintf(fid2, '\n%s\t%d\t%s\t%s\t%0.4f\t%0.4f', subStr, iPt, tract_info.Hemisphere{iTrk}, tract_info.Tract{iTrk}, scalar_mean(iPt), scalar_sd(iPt));
end
% Determine the mean streamline geometry for display in QC figures
track_mean = mean(tracks_interp, 3);
track_mean_sc = [track_mean scalar_mean];
track_means = [track_means struct('Subject', subStr, 'Tract', tract_info.Tract{iTrk}, 'Hemisphere', tract_info.Hemisphere{iTrk}, 'track_mean_sc_str', trk_restruc(track_mean_sc))];
% Open a new QC figure if needed
if newpage
set(gcf, 'PaperSize', [10.5 8])
set(gcf, 'PaperPosition', [0 0 10.5 8])
print(gcf, '-dpdf', fullfile(outDir, sprintf('Tracking_QC_%s_%d.pdf', trkName, ipage)), '-r300')
close(fh)
fh = figure; hold on
ipage = ipage+1;
newpage = 0;
end
% Draw QC figure
figure(fh)
subplot(5,5,i-25*(ipage-1))
trk_plot(header, trk_restruc(track_mean_sc), volume, [])
view(tract_info.view(iTrk,:))
title(subStr)
axis off
% Save the mean tract geometry if desired
if saveTrk
% x y z sc1
track_mean_sc(1:2,:,2) = [0 0 0 0; %min
0 0 0.1 1]; %max
track_mean_sc_str = trk_restruc(track_mean_sc);
header_mean_sc = header_sc;
header_mean_sc.n_count = 2;
trk_write(header_mean_sc, track_mean_sc_str, fullfile(outDir, sprintf('%s_%s_mean.trk', subStr, trkName)))
end
% Save raw streamlines to ASCII if desired
if saveASCII
tracks_sc_mat = trk_restruc(tracks_sc);
trk_write_ascii(tracks_sc_mat, fullfile(outDir, sprintf('%s_%s.txt', subStr, trkName)))
end
catch me % No streamlines
fprintf(fid1, '\n%s\t%s\t%s\t0', subStr, tract_info.Hemisphere{iTrk}, tract_info.Tract{iTrk});
fprintf('Failed to process subject %s %s\n', subStr, trkName)
warning(me.message)
end
end
% Save QC figure
set(gcf, 'PaperSize', [10.5 8])
set(gcf, 'PaperPosition', [0 0 10.5 8])
print(gcf, '-dpdf', fullfile(outDir, sprintf('Tracking_QC_%s_%d.pdf', trkName, ipage)), '-r300')
close(fh)
end
%% Clean up
export(starting_pts_out, 'file', fullfile(outDir, 'starting_pts_out.txt'))
fclose all;