forked from hMRI-group/hMRI-toolbox
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhmri_create_MTProt.m
1628 lines (1445 loc) · 72.3 KB
/
hmri_create_MTProt.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
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [fR1, fR2s, fMT, fA, PPDw, PT1w, PMTw] = hmri_create_MTProt(jobsubj) %#ok<*STOUT>
%==========================================================================
% This is hmri_create_MTProt, part of the hMRI-Toolbox.
%
% PURPOSE
% Estimation of quantitative parameters (R1, R2*, PD, MT) from
% multi-contrast multi-echo FLASH protocol
%
% FORMAT
% [fR1, fR2s, fMT, fA, PPDw, PT1w, PMTw] = hmri_create_MTProt(jobsubj)
%
% INPUTS
% jobsubj parameters for one subject out of the job list.
% NB: ONE SINGLE DATA SET FROM ONE SINGLE SUBJECT IS
% PROCESSED HERE, LOOP OVER SUBJECTS DONE AT HIGHER LEVEL.
% P_trans filenames of a pair of magnitude image + transmit bias map
% [p.u.] of B1+ field (not mandatory), saved as
% "b1_trans_input" field in jobsubj.
%
% OUTPUTS
% fR1 R1 map output filename (only quantitative if B1+ map provided)
% fR2s R2* map output filename (simple fit on PD data or OLS across
% all contrasts, according to defaults settings)
% fMT Magnetization transfer (MT) map output filename
% fA Proton density map output filename (free water concentration
% (PD) [%] or signal amplitude (A) [a.u.] depending on job and
% defaults settings (PDproc & RF sensitivity bias correction)).
% PPDw average PD-weighted image filename (or OLS fit at TE=0 if fullOLS = true)
% PT1w average T1-weighted image filename (or OLS fit at TE=0 if fullOLS = true)
% PMTw average MT-weighted image filename (or OLS fit at TE=0 if fullOLS = true)
%
% OTHER USEFUL VARIABLES EXPLAINED
% P_mtw, P_pdw, P_t1w (from jobsubj.raw_mpm) are MTw, PDw, T1w series of
% images respectively (multiple echoes)
% TE_mtw, TE_pdw, TE_t1w: echo times of the first echo (volume) of each
% contrast, respectively.
% TR_mtw, TR_pdw, TR_t1w: repetition time for each contrast,
% respectively.
% fa_mtw, fa_pdw, fa_t1w: excitation flip angles for each contrast,
% respectively.
%
%==========================================================================
% FEATURES AND REFERENCES
% Estimation of R1 and MT maps
% Helms et al., Magnetic Resonance in Medicine 60:1396-1407 (2008)
% Helms et al., Magnetic Resonance in Medicine 59:667-672 (2008)
%
% B1 correction of MT maps
% Weiskopf et al., Front. Neurosci. 2013, doi:
% 10.3389/fnins.2013.00095
%
% Imperfect spoiling correction
% Preibisch and Deichmann, MRM 61:125-135 (2009)
% Corrects for residual transverse magnetization coherences that
% remain despite RF & gradient spoiling and result in deviations from
% the Ernst equation (ideal expected signal in a FLASH acquisition).
% Modelling data and P2_a, P2_b correction factors calculation based
% on code supplied by Ralf Deichmann. Only a small subset of
% experimental parameters have been used and imperfect spoiling
% correction can only be applied if these correction factors are
% defined specifically for the acquisition protocol used.
%
% OLS R2* map calculation
% Ordinary least square fit of R2* estimated from all echoes from all
% three contrasts instead of the PD-weighted image only. This
% increases the SNR in R2* maps. It is noted that it potentially
% mixes in additional MT and T1 contrast, as discussed in:
% Weiskopf et al. Front. Neurosci. 2014 DOI: 10.3389/fnins.2014.00278
% This reference should be cited if you use this output.
%
%==========================================================================
% Written by
% Gunther Helms, MR-Research in Neurology and Psychiatry, University
% of Goettingen
% Martina Callaghan, John Ashburner, Wellcome Trust Centre for
% Neuroimaging at UCL, London
% Antoine Lutti, CHUV, Lausanne, Switzerland
% Nikolaus Weiskopf, Department of Neurophysics, Max Planck Institute
% for Human Cognitive and Brain Sciences, Leipzig, Germany
%
%==========================================================================
%% =======================================================================%
% Initiate processing
%=========================================================================%
flags = jobsubj.log.flags;
flags.PopUp = false;
hmri_log(sprintf('\t============ MPM PROCESSING - %s.m (%s) ============', mfilename, datestr(now)),flags);
% retrieves all required parameters for MPM processing
mpm_params = get_mpm_params(jobsubj);
% for convenience, define a few parameters to make formulae more readable
% and avoid number of repetitions:
% B1 transmit input:
% P_trans(1,:) = magnitude image (anatomical reference for coregistration)
% P_trans(2,:) = B1 map (p.u.)
P_trans = jobsubj.b1_trans_input;
% index number for each contrast - zero index means no images available
PDwidx = mpm_params.PDwidx;
T1widx = mpm_params.T1widx;
MTwidx = mpm_params.MTwidx;
% TE/TR/FA for each contrast
% T1w and MTw contrasts are optional.
% PDw must be present or nothing can be calculated. The script will abort
% at the next line if no PDw echoes available. In that case, a warning has
% been thrown earlier (in get_mpm_params).
TE_pdw = mpm_params.input(PDwidx).TE;
TR_pdw = mpm_params.input(PDwidx).TR;
fa_pdw = mpm_params.input(PDwidx).fa;
if T1widx
TE_t1w = mpm_params.input(T1widx).TE; %#ok<NASGU>
TR_t1w = mpm_params.input(T1widx).TR;
fa_t1w = mpm_params.input(T1widx).fa;
end
if MTwidx
TE_mtw = mpm_params.input(MTwidx).TE; %#ok<NASGU>
TR_mtw = mpm_params.input(MTwidx).TR;
fa_mtw = mpm_params.input(MTwidx).fa;
end
% other parameters
threshall = mpm_params.proc.threshall;
PDproc = mpm_params.proc.PD;
ISC = mpm_params.proc.ISC;
dt = [spm_type('float32'),spm_platform('bigend')]; % for nifti output
outbasename = spm_file(mpm_params.input(PDwidx).fnam(1,:),'basename'); % for all output files
calcpath = mpm_params.calcpath;
mpm_params.outbasename = outbasename;
respath = mpm_params.respath;
supplpath = mpm_params.supplpath;
% Number of echoes averaged (maximum number or echoes available for ALL
% contrasts AND under TE_limit (+1) - see get_mpm_params)
avg_nr = mpm_params.nr_echoes4avg;
% load PDw images
% PDw images are the reference space for all results. Therefore, the matrix
% dimension defined below (dm) is used across the whole script. It must not
% been redefined.
V_pdw = spm_vol(mpm_params.input(PDwidx).fnam);
dm = V_pdw(1).dim;
%% =======================================================================%
% Calculate R2* map from PDw echoes
%=========================================================================%
fR2s = '';
hmri_log(sprintf('\t-------- R2* map calculation (basic exponential decay) --------'),mpm_params.nopuflags);
hmri_log(sprintf(['INFO: minimum number of echoes required for R2* map calculation is %d.' ...
'\nNumber of PDw echoes available: %d'], mpm_params.neco4R2sfit, length(mpm_params.input(PDwidx).TE)), mpm_params.nopuflags);
if mpm_params.basicR2s
if length(mpm_params.input(PDwidx).TE)<6
hmri_log(sprintf(['GENERAL WARNING: deriving R2* map from an echo train including ' ...
'\nfewer than 6 echoes has not been validated nor investigated. ' ...
'\nFor a robust estimation, the minimum number of echoes required ' ...
'\ndepends on many factors, amongst which: ' ...
'\n\t- SNR/resolution' ...
'\n\t- distribution/spacing between TEs: note that early echoes are more' ...
'\n\t\taffected by the PDw contrast.' ...
'\nInterpret results carefully, with in mind a possible lack of robustness ' ...
'\nand reliability in the R2* estimation.']),mpm_params.defflags);
end
spm_progress_bar('Init',dm(3),'R2* fit','planes completed');
% create nifti object for output R2* map
fR2s = fullfile(calcpath,[outbasename '_' mpm_params.output(mpm_params.qR2s).suffix '.nii']);
Ni = hmri_create_nifti(fR2s, V_pdw(1), dt, mpm_params.output(mpm_params.qR2s).descrip{1});
% Fit R2* decay
for p = 1:dm(3)
data = zeros([dm(1:2),numel(V_pdw)]);
for cecho = 1:numel(V_pdw)
% Allows for reslicing across TE
data(:,:,cecho) = hmri_read_vols(V_pdw(cecho),V_pdw(1),p,mpm_params.interp);
end
% Don't want log of < 1 => assuming typical dynamic range of Dicom
% data, e.g. 12bit @TODO
data = max(data, 1);
R2s = hmri_calc_R2s(struct('data',data,'TE',TE_pdw),mpm_params.R2s_fit_method);
Ni.dat(:,:,p) = max(min(R2s,threshall.R2s),-threshall.R2s)*1000; % threshold T2* at +/- 0.1ms or R2* at +/- 10000 *(1/sec), negative values are allowed to preserve Gaussian distribution
spm_progress_bar('Set',p);
end
spm_progress_bar('Clear');
% Set and write metadata
input_files = mpm_params.input(PDwidx).fnam;
Output_hdr = init_mpm_output_metadata(input_files, mpm_params);
Output_hdr.history.output.imtype = 'Basic R2* map';
Output_hdr.history.output.units = 's-1';
set_metadata(fR2s,Output_hdr,mpm_params.json);
else
hmri_log(sprintf('INFO: No (basic) R2* map will be calculated.\nInsufficient number of PDw echoes.'), mpm_params.defflags);
end
%% =======================================================================%
% Reading and averaging the images
%=========================================================================%
hmri_log(sprintf('\t-------- Reading and averaging the images --------'),mpm_params.nopuflags);
% Average only first few echoes for increased SNR and fit T2*
Pavg = cell(1,mpm_params.ncon);
for ccon=1:mpm_params.ncon % loop over available contrasts
Pavg{ccon} = fullfile(calcpath,[outbasename '_' mpm_params.input(ccon).tag 'w.nii']);
V = spm_vol(mpm_params.input(ccon).fnam);
Ni = hmri_create_nifti(Pavg{ccon}, V(1), dt, ...
sprintf('Averaged %sw images - %d echoes', mpm_params.input(ccon).tag, avg_nr));
spm_progress_bar('Init',dm(3),Ni.descrip,'planes completed');
% sm = 0;
for p = 1:dm(3)
M = spm_matrix([0 0 p]);
Y = zeros(dm(1:2));
for nr = 1:avg_nr
M1 = V(nr).mat\V(1).mat*M;
Y = Y + spm_slice_vol(V(nr),M1,dm(1:2),mpm_params.interp);
end
Ni.dat(:,:,p) = Y/avg_nr;
% sm = sm + sum(Y(:))/avg_nr;
spm_progress_bar('Set',p);
end
% avg = sm/prod(dm);
spm_progress_bar('Clear');
input_files = mpm_params.input(ccon).fnam;
Output_hdr = init_mpm_output_metadata(input_files, mpm_params);
Output_hdr.history.output.imtype = Ni.descrip;
Output_hdr.history.output.units = 'a.u.';
set_metadata(Pavg{ccon},Output_hdr,mpm_params.json);
end
% Average T1w image for PD calculation
% (average over PDproc.nr_echoes_forA echoes, see hmri_defaults):
% @TODO - ensure reference volume is always consistent, e.g. V v's V_pdw to
% ensure we are always in the same reference space of the PDw volume.
if (PDwidx && T1widx)
PT1w_forA = fullfile(calcpath,[outbasename '_T1w_forA.nii']);
V = spm_vol(mpm_params.input(T1widx).fnam);
Ni = hmri_create_nifti(PT1w_forA, V(1), dt, ...
sprintf('Averaged T1w images for PD calculation - %d echoes',PDproc.nr_echoes_forA));
spm_progress_bar('Init',dm(3),Ni.descrip,'planes completed');
for p = 1:dm(3)
Y = zeros(dm(1:2));
for nr = 1:PDproc.nr_echoes_forA
Y = Y + hmri_read_vols(V(nr),V(1),p,mpm_params.interp);
end
Ni.dat(:,:,p) = Y./PDproc.nr_echoes_forA;
spm_progress_bar('Set',p);
end
spm_progress_bar('Clear');
end
%% =======================================================================%
% Coregistering the images
%=========================================================================%
hmri_log(sprintf('\t-------- Coregistering the images --------'),mpm_params.nopuflags);
% NOTE: coregistration can be disabled using the hmri_def.coreg2PDw flag
contrastCoregParams = zeros(mpm_params.ncon, 6);
if mpm_params.coreg
if MTwidx
contrastCoregParams(MTwidx,:) = hmri_coreg(Pavg{PDwidx}, Pavg{MTwidx}, mpm_params.coreg_flags);
end
if T1widx
contrastCoregParams(T1widx,:) = hmri_coreg(Pavg{PDwidx}, Pavg{T1widx}, mpm_params.coreg_flags);
hmri_coreg(Pavg{PDwidx}, PT1w_forA, mpm_params.coreg_flags);
end
end
V_trans = [];
if ~isempty(P_trans)
% Load B1 mapping data if available and coregister to PDw
% P_trans(1,:) = magnitude image (anatomical reference for coregistration)
% P_trans(2,:) = B1 map (p.u.)
if mpm_params.coreg
hmri_coreg(Pavg{PDwidx}, P_trans, mpm_params.coreg_bias_flags);
end
V_trans = spm_vol(P_trans);
end
% parameters saved for quality assessment
if mpm_params.QA.enable
if exist(mpm_params.QA.fnam,'file')
mpm_params.QA = spm_jsonread(mpm_params.QA.fnam);
end
if MTwidx
mpm_params.QA.ContrastCoreg.MT2PD = contrastCoregParams(MTwidx,:);
end
if T1widx
mpm_params.QA.ContrastCoreg.T12PD = contrastCoregParams(T1widx,:);
end
spm_jsonwrite(mpm_params.QA.fnam, mpm_params.QA, struct('indent','\t'));
end
% load averaged images
Vavg(PDwidx) = spm_vol(Pavg{PDwidx});
if MTwidx; Vavg(MTwidx) = spm_vol(Pavg{MTwidx}); end
if T1widx
Vavg(T1widx) = spm_vol(Pavg{T1widx});
VT1w_forA = spm_vol(PT1w_forA);
end
%% =======================================================================%
% multi-contrast R2* map calculation for quality assessment by AL
% METHOD: R2s map calculated for each contrast separately (PDw, MTw, T1w)
% to evaluate the SD of each R2* map in white matter (i.e. intra-run motion
% measurement).
%=========================================================================%
if mpm_params.QA.enable
hmri_log(sprintf('\t-------- Multi-contrast R2* map calculation for QA --------'),mpm_params.nopuflags);
fR2sQA = cell(1,mpm_params.ncon);
for ccon = 1:mpm_params.ncon
% only if enough echoes
if mpm_params.estaticsR2s(ccon)
fR2sQA{ccon} = fullfile(calcpath,[outbasename '_R2s_' mpm_params.input(ccon).tag 'w.nii']);
Ni = hmri_create_nifti(fR2sQA{ccon}, V_pdw(1), dt, ...
'OLS R2* map [s-1]');
Vcon = spm_vol(mpm_params.input(ccon).fnam);
% The assumption is that the result of co-registering the average
% weighted volumes is applicable for each of the echoes of that
% contrast => Replicate the mat field across contrasts for all echoes.
% % matField = cat(3, repmat(VPDw.mat, [1, 1, nPD]), ...
% % repmat(VMTw.mat, [1, 1, nMT]), repmat(VT1w.mat, [1, 1, nT1]));
spm_progress_bar('Init',dm(3),'multi-contrast R2* fit','planes completed');
for p = 1:dm(3)
data = zeros([dm(1:2),numel(Vcon)]);
for cecho = 1:numel(mpm_params.input(ccon).TE)
% Allows for reslicing across TE
data(:,:,cecho) = hmri_read_vols(Vcon(cecho),V_pdw(1),p,mpm_params.interp);
end
% Don't want log of < 1 => assuming typical dynamic range of Dicom
% data, e.g. 12bit @TODO
data = max(data, 1);
R2s = hmri_calc_R2s(struct('data',data,'TE',mpm_params.input(ccon).TE),'ols');
Ni.dat(:,:,p) = max(min(R2s,threshall.R2s),-threshall.R2s)*1000; % threshold T2* at +/- 0.1ms or R2* at +/- 10000 *(1/sec), negative values are allowed to preserve Gaussian distribution
spm_progress_bar('Set',p);
end
spm_progress_bar('Clear');
end
end
end
%% =======================================================================%
% OLS R2* map calculation by MFC
% [Reference: ESTATICS, Weiskopf et al. 2014]
%=========================================================================%
fR2s_OLS = '';
hmri_log(sprintf('\t-------- %s R2* map calculation (ESTATICS) --------',mpm_params.R2s_fit_method),mpm_params.nopuflags);
LogMsg = sprintf(['INFO: minimum number of echoes required for R2* map calculation is %d.' ...
'\nNumber of echoes available: '], mpm_params.neco4R2sfit);
for ccon = 1:mpm_params.ncon
LogMsg = sprintf('%s %sw = %d.',LogMsg, mpm_params.input(ccon).tag, length(mpm_params.input(ccon).TE));
end
hmri_log(LogMsg, mpm_params.nopuflags);
if mpm_params.proc.R2sOLS && any(mpm_params.estaticsR2s)
if length(mpm_params.input(PDwidx).TE)<6
hmri_log(sprintf(['GENERAL WARNING: deriving R2* map from echo trains including ' ...
'\nfewer than 6 echoes has not been validated nor investigated. ' ...
'\nFor a robust estimation, the minimum number of echoes required ' ...
'\ndepends on many factors, amongst which: ' ...
'\n\t- SNR/resolution' ...
'\n\t- distribution/spacing between TEs: note that early echoes are more' ...
'\n\t\taffected by the specific contrast, violating the ESTATICS assumption ' ...
'\n\t\tof a common decay between contrasts.' ...
'\n\t- number of contrasts available (fewer echoes per contrast required ' ...
'\n\t\tfor 3 (PDw, T1w, MTw) contrasts as compared to 2 or even 1) ' ...
'\nInterpret results carefully, with in mind a possible lack of robustness ' ...
'\nand reliability in the R2* estimation.']),mpm_params.defflags);
end
% OLS fit at TE=0: to be used instead of averaged echoes of each
% contrast if "fullOLS" option is enabled
if mpm_params.fullOLS
hmri_log(sprintf('\t-------- and fit to TE=0 for each contrast --------'),mpm_params.nopuflags);
Nmap = nifti;
Pte0 = cell(1,mpm_params.ncon);
for ccon = 1:mpm_params.ncon
% requires a minimum of neco4R2sfit echoes for a robust fit
if mpm_params.estaticsR2s(ccon)
Pte0{ccon} = fullfile(calcpath,[outbasename '_' mpm_params.input(ccon).tag 'w_' mpm_params.R2s_fit_method 'fit_TEzero.nii']);
Ni = hmri_create_nifti(Pte0{ccon}, V_pdw(1), dt, ...
sprintf('%s fit to TE=0 for %sw images - %d echoes', mpm_params.R2s_fit_method, mpm_params.input(ccon).tag, length(mpm_params.input(ccon).TE)));
% set metadata
input_files = mpm_params.input(ccon).fnam;
Output_hdr = init_mpm_output_metadata(input_files, mpm_params);
Output_hdr.history.output.imtype = Ni.descrip;
Output_hdr.history.output.units = 'a.u.';
set_metadata(Pte0{ccon},Output_hdr,mpm_params.json);
% re-load the updated NIFTI file (in case extended header has
% been added, the offset has changed and must be updated before
% writing the data to the file!)
Nmap(ccon) = nifti(Pte0{ccon});
else
Pte0{ccon} = '';
end
end
end % init nifti objects for fullOLS case
fR2s_OLS = fullfile(calcpath,[outbasename '_' mpm_params.output(mpm_params.qR2s).suffix '_' mpm_params.R2s_fit_method '.nii']);
Ni = hmri_create_nifti(fR2s_OLS, V_pdw(1), dt, ...
[mpm_params.R2s_fit_method ' R2* map [s-1]']);
% Combine the data and echo times:
nechoes = zeros(1,mpm_params.ncon);
for ccon = 1:mpm_params.ncon
% only if enough echoes
if mpm_params.estaticsR2s(ccon)
nechoes(ccon) = size(mpm_params.input(ccon).fnam,1);
end
end
spm_progress_bar('Init',dm(3),[mpm_params.R2s_fit_method ' R2* fit','planes completed']);
for p = 1:dm(3)
for ccon = 1:mpm_params.ncon
% only if enough echoes
if mpm_params.estaticsR2s(ccon)
data = zeros([dm(1:2), nechoes(ccon)]);
Vcon = spm_vol(mpm_params.input(ccon).fnam);
for cecho = 1:nechoes(ccon)
% Additionally pass the result of co-registration since
% want output (e.g. intercept) in Vavg space.
% Don't want log of < 1 => assuming typical dynamic range of Dicom
% data, e.g. 12bit @TODO
data(:,:,cecho) = max(hmri_read_vols(Vcon(cecho),V_pdw(1),p,mpm_params.interp, contrastCoregParams(ccon,:)), 1);
end
dataToFit(ccon).data = data;
dataToFit(ccon).TE = mpm_params.input(ccon).TE;
end
end
[R2s, intercepts] = hmri_calc_R2s(dataToFit,mpm_params.R2s_fit_method);
% Writes "fullOLS" images (OLS fit to TE=0 for each contrast)
if mpm_params.fullOLS
for ccon = 1:mpm_params.ncon
% only if enough echoes
if mpm_params.estaticsR2s(ccon)
Nmap(ccon).dat(:,:,p) = intercepts{ccon};
end
end
end
% NB: mat field defined by V_pdw => first PDw echo
% threshold T2* at +/- 0.1ms or R2* at +/- 10000 *(1/sec),
% negative values are allowed to preserve Gaussian distribution.
Ni.dat(:,:,p) = max(min(R2s,threshall.R2s),-threshall.R2s)*1000;
spm_progress_bar('Set',p);
end
spm_progress_bar('Clear');
% Set metadata (R2S_OLS)
input_files = mpm_params.input(PDwidx).fnam;
if (T1widx); input_files = char(input_files, mpm_params.input(T1widx).fnam); end
if (MTwidx); input_files = char(input_files, mpm_params.input(MTwidx).fnam); end
Output_hdr = init_mpm_output_metadata(input_files, mpm_params);
Output_hdr.history.output.imtype = mpm_params.output(mpm_params.qR2s).descrip;
Output_hdr.history.output.units = mpm_params.output(mpm_params.qR2s).units;
set_metadata(fR2s_OLS,Output_hdr,mpm_params.json);
else
hmri_log(sprintf('INFO: No R2* map will be calculated using the ESTATICS model. \nInsufficient number of echoes.'), mpm_params.defflags);
end % OLS code
% if "fullOLS" option enabled AND could be applied to every contrast
% (enough echoes for every contrast), Pte0 images replace Pavg images in
% the rest of the code. We need to apply the substitution and reload the
% images...
if mpm_params.fullOLS && all(mpm_params.estaticsR2s)
for ccon = 1:mpm_params.ncon
Pavg{ccon} = Pte0{ccon};
Vavg(ccon) = spm_vol(Pavg{ccon});
end
end
%% =======================================================================%
% Prepare output for R1, PD and MT maps
%=========================================================================%
Nmap = nifti;
for ii=1:length(mpm_params.output)-1*(~isempty(fR2s)||~isempty(fR2s_OLS)) % R2s output already done
% define NIFTI objects for output images
outputFn = fullfile(calcpath,[outbasename '_' mpm_params.output(ii).suffix '.nii']);
Nmap(ii) = hmri_create_nifti(outputFn, V_pdw(1), dt, ...
fullfile(calcpath,[outbasename '_' mpm_params.output(ii).suffix '.nii']));
end
fR1 = '';
fA = '';
fMT = '';
if (PDwidx && T1widx)
fR1 = fullfile(calcpath,[outbasename '_' mpm_params.output(mpm_params.qR1).suffix '.nii']);
fA = fullfile(calcpath,[outbasename '_' mpm_params.output(mpm_params.qPD).suffix '.nii']);
if MTwidx
fMT = fullfile(calcpath,[outbasename '_' mpm_params.output(mpm_params.qMT).suffix '.nii']);
end
end
%% =======================================================================%
% Map calculation continued (R1, PD, MT)
%=========================================================================%
hmri_log(sprintf('\t-------- Map calculation continued (R1, MTR) --------'), mpm_params.nopuflags);
M0 = Ni.mat;
fa_pdw_rad = fa_pdw * pi / 180;
if MTwidx; fa_mtw_rad = fa_mtw * pi / 180; end
if T1widx; fa_t1w_rad = fa_t1w * pi / 180; end
spm_progress_bar('Init',dm(3),'Calculating maps','planes completed');
% First calculate R1 & MTR
for p = 1:dm(3)
M = M0*spm_matrix([0 0 p]);
% PDw images are always available, so this bit is always loaded:
PDw = spm_slice_vol(Vavg(PDwidx),Vavg(PDwidx).mat\M,dm(1:2),mpm_params.interp);
if ~isempty(V_trans)
f_T = spm_slice_vol(V_trans(2,:),V_trans(2,:).mat\M,dm(1:2),mpm_params.interp)/100; % divide by 100, since p.u. maps
else
f_T = [];
end
% Standard magnetization transfer ratio (MTR) in percent units [p.u.]
% only if trpd = trmt and fapd = famt and if PDw and MTw images are
% available
if (MTwidx && PDwidx)
MTw = spm_slice_vol(Vavg(MTwidx),Vavg(MTwidx).mat\M,dm(1:2),mpm_params.interp);
if (TR_mtw == TR_pdw) && (fa_mtw == fa_pdw) % additional MTR image...
MTR = (PDw-MTw)./(PDw+eps) * 100;
% write MTR image
Nmap(mpm_params.qMTR).dat(:,:,p) = max(min(MTR,threshall.MTR),-threshall.MTR);
end
end
% T1 map and A/PD maps can only be calculated if T1w images are
% available:
if T1widx
T1w = spm_slice_vol(Vavg(T1widx),Vavg(T1widx).mat\M,dm(1:2),mpm_params.interp);
% Transmit bias corrected quantitative T1 values; if f_T empty then semi-quantitative
R1=hmri_calc_R1(struct('data',PDw,'fa',fa_pdw_rad,'TR',TR_pdw,'B1',f_T),...
struct('data',T1w,'fa',fa_t1w_rad,'TR',TR_t1w,'B1',f_T),...
mpm_params.small_angle_approx);
if ISC.enabled&&~isempty(f_T)
% MFC: We do have P2_a and P2_b parameters for this sequence
% => T1 = A(B1) + B(B1)*T1app (see Preibisch 2009)
R1 = R1./(...
+(ISC.P2_a(1)*f_T.^2 + ISC.P2_a(2)*f_T + ISC.P2_a(3)).*R1 ...
+(ISC.P2_b(1)*f_T.^2 + ISC.P2_b(2)*f_T + ISC.P2_b(3)) ...
);
end
R1 = R1*1e6;
R1(R1<0) = 0;
tmp = R1;
tmp(isnan(tmp)) = 0;
Nmap(mpm_params.qR1).dat(:,:,p) = min(max(tmp,-threshall.R1),threshall.R1)*1e-3; % truncating images
end
spm_progress_bar('Set',p);
end
spm_progress_bar('Clear');
% apply UNICORT if required:
% NOTE: the (unicort) B1 map will then be used to calculate PD, but
% not to correct R1 map from RF spoiling effects (unsuitable
% interfering and higher order corrections!) nor to correct MT map
% from higher order transmit bias effects.
V_trans_unicort = [];
if (mpm_params.UNICORT.R1)
out_unicort = hmri_create_unicort(Pavg{PDwidx}, fR1, jobsubj);
fR1 = out_unicort.R1u{1};
P_trans_unicort = out_unicort.B1u{1};
V_trans_unicort = spm_vol(P_trans_unicort);
Output_hdr = get_metadata(fR1);
Output_hdr{1}.history.output.imtype = mpm_params.output(mpm_params.qR1).descrip;
set_metadata(fR1,Output_hdr{1},mpm_params.json);
end
V_R1 = spm_vol(fR1);
hmri_log(sprintf('\t-------- Map calculation continued (MT, A/PD) --------'), mpm_params.nopuflags);
spm_progress_bar('Init',dm(3),'Calculating maps (continued)','planes completed');
% Then calculate MT & PD
for p = 1:dm(3)
M = M0*spm_matrix([0 0 p]);
if ~isempty(V_trans)
f_T = spm_slice_vol(V_trans(2,:),V_trans(2,:).mat\M,dm(1:2),mpm_params.interp)/100; % divide by 100, since p.u. maps
elseif ~isempty(V_trans_unicort)
f_T = spm_slice_vol(V_trans_unicort(1,:),V_trans_unicort(1,:).mat\M,dm(1:2),mpm_params.interp)/100; % divide by 100, since p.u. maps
else
f_T = [];
end
% PDw images are always available, so this bit is always loaded:
PDw = spm_slice_vol(Vavg(PDwidx),Vavg(PDwidx).mat\M,dm(1:2),mpm_params.interp);
% T1 map and A/PD maps can only be calculated if T1w images are
% available:
if T1widx
T1w = spm_slice_vol(Vavg(T1widx),Vavg(T1widx).mat\M,dm(1:2),mpm_params.interp);
% if "fullOLS" option enabled, use the OLS fit at TE=0 as
% "T1w_forA"; otherwise use the average calculated earlier (by
% default, corresponds to the first echo to reduce R2* bias)
if mpm_params.fullOLS
T1w_forA = T1w;
else
T1w_forA = spm_slice_vol(VT1w_forA,VT1w_forA.mat\M,dm(1:2),mpm_params.interp);
end
% A values proportional to PD
% f_T correction is applied either if:
% - f_T has been provided as separate B1 mapping measurement (not
% UNICORT!) or
% - f_T has been calculated using UNICORT *AND* the UNICORT.PD flag
% is enabled (advanced user only! method not validated yet!)
if(~isempty(f_T))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.PD)
A = hmri_calc_A(struct('data',PDw,'fa',fa_pdw_rad,'TR',TR_pdw,'B1',f_T),...
struct('data',T1w_forA,'fa',fa_t1w_rad,'TR',TR_t1w,'B1',f_T),...
mpm_params.small_angle_approx);
else
% semi-quantitative A
A = hmri_calc_A(struct('data',PDw,'fa',fa_pdw_rad,'TR',TR_pdw,'B1',1),...
struct('data',T1w_forA,'fa',fa_t1w_rad,'TR',TR_t1w,'B1',1),...
mpm_params.small_angle_approx);
end
tmp = A;
tmp(isinf(tmp)) = 0;
tmp(isnan(tmp)) = 0;
Nmap(mpm_params.qPD).dat(:,:,p) = max(min(tmp,threshall.A),-threshall.A);
% dynamic range increased to 10^5 to accommodate phased-array coils and symmetrical for noise distribution
% for MT maps calculation, one needs MTw images on top of the T1w
% and PDw ones...
if MTwidx
MTw = spm_slice_vol(Vavg(MTwidx),Vavg(MTwidx).mat\M,dm(1:2),3);
R1_forMT = hmri_calc_R1(struct('data',PDw,'fa',fa_pdw_rad,'TR',TR_pdw,'B1',1),...
struct('data',T1w,'fa',fa_t1w_rad,'TR',TR_t1w,'B1',1),...
mpm_params.small_angle_approx);
R1_forMT(isnan(R1_forMT))=0;
A_forMT = hmri_calc_A(struct('data',PDw,'fa',fa_pdw_rad,'TR',TR_pdw,'B1',1),...
struct('data',T1w_forA,'fa',fa_t1w_rad,'TR',TR_t1w,'B1',1),...
mpm_params.small_angle_approx);
A_forMT(isnan(A_forMT))=0;
% MT in [p.u.]; offset by - famt * famt / 2 * 100 where MT_w = 0 (outside mask)
MT = ( (A_forMT * fa_mtw_rad - MTw) ./ (MTw+eps) .* R1_forMT*TR_mtw - fa_mtw_rad^2 / 2 ) * 100;
% f_T correction is applied either if:
% - f_T has been provided as separate B1 mapping measurement (not
% UNICORT!) or
% - f_T has been calculated using UNICORT *AND* the UNICORT.MT flag
% is enabled (advanced user only! method not validated yet!)
if (~isempty(f_T))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.MT)
MT = MT .* (1 - 0.4) ./ (1 - 0.4 * f_T);
end
tmp = MT;
Nmap(mpm_params.qMT).dat(:,:,p) = max(min(tmp,threshall.MT),-threshall.MT);
end
end
spm_progress_bar('Set',p);
end
spm_progress_bar('Clear');
% set metadata for all output images
input_files = mpm_params.input(PDwidx).fnam;
if (T1widx); input_files = char(input_files, mpm_params.input(T1widx).fnam); end
if (MTwidx); input_files = char(input_files, mpm_params.input(MTwidx).fnam); end
Output_hdr = init_mpm_output_metadata(input_files, mpm_params);
for ctr = 1:length(mpm_params.output)-1
Output_hdr.history.output.imtype = mpm_params.output(ctr).descrip;
Output_hdr.history.output.units = mpm_params.output(ctr).units;
set_metadata(fullfile(calcpath,[outbasename '_' mpm_params.output(ctr).suffix '.nii']),Output_hdr,mpm_params.json);
end
%% =======================================================================%
% ACPC Realign all images - only if MT map created
%=========================================================================%
if mpm_params.ACPCrealign
if (MTwidx && PDwidx && T1widx)
hmri_log(sprintf('\t-------- ACPC Realign all images --------'), mpm_params.nopuflags);
% Define and calculate masked MT image
% Load MT image
V_MT = spm_vol(fMT);
MTimage = spm_read_vols(V_MT);
% Define new file name for masked MT image
V_MT.fname = fullfile(calcpath,['masked_' spm_str_manip(fMT,'t')]);
% Load average PDw image (mask based on averaged PDw image)
PDWimage = spm_read_vols(Vavg(PDwidx));
% Mask MT image and save the masked MT image ('masked_..._MT.nii')
MTimage(PDWimage<0.6*mean(PDWimage(:)))=0;
spm_write_vol(V_MT,MTimage);
% Use masked MT image to calculate transformation for ACPC realignment
% (to increase robustness in segmentation):
[~,R] = hmri_create_comm_adjust(1,V_MT.fname,V_MT.fname,8,0,fullfile(spm('Dir'),'canonical','avg152T1.nii'));
% Collect all images from all defined output directories:
allpaths = jobsubj.path;
ACPC_images = [];
f = fieldnames(allpaths);
for cfield = 1:length(f)
tmpfiles = spm_select('FPList',allpaths.(f{cfield}),'.*.(img|nii)$');
if ~isempty(tmpfiles)
if ~isempty(ACPC_images)
ACPC_images = char(ACPC_images,spm_select('FPList',allpaths.(f{cfield}),'.*.(img|nii)$'));
else
ACPC_images = tmpfiles;
end
end
end
% Apply transform to ALL images
for i=1:size(ACPC_images,1)
spm_get_space(deblank(ACPC_images(i,:)),...
R*spm_get_space(deblank(ACPC_images(i,:))));
Vsave = spm_vol(ACPC_images(i,:));
Vsave.descrip = [Vsave.descrip ' - AC-PC realigned'];
spm_write_vol(Vsave,spm_read_vols(spm_vol(ACPC_images(i,:))));
end
% Save transformation matrix
spm_jsonwrite(fullfile(supplpath,'hMRI_map_creation_ACPCrealign_transformation_matrix.json'),R,struct('indent','\t'));
else
hmri_log(sprintf(['WARNING: ACPC Realign was enabled, but no MT map was available \n' ...
'to proceed. ACPC realignment must be done separately, e.g. you can \n'...
'run [hMRI tools > Auto-reorient] before calculating the maps.\n' ...
'NOTE: segmentation might crash if no initial reorientation.']), mpm_params.defflags);
end
end
% for quality assessment and/or PD map calculation
% segmentation preferentially performed on MT map but can be done on R1 map
% if no MT map available. Therefore, we must at least have R1 available,
% i.e. both PDw and T1w inputs...
if (mpm_params.QA.enable||(PDproc.calibr)) && (PDwidx && T1widx)
if ~isempty(fMT)
Vsave = spm_vol(fMT);
threshMT=threshall.MT;
else % ~isempty(fR1);
Vsave = spm_vol(fR1);
threshMT=threshall.R1*1e3;
end
MTtemp = spm_read_vols(Vsave);
% The 5 outer voxels in all directions are nulled in order to remove
% artefactual effects from the MT map on segmentation:
MTtemp(1:5,:,:)=0; MTtemp(end-5:end,:,:)=0;
MTtemp(:,1:5,:)=0; MTtemp(:,end-5:end,:)=0;
MTtemp(:,:,1:5)=0; MTtemp(:,:,end-5:end)=0;
% Null very bright and negative voxels
MTtemp(abs(MTtemp)==threshMT)=0;
MTtemp(isinf(MTtemp))=0;
MTtemp(isnan(MTtemp))=0;
MTtemp(MTtemp<0)=0;
Vsave.fname = spm_file(Vsave.fname,'suffix','_outer_suppressed');
spm_write_vol(Vsave,MTtemp);
% use unified segmentation with uniform defaults across the toobox:
job_brainmask = hmri_get_defaults('segment');
job_brainmask.channel.vols = {Vsave.fname};
job_brainmask.channel.write = [0 0]; % no need to write BiasField nor BiasCorrected image
output_list = spm_preproc_run(job_brainmask);
fTPM = char(cat(1,output_list.tiss.c));
end
% for quality assessment - the above segmentation must have run
if mpm_params.QA.enable && exist('fTPM','var') && any(mpm_params.estaticsR2s)
% Load existing QA results
if exist(mpm_params.QA.fnam,'file')
mpm_params.QA = spm_jsonread(mpm_params.QA.fnam);
end
% Calculate WM mask
TPMs = spm_read_vols(spm_vol(fTPM));
WMmask = zeros(size(squeeze(TPMs(:,:,:,2))));
WMmask(squeeze(TPMs(:,:,:,2))>=PDproc.WMMaskTh) = 1;
WMmask = spm_erode(spm_erode(double(WMmask)));
% Load OLS R2s maps calculated for each contrast, mask them and
% calculate SD within the WM mask (measure of the intra-run motion for
% each contrast)
for ccon = 1:mpm_params.ncon
R2s = spm_read_vols(spm_vol(fR2sQA{ccon}));
MaskedR2s = squeeze(R2s.*WMmask);
SDR2s = std(MaskedR2s(MaskedR2s~=0),[],1);
mpm_params.QA.SDR2s.([mpm_params.input(ccon).tag 'w']) = SDR2s;
end
spm_jsonwrite(mpm_params.QA.fnam, mpm_params.QA, struct('indent','\t'));
end
% PD map calculation
% for correction of the R2s bias in the A map if that option is enabled...
if PDproc.T2scorr && (~isempty(fR2s)||~isempty(fR2s_OLS))
% uses OLS it if available - less noisy
if ~isempty(fR2s_OLS)
PR2s = fR2s_OLS;
else
PR2s = fR2s;
end
% calculate correction (expected to be between 1 and 1.5 approx)
R2s = spm_read_vols(spm_vol(PR2s));
R2scorr4A = zeros(size(R2s));
for cecho=1:mpm_params.proc.PD.nr_echoes_forA
TE = mpm_params.input(PDwidx).TE(cecho)*0.001; % in seconds
R2scorr4A = R2scorr4A + exp(-TE.*R2s);
end
R2scorr4A = R2scorr4A/mpm_params.proc.PD.nr_echoes_forA;
% save correction for inspection
fR2scorr4A = spm_file(PR2s,'suffix','_corr4A');
NiR2scorr4A = hmri_create_nifti(fR2scorr4A, V_pdw(1), dt, ...
'R2* bias correction factor for A map (T2scorr option)');
NiR2scorr4A.dat(:,:,:) = R2scorr4A;
% apply correction
fAcorr = spm_file(fA,'suffix','_R2scorr');
NiAcorr = hmri_create_nifti(fAcorr, V_pdw(1), dt, ...
'R2* bias corrected A map (T2scorr option)');
tmp = spm_read_vols(spm_vol(fA))./(R2scorr4A+eps);
tmp(isnan(tmp)|isinf(tmp)) = 0;
tmp = max(min(tmp,threshall.A),-threshall.A);
NiAcorr.dat(:,:,:) = tmp;
fA = fAcorr;
end
% PD map calculation continued
if ~isempty(f_T) && ~isempty(fA) && exist('fTPM','var') && (mpm_params.UNICORT.PD || ~mpm_params.UNICORT.R1)
% if calibration enabled, do the Unified Segmentation bias correction
% if required and calibrate the PD map
if PDproc.calibr
PDcalculation(fA,fTPM,mpm_params);
end
end
% copy final result files into Results directory
% NB: to avoid ambiguity for users, only the 4 final maps to be used in
% further analysis are in Results, all other files (additional maps,
% processing parameters, etc) are in Results/Supplementary.
if ~isempty(fR1)
fR1_final = fullfile(respath, spm_file(fR1,'filename'));
try copyfile(fR1,fR1_final); end
try copyfile([spm_str_manip(fR1,'r') '.json'],[spm_str_manip(fR1_final,'r') '.json']); end %#ok<*TRYNC>
fR1 = fR1_final;
end
if ~isempty(fR2s)
fR2s_final = fullfile(respath, spm_file(fR2s,'filename'));
copyfile(fR2s,fR2s_final);
try copyfile([spm_str_manip(fR2s,'r') '.json'],[spm_str_manip(fR2s_final,'r') '.json']); end
fR2s = fR2s_final;
end
% NB: if OLS calculation of R2s map has been done, the output file for R2s
% map is the OLS result. In that case, the basic R2s map is moved to
% Results/Supplementary while the R2s_OLS is copied into results directory:
if mpm_params.proc.R2sOLS && ~isempty(fR2s_OLS)
% move basic R2s map to Results/Supplementary
movefile(fR2s_final, fullfile(supplpath, spm_file(fR2s_final,'filename')));
try movefile([spm_str_manip(fR2s_final,'r') '.json'],fullfile(supplpath, [spm_file(fR2s_final,'basename') '.json'])); end
% copy OLS_R2s map to Results
fR2s_OLS_final = fullfile(respath, spm_file(fR2s_OLS,'filename'));
copyfile(fR2s_OLS,fR2s_OLS_final);
try copyfile([spm_str_manip(fR2s_OLS,'r') '.json'],[spm_str_manip(fR2s_OLS_final,'r') '.json']); end
% the hmri_create_MTProt fR2s output in now the R2s_OLS map
fR2s = fR2s_OLS_final;
end
if ~isempty(fMT)
fMT_final = fullfile(respath, spm_file(fMT,'filename'));
copyfile(fMT,fMT_final);
try copyfile([spm_str_manip(fMT,'r') '.json'],[spm_str_manip(fMT_final,'r') '.json']); end
fMT = fMT_final;
end
if ~isempty(fA)
fA_final = fullfile(respath, spm_file(fA,'filename'));
copyfile(fA,fA_final);
try copyfile([spm_str_manip(fA,'r') '.json'],[spm_str_manip(fA_final,'r') '.json']); end
fA = fA_final;
end
PPDw_final = fullfile(supplpath, spm_file(Pavg{PDwidx},'filename'));
copyfile(Pavg{PDwidx},PPDw_final);
try copyfile([spm_str_manip(Pavg{PDwidx},'r') '.json'],[spm_str_manip(PPDw_final,'r') '.json']); end
PPDw = PPDw_final;
if T1widx
PT1w_final = fullfile(supplpath, spm_file(Pavg{T1widx},'filename'));
copyfile(Pavg{T1widx},PT1w_final);
try copyfile([spm_str_manip(Pavg{T1widx},'r') '.json'],[spm_str_manip(PT1w_final,'r') '.json']); end
PT1w = PT1w_final;
else
PT1w = '';
end
if MTwidx
PMTw_final = fullfile(supplpath, spm_file(Pavg{MTwidx},'filename'));
copyfile(Pavg{MTwidx},PMTw_final);
try copyfile([spm_str_manip(Pavg{MTwidx},'r') '.json'],[spm_str_manip(PMTw_final,'r') '.json']); end
PMTw = PMTw_final;
else
PMTw = '';
end
% save processing params (mpm_params)
spm_jsonwrite(fullfile(supplpath,'hMRI_map_creation_mpm_params.json'),mpm_params,struct('indent','\t'));
spm_progress_bar('Clear');
hmri_log(sprintf('\t============ MPM PROCESSING: completed (%s) ============', datestr(now)),mpm_params.nopuflags);
end
%% =======================================================================%
% Proton density map calculation
%=========================================================================%
function PDcalculation(fA, fTPM, mpm_params)
% fA is the filename of the output A image
% (not yet properly quantitative PD map when it enters PDcalculation)
% fTMP is the list of TPMs generated from MT map
hmri_log(sprintf('\t-------- Proton density map calculation --------'), mpm_params.nopuflags);
PDproc = mpm_params.proc.PD;
threshA = mpm_params.proc.threshall.A;
calcpath = mpm_params.calcpath;
TPMs = spm_read_vols(spm_vol(fTPM));
WBmask = zeros(size(squeeze(TPMs(:,:,:,1))));
WBmask(sum(cat(4,TPMs(:,:,:,1:2),TPMs(:,:,:,end)),4)>=PDproc.WBMaskTh) = 1;
WMmask=zeros(size(squeeze(TPMs(:,:,:,1))));
WMmask(squeeze(TPMs(:,:,:,2))>=PDproc.WMMaskTh) = 1;
% Save masked A map for bias-field correction later
V_maskedA = spm_vol(fA);
V_maskedA.fname = fullfile(calcpath,['masked_' spm_str_manip(V_maskedA.fname,'t')]);
maskedA = spm_read_vols(spm_vol(fA)).*WBmask;
maskedA(isinf(maskedA)) = 0;
maskedA(isnan(maskedA)) = 0;
maskedA(maskedA==threshA) = 0;
maskedA(maskedA<0) = 0;
spm_write_vol(V_maskedA,maskedA);
% Bias-field correction of masked A map
% use unified segmentation with uniform defaults across the toolbox:
job_bfcorr = hmri_get_defaults('segment');
job_bfcorr.channel.vols = {V_maskedA.fname};
job_bfcorr.channel.biasreg = PDproc.biasreg;
job_bfcorr.channel.biasfwhm = PDproc.biasfwhm;
job_bfcorr.channel.write = [1 0]; % need the BiasField, obviously!
for ctis=1:length(job_bfcorr.tissue)
job_bfcorr.tissue(ctis).native = [0 0]; % no need to write c* volumes
end
output_list = spm_preproc_run(job_bfcorr);