-
Notifications
You must be signed in to change notification settings - Fork 2
/
FBWM.m
397 lines (321 loc) · 14.9 KB
/
FBWM.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
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
%% "FBWM": fiber ball white matter model
% input:
% - path_data: path to 4D nifti
% - path_gradient: path to gradient directions in same order as data
% - path_bval: path to bvalues in same order as data
% - degree: degree of spherical harmonics used for fODF
% - path_brain_mask: path to binary mask
% - path_output: path to output data
% - path_dke_parameters= path to template file DKEParameters.txt
% Authors: Emilie McKinnon (mckinnon@musc.edu), August 31, 2018
% Copyright (c) 2018 medical university of south carolina (MUSC)
%
% Permission is hereby granted, free of charge
% to use the Software solely for non-commercial research, including the rights to use
% and modify the Software, subject to the following conditions:
%
% 1. The above copyright notice and this permission notice shall be
% included by Recipient in all copies or substantial portions of the
% Software.
%
% 2. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND,
% EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIESOF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BELIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF ORIN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% 3. In no event shall MUSC be liable for direct, indirect, special,
% incidental or consequential damages in connection with the Software.
% Recipient will defend, indemnify and hold MUSC harmless from any claims or
% liability resulting from the use of the Software by recipient.
%
% 4. Neither anything contained herein nor the delivery of the Software to
% recipient shall be deemed to grant the Recipient any right or licenses
% under any patents or patent application owned by MUSC.
%
% 5. The Software may only be used for non-commercial research and may not
% be used for clinical care.
%
% 6. Any publication by Recipient of research involving the Software shall
% cite the references listed below.
%
%
% REFERENCES
% McKinnon, E. T., Helpern, J. A., & Jensen, J. H. (2018). Modeling white matter microstructure with fiber ball imaging.
% NeuroImage, 176, 11-21.
% Jensen, J. H., Glenn, G. R., & Helpern, J. A. (2016). Fiber ball imaging. Neuroimage, 124, 824-833.
%% input
path_data='/Users/Example/4D_coreg.nii';
path_gradient='/Users/Example/gradient_all.txt';
path_brain_mask='/Users/Example/brain_mask.nii';
degree=6;
path_bval='/Users/Example/bvals_all.txt';
path_output= '/Users/Example/output/';
path_dke_parameters='/Users/Example/DKEParameters.txt';
%%
% read in bvalues and gradient table
bt= textscan(fopen(fullfile(path_bval)),'%f');
gt= textscan(fopen(fullfile(path_gradient)),'%f%f%f');
% define bvalues used
bval=unique(bt{1});
bval_fbi=bval(end);
bval_dki=bval(2:end-1);
index_b1= find(bt{1}==bval_dki(1));
index_b2= find(bt{1}==bval_dki(2));
index_fbi= find(bt{1}==bval_fbi);
% number of gradient directions for each b-value
ndir_b1=length(index_b1) ;
ndir_b2=length(index_b2) ;
ndir_fbi=length(index_fbi);
ndir=[ndir_b1 ndir_b2 ndir_fbi];
% check units of bvalues and put it in ms/um2
if (bval_dki/1000)<1
else
bval_dki=bval_dki/1000;
bval_fbi=bval_fbi/1000;
end
bval=[bval_dki' bval_fbi'];
%Separate b0 from DWIs
idx_b0 = find(bt{1}==0);
% Read in data
S=spm_read_vols(spm_vol(fullfile(path_data)));
dim=size(S);
S_b1=zeros(1,prod(dim(1:3)));
S_b2=zeros(1,prod(dim(1:3)));
Sb6_reshape=zeros(1,prod(dim(1:3)));
S_b1_4D=S(:,:,:,index_b1);
S_b2_4D=S(:,:,:,index_b2);
for i=1:ndir_b1
S_b1(i,:)=reshape(S_b1_4D(:,:,:,i),[1,prod(dim(1:3))]);
end
for i=1:ndir_b2
S_b2(i,:)=reshape(S_b2_4D(:,:,:,i),[1,prod(dim(1:3))]);
end
S_b1_b6=S(:,:,:,index_fbi);
for i=1:length(index_fbi)
Sb6_reshape(i,:)=reshape(S_b1_b6(:,:,:,i),[1,prod(dim(1:3))]);
end
S_b0=S(:,:,:,idx_b0);
S_b0_mean=mean(S_b0,4);
s0=reshape(S_b0_mean,[1, prod(dim(1:3))]);
% brain mask
hdr=spm_vol(path_brain_mask);
brainmask=spm_read_vols(hdr);
brain_mask_reshape=reshape(brainmask,[1,prod(dim(1:3))]);
% seperate gradient tables
GT_fbi = [gt{1}(index_fbi) gt{2}(index_fbi) gt{3}(index_fbi)];
GT_fbi=GT_fbi./sqrt(GT_fbi(:,1).^2+GT_fbi(:,2).^2+GT_fbi(:,3).^2);
GT_b1 = [gt{1}(index_b1) gt{2}(index_b1) gt{3}(index_b1)];
GT_b1=GT_b1./sqrt(GT_b1(:,1).^2+GT_b1(:,2).^2+GT_b1(:,3).^2);
GT_b1=GT_b1(1:ndir_b1,:);
GT_b2 = [gt{1}(index_b2) gt{2}(index_b2) gt{3}(index_b2)];
GT_b2=GT_b2./sqrt(GT_b2(:,1).^2+GT_b2(:,2).^2+GT_b2(:,3).^2);
GT_b2=GT_b2(1:ndir_b2,:);
GT={GT_b1,GT_b2,GT_fbi};
%GET SH BASIS FUNCTIONS
B_fbi=getSH(degree,[ atan2(GT_fbi(:,2),GT_fbi(:,1)) acos(GT_fbi(:,3))],'complex');
B_b1=getSH(degree,[ atan2(GT_b1(:,2),GT_b1(:,1)) acos(GT_b1(:,3))],'complex');
B_b2=getSH(degree,[ atan2(GT_b2(:,2),GT_b2(:,1)) acos(GT_b2(:,3))],'complex');
Harm_id={1,5:9,17:25,37:49,65:81}; % max degree = 8
B_fbi=B_fbi(:,cell2mat(Harm_id(1:degree/2+1)));
B_b1=B_b1(:,cell2mat(Harm_id(1:degree/2+1)));
B_b2=B_b2(:,cell2mat(Harm_id(1:degree/2+1)));
B={B_b1,B_b2,B_fbi};
%% Calculate diffusion tensor
mkdir([path_output '/DKE/'])
S_DKE=cat(4,S_b0_mean,S_b1_4D,S_b2_4D);
hdr.fname=[path_output '/DKE/4D_dki.nii'];
hdr_4D=repmat(hdr,[1,ndir_b1+ndir_b2+1]);
% write the 4D nifti file for DKE
for ii=1:(ndir_b1+ndir_b2+1)
hdr_4D(ii).n=[ii 1];
spm_write_vol(hdr_4D(ii),S_DKE(:,:,:,ii).*brainmask);
end
% write gradient file for DKE
fn = fullfile(path_output,'gradient_b1.txt');
fid = fopen(fn,'w');
fprintf(fid,'%18.15f\t%18.15f\t%18.15f\n',GT_b1');
fclose(fid);
fn = fullfile(path_output,'gradient_b2.txt');
fid = fopen(fn,'w');
fprintf(fid,'%18.15f\t%18.15f\t%18.15f\n',GT_b2');
fclose(fid);
% adjust dkeparameters template file
fid=fopen(path_dke_parameters); % original DKE_parameters
fout=fullfile(path_output,'DKE','DKEParameters.txt'); % new ft_parameters
fidout=fopen(fout,'w');
while(~feof(fid))
s=fgetl(fid);
s=strrep(s,'''/Users/Example/output/gradient.txt''',['{' '''' path_output 'gradient_b1.txt'',''' path_output 'gradient_b2.txt''}']);
s=strrep(s,'/Users/Example/output/',path_output);
s=strrep(s,'ndir = 30',['ndir=[' num2str(ndir_b1) ',' num2str(ndir_b2) ']']);
s=strrep(s,'bval=[0 1000 2000]',['bval=[0' num2str(bval_dki'*1000)]);
fprintf(fidout,'%s\n',s);
disp(s)
end
fclose(fid);
fclose(fidout);
% DKE source code
poolobj = gcp('nocreate');
delete(poolobj);
dke([path_output '/DKE/DKEParameters.txt'])
% DKE command line windows executables
%command=['dke ' path_output '/DKE/DKEParameters.txt'];
%system(command);
load([path_output '/DKE/DT.mat'])
%% Core calculations (can take about 10 min)
% initialisation
tic
De=zeros(3,3,prod(prod(dim(1:3))));
intermediate_A_img=zeros(3,3,prod(prod(dim(1:3))));
Cost=zeros(100,prod(dim(1:3)));
intermediate_DT_img=zeros(3,3,prod(prod(dim(1:3))));
zeta_img=zeros(1,prod(dim(1:3)));
faa_f=zeros(1,prod(dim(1:3)));
% calculate legendre functions and g2l
D0=3;
idx_Y = 0;
for l = 0:degree/2
P2l0(idx_Y+1:idx_Y+(4*l+1), :) = ((-1)^(l) * factorial(2*l))/(4^(l)*factorial(l)^2)*ones(4*l+1,1); %equation 12 FBI paper
G2l(idx_Y+1:idx_Y+(4*l+1), :) = (factorial(l)*(bval_fbi*D0)^((l+0.5))/gamma(2*l+3/2)*hypergeom((l+0.5),2*l+3/2,-bval_fbi*D0))*ones(4*l+1,1); %equation18 FBI paper
idx_Y = idx_Y + 4*l+1;
end
% voxel-wise estimation of other parameters
parfor vox=1:prod(dim(1:3))
G2l_fa_R_b=zeros(length(bval),100, length(cell2mat(Harm_id(1:degree/2+1))));
G2l_fa_R=zeros(length(cell2mat(Harm_id(1:degree/2+1))),100);
G2l_fa_R_large=zeros(length(cell2mat(Harm_id(1:degree/2+1))),100);
%voxelwise signal values
S0_vox=s0(vox);
Sfbi_vox=Sb6_reshape(:,vox);
Sb1_vox=squeeze(S_b1(:,vox));
Sb2_vox=squeeze(S_b2(:,vox));
S={Sb1_vox, Sb2_vox, Sfbi_vox};
if brain_mask_reshape(vox)==1
A2l=(B_fbi'*B_fbi)^-1*B_fbi'*Sfbi_vox/s0(vox); %equation 4 FBWM paper
C2l=A2l*G2l(1).*(sqrt(4*pi)*P2l0*A2l(1).*G2l).^-1; %equation 8 FBWM paper
C2l = bsxfun(@rdivide,C2l,C2l(1)); %normalization
C2l= bsxfun(@times,C2l,1/(sqrt(4*pi))); %normalization
zeta=A2l(1)*sqrt(bval_fbi)./pi;
scale=repmat((((C2l(1)*sqrt(30))).^-1),[1 6]); % scale factor A
A=scale.*[sqrt(30)/3*C2l(1)-sqrt(6)/3*C2l(4)+C2l(6)+C2l(2),...
sqrt(30)/3*C2l(1)-sqrt(6)/3*C2l(4)-C2l(6)-C2l(2),...
sqrt(30)/3*C2l(1)+2*sqrt(6)/3*C2l(4),...
1i*C2l(6)-1i*C2l(2),...
-C2l(5)+C2l(3),...
-1i*C2l(5)-1i*C2l(3)]; %equation 12 FBWM paper
intermediate_DT=[DT(1,vox), DT(4,vox), DT(5,vox);
DT(4,vox), DT(2,vox), DT(6,vox);
DT(5,vox), DT(6,vox), DT(3,vox)];
intermediate_A=[A(1), A(4), A(5);
A(4), A(2), A(6);
A(5), A(6), A(3)];
faa_f(vox)= sqrt(3*sum(abs(C2l(2:6)').^2)./(5*abs(C2l(1)').^2 + 2*sum(abs(C2l(2:6)').^2))); % equation 13 FBWM
awf_grid=linspace(0,1,100);
for b=1:length(bval)
idx_hyper=bval(b)*(awf_grid.^2)* zeta^-2 <20; % when bda becomes large use equation 19 FBI paper
idx_Y = 0;
for l = 0:degree/2
hypergeom_opt=sum((gamma((l+0.5)+(0:100)).*gamma(2*l+3/2).*((-bval(b)*(awf_grid(idx_hyper).^2*zeta^-2))'.^(0:100))./(factorial((0:100)).*gamma(2*l+3/2+(0:100)).*gamma((l+0.5)))),2)';
G2l_fa_R(idx_Y+1:idx_Y+(4*l+1), idx_hyper)=repmat((factorial(l)*(bval(b).*(awf_grid(idx_hyper).^2*zeta^-2)).^(l+0.5)/gamma(2*l+3/2).*hypergeom_opt),4*l+1,1); % equation 9 FBWM paper
idx_Y = idx_Y + 4*l+1;
end
G2l_fa_R_b(b,idx_hyper,:)=G2l_fa_R(:, idx_hyper)';
idx_Y = 0;
for l = 0:degree/2
G2l_fa_R_large(idx_Y+1:idx_Y+(4*l+1), ~idx_hyper)=repmat((exp(-l*(2*l+1)./((2*bval(b).*(awf_grid(~idx_hyper).^2*zeta^-2))))),4*l+1,1); %equation 19 FBI paper
idx_Y = idx_Y + 4*l+1;
end
G2l_fa_R_b(b,~idx_hyper,:)=G2l_fa_R_large(:, ~idx_hyper)';
end
% grid search
for grid=1:100
for b=1:length(bval)
awf=awf_grid(grid);
Se=(S0_vox.*exp((-bval(b)*(1-awf).^-1).*diag((GT{b}*(intermediate_DT-(awf.^3*zeta.^-2)*intermediate_A)*GT{b}'))))*(1-awf); %equation 3 FBWM
Sa=(ones(length(GT{b}),1).*(2*pi*S0_vox*zeta*sqrt(pi./bval(b)))).*(B{b}*(P2l0.*squeeze(G2l_fa_R_b(b,grid,:)).*C2l)); %equation 15 FBWM
Cost(grid,vox)= Cost(grid,vox)+ndir(b)^-1*sum((S{b}-Se-Sa).^2);
end
Cost(grid,vox)=S0_vox^-1*sqrt(length(bval)^-1*Cost(grid,vox)); % equation 21
end
intermediate_DT_img(:,:,vox)=intermediate_DT;
zeta_img(vox)=zeta;
intermediate_A_img(:,:,vox)=intermediate_A;
end
end
toc
%% Calculate other tissue microstructure parameters
De_par=zeros(1,prod(dim(1:3)));
De_rad=zeros(1,prod(dim(1:3)));
De_mean=zeros(1,prod(dim(1:3)));
De_fa=zeros(1,prod(dim(1:3)));
De(isnan(De))=0;
De(isinf(De))=0;
L_1=zeros(1,prod(dim(1:3)));
L_2=zeros(1,prod(dim(1:3)));
L_3=zeros(1,prod(dim(1:3)));
[min_cost,idx]=sort(Cost); % find index where cost function is minimal
awf_grid=linspace(0,1,100);
min_awf=awf_grid(idx(1,:)); % which grid point is this?
Da=min_awf.^2./zeta_img.^2; % Equation 22 FBWM
for i=1:prod(dim(1:3))
if brain_mask_reshape(i)==1
De(:,:,i)=(intermediate_DT_img(:,:,i)-(min_awf(i).^3.*zeta_img(i).^-2).*intermediate_A_img(:,:,i))./(1-min_awf(i)); %equation 19 FBWM
int=De(:,:,i);
int(isnan(De(:,:,i)))=0; %get rid of NaNs
int(isinf(De(:,:,i)))=0; %get rid of INFs
[~, L] = eig(int);
L = diag(L);
L = sort(L, 'descend');
n=1;
while (L(1)<0) || (L(2)<0) || (L(3)<0) % Find different AWF value if the eigenvalues are < 0
n=n+1;
if n<100
min_awf(i)=awf_grid(idx(n,i));
else % if all awf values result in negative eigenvalues set AWF to 0
min_awf(i)=0;
De(:,:,i)=(intermediate_DT_img(:,:,i)-(min_awf(i).^3.*zeta_img(i).^-2).*intermediate_A_img(:,:,i))./(1-min_awf(i));
Da(i)=min_awf(i).^2./zeta_img(i).^2;
break
end
De(:,:,i)=(intermediate_DT_img(:,:,i)-(min_awf(i).^3.*zeta_img(i).^-2).*intermediate_A_img(:,:,i))./(1-min_awf(i));
Da(i)=min_awf(i).^2./zeta_img(i).^2;
end
% calculate eigenvalues again with the correct AWF value
int=De(:,:,i);
int(isnan(De(:,:,i)))=0; %get rid of NaNs
int(isinf(De(:,:,i)))=0; %get rid of INFs
[~, L] = eig(int);
L = diag(L);
L = sort(L, 'descend');
De_par(i)=L(1); % equation 24 FBWM
De_rad(i)=sum(L(2:3))*0.5; %equation 25 FBWM
De_fa(i)=sqrt(((L(1) - L(2)) ^ 2 + (L(1) - L(3)) ^ 2 + (L(2) - L(3)) ^ 2) / (2 * sum(L .^ 2)));
De_mean(i)=1/3*(2*De_rad(i)+De_par(i)); % equation 23 FBWM ( sort of )
L_1(i)=L(1);
L_2(i)=L(2);
L_3(i)=L(3);
end
end
%% write outputs
mkdir([path_output '/FBWM/']);
hdr.dt=[16 0];
zeta_reshape=reshape(zeta_img,dim(1:3));hdr.fname=[path_output '/FBWM/zeta_fbwm.nii'];spm_write_vol(hdr,zeta_reshape);
min_awf_reshape=reshape(min_awf,dim(1:3));hdr.fname=[path_output '/FBWM/awf_fbwm.nii'];spm_write_vol(hdr,min_awf_reshape);
Da_reshape=reshape(Da,dim(1:3));hdr.fname=[path_output '/FBWM/da_fbwm.nii']; spm_write_vol(hdr,Da_reshape);
De_par_reshape=reshape(De_par,dim(1:3));hdr.fname=[path_output '/FBWM/De_par_fbwm.nii']; spm_write_vol(hdr,De_par_reshape);
De_rad_reshape=reshape(De_rad,dim(1:3));hdr.fname=[path_output '/FBWM/De_rad_fbwm.nii']; spm_write_vol(hdr,De_rad_reshape);
De_fa_reshape=reshape(De_fa,dim(1:3));hdr.fname=[path_output '/FBWM/De_fa_fbwm.nii']; spm_write_vol(hdr,De_fa_reshape);
De_mean_reshape=reshape(De_mean,dim(1:3));hdr.fname=[path_output '/FBWM/De_mean_fbwm.nii']; spm_write_vol(hdr,De_mean_reshape);
min_cost_reshape=reshape(min_cost(1,:),dim(1:3));hdr.fname=[path_output '/FBWM/min_cost_fbwm.nii'];spm_write_vol(hdr,min_cost_reshape);
faa_f_reshape=reshape(faa_f,dim(1:3));hdr.fname=[path_output '/FBWM/faa_f_fbwm.nii'];spm_write_vol(hdr,faa_f_reshape);
L_1_reshape=reshape(L_1,dim(1:3));hdr.fname=[path_output '/FBWM/L_1_De_fbwm.nii'];spm_write_vol(hdr,L_1_reshape);
L_2_reshape=reshape(L_2,dim(1:3));hdr.fname=[path_output '/FBWM/L_2_De_fbwm.nii'];spm_write_vol(hdr,L_2_reshape);
L_3_reshape=reshape(L_3,dim(1:3));hdr.fname=[path_output '/FBWM/L_3_De_fbwm.nii'];spm_write_vol(hdr,L_3_reshape);
save([path_output '/FBWM/de.mat'],'De');
save([path_output '/FBWM/A.mat'],'intermediate_A_img');
save([path_output '/FBWM/Cost.mat'],'Cost');
save([path_output '/FBWM/intermediate_DT.mat'],'intermediate_DT_img');
save([path_output '/FBWM/intermediate_A_img.mat'],'intermediate_A_img');