forked from FernandoPalazuelos/Displasias
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdisplasia_tckfixelsample_CSD.m
324 lines (261 loc) · 11.6 KB
/
displasia_tckfixelsample_CSD.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
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
function VALUES = displasia_tckfixelsample(f_tck, f_PDD, f_nComp, ff_values_in, f_prefix)
% VALUES = displasia_tckfixelsample(f_tck, f_PDD, f_nComp, ff_values_in, f_prefix)
%
% f_tck : Filename for the streamlines tck
% f_PDD : Filename for the Principal Diffusion Directions file (MRDS, 4D).
% f_nComp : Filename for the number of components (MRDS, 3D).
% ff_values_in : Cell array of filenames of MRDS metrics to sample.
% Each file should be MRDS, 4D.
% f_prefix : Prefix for the output file names.
%
% Consider:
% addpath('/home/lconcha/software/mrtrix_matlab/matlab');
% addpath(genpath('/home/lconcha/software/dicm2nii-master'))
% addpath /home/lconcha/software/Displasias/
%
% __________________________________________________________________________________
% EXAMPLE:
% f_tck = 'dwi/15/tck/dwi_l_out_resampled_native.tck';
% f_PDD = 'dwi/dwi_MRDS_Diff_BIC_PDDs_CARTESIAN.nii.gz';
% f_MRDS_ncomp = 'dwi/dwi_MRDS_Diff_BIC_NUM_COMP.nii.gz';
% f_MRDS_FA = 'dwi/dwi_MRDS_Diff_BIC_FA.nii.gz';
% f_MRDS_MD = 'dwi/dwi_MRDS_Diff_BIC_MD.nii.gz';
% ff_values = {f_MRDS_FA, f_MRDS_MD};
% f_prefix = '/tmp/prefix';
%
% VALUES = displasia_tckfixelsample(f_tck, f_PDD, f_MRDS_ncomp, ff_values, f_prefix);
% __________________________________________________________________________________
%
% LU15 (0N(H4
% INB-UNAM
% Feb 2023
%% Load tck
% tck_world = read_mrtrix_tracks(f_tck);
% tmptck = [tempname '.tck'];
% fprintf(1, '[INFO] Converting tck to voxel coordinates.\n')
% systemcommand = ['export LD_LIBRARY_PATH="";tckconvert -scanner2voxel ' f_nComp ' ' f_tck ' ' tmptck ' -force -quiet'];
% fprintf(1,' executing: %s\n',systemcommand);
% fprintf('Loading %s\n',f_tck);
% [status,result] = system(systemcommand);
% tck = read_mrtrix_tracks(tmptck);
% [status,result] = system(['rm -f ' tmptck]);
[tck,tck_world] = displasia_load_tck_voxelcoords(f_tck,f_nComp);
%% Load voxel data for PDD and nComp
fprintf('Loading %s\n',f_PDD);
[info,PDD] = displasia_load_nii(f_PDD);
if ndims(PDD) ~= 4
fprintf(1,'ERROR. %s does not have 4 dimensions. It should be a 4D volume with nvolumes = 3, 6, 9 or 12. Bye.\n',f_PDD);
VALUES = [];
return
end
fprintf('Loading %s\n',f_nComp);
[info,nComp] = displasia_load_nii(f_nComp);
if ndims(nComp) ~= 3
fprintf(1,'ERROR. %s should have three dimensions Bye.\n',f_nComp);
VALUES = [];
return
end
%% Prepare tsfs
%nFixels = size(PDD,4) ./ 3;
nFixels = 3; % forcing 3 pixels
tsf_par = tck_world;
tsf_perp = tck_world;
tsf_index_par = tck_world;
tsf_ncomp = tck_world;
tsf_dot_parallel2streamline = tck_world;
tsf_dot_perp2slicenormal = tck_world;
%% Identify parallel/perpendicular
fprintf(1,'Identifying par/perp... ')
for s = 1 : length(tck.data)
if mod(s,10) == 0
fprintf (1,'%d ',length(tck.data)-s);
end
this_streamline = tck.data{s};
this_index_par = zeros(size(this_streamline,1),1);
this_index_perp = zeros(size(this_streamline,1),1);
this_nComp = zeros(size(this_streamline,1),1);
this_dot_parallel2streamline = zeros(size(this_streamline,1),1);
this_dot_perp2slicenormal = zeros(size(this_streamline,1),1);
Rxyz1 = this_streamline(1,:);
origin = [0 0 0];
Rxyz3 = this_streamline(end,:);
PLANE = createPlane(normalizeVector3d(Rxyz1), origin ,normalizeVector3d(Rxyz3)); % create a plane centered at origin
NORMAL = planeNormal(PLANE);
for p = 1 : size(this_streamline,1);
Axyz = this_streamline(p,:);
if p == size(this_streamline,1)
Bxyz = this_streamline(p-1,:);
else
Bxyz = this_streamline(p+1,:);
end
normSegment = (Axyz-Bxyz) ./ norm(Axyz-Bxyz);
mindices = Axyz +1; % matlab offset
PDD1(1) = interp3(PDD(:,:,:,1),mindices(2), mindices(1), mindices(3)); % I cannot get interpn to work, so I do this stupid thing.
PDD1(2) = interp3(PDD(:,:,:,2),mindices(2), mindices(1), mindices(3));
PDD1(3) = interp3(PDD(:,:,:,3),mindices(2), mindices(1), mindices(3));
PDD2(1) = interp3(PDD(:,:,:,4),mindices(2), mindices(1), mindices(3));
PDD2(2) = interp3(PDD(:,:,:,5),mindices(2), mindices(1), mindices(3));
PDD2(3) = interp3(PDD(:,:,:,6),mindices(2), mindices(1), mindices(3));
PDD3(1) = interp3(PDD(:,:,:,7),mindices(2), mindices(1), mindices(3));
PDD3(2) = interp3(PDD(:,:,:,8),mindices(2), mindices(1), mindices(3));
PDD3(3) = interp3(PDD(:,:,:,9),mindices(2), mindices(1), mindices(3));
normPDD1= PDD1./norm(PDD1);
normPDD2= PDD2./norm(PDD2);
normPDD3= PDD3./norm(PDD3);
normPDDs = [normPDD1;normPDD2;normPDD3];
dots(1) = dot(normSegment,normPDD1);
dots(2) = dot(normSegment,normPDD2);
dots(3) = dot(normSegment,normPDD3);
thisnComp = interp3(nComp,mindices(2), mindices(1), mindices(3), 'nearest');
if s==50 & p==4
fprintf('Checking time... ncomp for streamline %d, point %d, is %d\n',s,p,thisnComp);
end
if thisnComp < 3
dots(thisnComp+1:end) = NaN; % Remove PDDs if nCom does not support them.
end
if thisnComp > 1
[themax,indexpar] = max(abs(dots));
[themin,indexperp] = min(abs(dots));
else
[themax,indexpar] = max(abs(dots));
themin = NaN;
indexperp = 3;
end
this_index_par(p,1) = indexpar;
this_index_perp(p,1) = indexperp;
this_nComp(p,1) = thisnComp;
% calculate the absolute dot products between:
% Streamline to parallel tensor
this_dot_parallel2streamline(p,1) = abs(dots(indexpar));
% Slice normal to perpendicular tensor
if thisnComp > 1
this_dot_perp2slicenormal(p,1) = abs(dot(normPDDs(indexperp,:),NORMAL));
else
this_dot_perp2slicenormal(p,1) = -999; % cannot calculate this value if we only found one tensor. -999 is a placeholder for trash.
end
end
try
tsf_index_par.data{s} = this_index_par;
tsf_index_perp.data{s} = this_index_perp;
tsf_ncomp.data{s} = this_nComp;
tsf_dot_parallel2streamline.data{s} = this_dot_parallel2streamline;
tsf_dot_perp2slicenormal.data{s} = this_dot_perp2slicenormal;
VALUES.dot_parallel2streamline(s,:) = this_dot_parallel2streamline;
VALUES.dot_perp2slicenormal(s,:) = this_dot_perp2slicenormal;
VALUES.ncomp{s} = this_nComp;
catch
fprintf(1,'Hey!')
end
end
fprintf (1,'\nFinished identifying par/perp\n',s);
%% Do the sampling
for i = 1 : length(ff_values_in)
f_values_in = ff_values_in{i};
fprintf('Loading %s ... ',f_values_in);
%V = niftiread(f_values_in);
[info,V] = displasia_load_nii(f_values_in);
fprintf(1,'\n');
info = niftiinfo(f_values_in);
[fold,fname,ext] = fileparts(info.Filename);
varName = strrep(fname,'.nii','');
if ndims(V) ~= 4
fprintf(1,'ERROR. %s does not have 4 dimensions. This script can only handle 4D. Bye.\n',f_values_in);
VALUES = [];
return
end
fprintf(1,'[INFO] Sampling %s \n', f_values_in)
for s = 1 : length(tck.data)
if mod(s,10) == 0
fprintf (1,'%d ',length(tck.data)-s);
end
this_streamline = tck.data{s};
this_data_par = zeros(size(this_streamline,1),1);
this_data_perp = zeros(size(this_streamline,1),1);
for p = 1 : size(this_streamline,1);
xyz = this_streamline(p,:);
% vox_indices = [xyz 1] * inv(info.Transform.T);
% vox_indices = vox_indices(1:3);
% mindices = vox_indices + 1;
% matlab_indices = uint8(vox_indices + 1);
mindices = xyz +1;
%thisnComp = interp3(nComp,mindices(2), mindices(1), mindices(3), 'nearest');
thisnComp = tsf_ncomp.data{s}(p);
indexpar = tsf_index_par.data{s}(p);
indexperp = tsf_index_perp.data{s}(p);
vals(1) = interp3(V(:,:,:,1),mindices(2), mindices(1), mindices(3));
vals(2) = interp3(V(:,:,:,2),mindices(2), mindices(1), mindices(3));
vals(3) = interp3(V(:,:,:,3),mindices(2), mindices(1), mindices(3));
if thisnComp < 3
vals(thisnComp+1:end) = -1; % remove values if nComp does not support them.
end
if max(vals) < 0 && thisnComp > 0
fprintf(1,'WTF? All values are invalid!')
fprintf(1,'Streamline %d, point %d\n',s,p);
disp(vals)
end
val_par = vals(indexpar);
val_perp = vals(indexperp);
this_data_par(p,1) = val_par;
this_data_perp(p,1) = val_perp;
if s==3 & p==1
fprintf(1,'checking time\n');
fprintf(1,' Streamline %d, point %d. Parallel value is %1.3g, Perp value is %1.3g\n',s,p,val_par,val_perp);
end
end
tsf_par.data{s} = this_data_par;
tsf_perp.data{s} = this_data_perp;
end
fprintf (1,'\nFinished sampling %s\n',fname);
%%%%% write per-value tsf files
f_tsf_par_out = [f_prefix '_' varName '_par.tsf'];
f_tsf_perp_out = [f_prefix '_' varName '_perp.tsf'];
fprintf(1,' [INFO] Writing tsf_par: %s\n',f_tsf_par_out);
write_mrtrix_tsf(tsf_par,f_tsf_par_out);
fprintf(1,' [INFO] Writing tsf_perp: %s\n',f_tsf_perp_out);
write_mrtrix_tsf(tsf_perp,f_tsf_perp_out);
if regexp(varName,'^[0-9]')
varName = ['x_' varName];
end
VALUES.par.(varName) = tsf_par.data;
VALUES.perp.(varName) = tsf_perp.data;
end
%%%%% writer overall tsf files
fprintf(1,'[INFO] Writing tsf files\n');
f_tsf_dot_parallel2streamline = [f_prefix '_dot_parallel2streamline.tsf'];
f_tsf_dot_perp2slicenormal = [f_prefix '_dot_perp2slicenormal.tsf'];
f_tsf_ncomp = [f_prefix '_ncomp.tsf'];
fprintf(1,' [INFO] Writing tsf_dot_parallel2streamline: %s\n',f_tsf_dot_parallel2streamline);
write_mrtrix_tsf(tsf_dot_parallel2streamline,f_tsf_dot_parallel2streamline);
fprintf(1,' [INFO] Writing tsf_dot_perp2slicenormal: %s\n',f_tsf_dot_perp2slicenormal);
write_mrtrix_tsf(tsf_dot_perp2slicenormal,f_tsf_dot_perp2slicenormal);
fprintf(1,' [INFO] Writing tsf_ncomp: %s\n',f_tsf_ncomp);
write_mrtrix_tsf(tsf_ncomp,f_tsf_ncomp);
f_tsf_par_index_out = [f_prefix '_par_index.tsf'];
fprintf(1,' [INFO] Writing tsf_index_par: %s\n',f_tsf_par_index_out);
write_mrtrix_tsf(tsf_index_par,f_tsf_par_index_out)
fprintf(1,'[INFO] Writing text files\n');
varNames = fieldnames(VALUES.par);
for n = 1 : length(varNames)
thisVarName = varNames{n};
f_txt = [f_prefix '_' thisVarName '_par.txt'];
thismat = cell2mat(VALUES.par.(thisVarName));
fprintf(1,' [INFO] Writing %s\n',f_txt);
save(f_txt,'thismat','-ascii');
f_txt = [f_prefix '_' thisVarName '_perp.txt'];
thismat = cell2mat(VALUES.perp.(thisVarName));
fprintf(1,' [INFO] Writing %s\n',f_txt);
save(f_txt,'thismat','-ascii');
end
f_txt = [f_prefix '_dot_parallel2streamline.txt'];
thismat = VALUES.dot_parallel2streamline;
fprintf(1,' [INFO] Writing %s\n',f_txt);
save(f_txt,'thismat','-ascii');
f_txt = [f_prefix '_dot_perp2slicenormal.txt'];
thismat = VALUES.dot_perp2slicenormal;
fprintf(1,' [INFO] Writing %s\n',f_txt);
save(f_txt,'thismat','-ascii');
f_txt = [f_prefix '_nComp.txt'];
thismat = cell2mat(VALUES.ncomp);
fprintf(1,' [INFO] Writing %s\n',f_txt);
save(f_txt,'thismat','-ascii');