forked from rordenlab/spmScripts
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnii_dtibatch.m
157 lines (149 loc) · 6.54 KB
/
nii_dtibatch.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
function nii_dtibatch (dtiBvecNames, isEddyCorrect)
%quick processing of DTI - uses faster eddy_correct rather than eddy/topup
%assumes angulations have been correctly adjusted
% dtiNii: name of bvec file(s), e.g. img.bvec
% isEddyCorrect : if true simple undistortion applied, if false than quick and dirty
%Examples
% nii_dtibatch
% nii_dtibatch('DTI_axial_04.bvec');
% nii_dtibatch(strvcat('DTI_axial_04.bvec','DTI_30_deg_AP_05.bvec'));
%Chris Rorden 4 July 2014, BSD License
%Subsequently you can view the processed images from the command line
% fslview DTI_axial_04_FA.nii.gz DTI_axial_04_V1.nii.gz
if ~exist('dtiBvecNames','var') %file not specified
[A,Apth] = uigetfile({'*.bvec';'*.*'},'Select b-vector file(s)', 'MultiSelect', 'on');
dtiBvecNames = strcat(Apth,char(A));
end;
if ~exist('isEddyCorrect','var') %file not specified
isEddyCorrect = false;
end;
fsldir= '/usr/local/fsl';
if ~exist(fsldir,'dir')
error('Unable to find %s', fsldir);
end
for i=1:size(dtiBvecNames,1)
dtiBvec = deblank(dtiBvecNames(i,:)); %positive image
[pth,nam,ext] = fileparts(dtiBvec);
imgNam = fullfile(pth, [nam '.nii']); %img.nii
if ~exist(imgNam,'file')
imgNam = fullfile(pth, [nam '.nii.gz']); %img.nii.gz
end;
if ~exist(imgNam,'file'), warning('Unable to find %s\n', imgNam); continue; end;
refVol = refVolSub(dtiBvec);
%next: permute all possible b-vector alternatives
dtiSub(fsldir,imgNam,dtiBvec,refVol, isEddyCorrect);
%tracktionSub(dtiBvec);
end
viewSub(fsldir,dtiBvec) %display results
%end main loop.... subroutines follow
function tracktionSub(dtiBvec)
if strcmpi(computer, 'GLNXA64')
exeNam = 'tracktionLX'
else
exeNam = 'tracktion';
end
exeNam = fullfile(fileparts(which(mfilename)), exeNam);
p = fileparts(exeNam);
if isempty(p)
exeNam = fullfile(pwd, exeNam);
end
if ~exist(exeNam,'file')
fprintf('Skipped tractography: unable to find %s\n', exeNam);
return;
end
command=sprintf('%s "%s"\n',exeNam, dtiBvec);
system(command);
%end tracktionSub()
function ref = refVolSub(vNam) %find first B0 volume
[pth,nam,ext] = fileparts(vNam);
bNam = fullfile(pth, [ nam '.bval'] ); %Eddy corrected data
if ~exist(bNam, 'file'), error('Unable to find file %s',bNam); end;
b = textread(bNam);
ref = find(b==0, 1, 'first');
if isempty(ref)
m = min(b);
if m < 11 %some systems report b=~5 for the "b=0" volumes
ref = find(b==m, 1, 'first');
end
end;
if isempty(ref), error('No b-zero volume in file %s', bNam); end;
ref = ref - 1;%fsl indexes volumes from 0
%end refVolSub()
function maskNam = betSub(fsldir,imgNam, refVol) %brain extract
setenv('FSLDIR', fsldir);
[pth,nam,ext] = fileparts(imgNam);
if strcmpi(ext,'.gz'), [~,nam] = fileparts(nam); end;
inNam = imgNam;
maskNam = fullfile(pth, nam ); %will generate image "dti_mask.nii.gz"
if (exist('refVol','var')) && (refVol > 0)
%fslroi WIPDiffHR2.2SENSE4.2_ID123_12 msk 32 1
refNam = fullfile(pth, [nam '_' num2str(refVol)] ); %will generate image "dti_mask.nii.gz"
command=sprintf('sh -c ". ${FSLDIR}/etc/fslconf/fsl.sh; ${FSLDIR}/bin/fslroi %s %s %d 1"\n',imgNam,refNam, refVol);
system(command);
inNam = refNam;
end
command=sprintf('sh -c ". ${FSLDIR}/etc/fslconf/fsl.sh; ${FSLDIR}/bin/bet %s %s -f 0.3 -g 0 -n -m"\n',inNam,maskNam);
maskNam = fullfile(pth, [nam '_mask.nii.gz']); %will generate image "dti_mask.nii.gz"
system(command);
fprintf(command);
%end betSub()
function eccNam = eddySub(fsldir,imgNam,refVol) %eddy current correct and brain extract
fprintf('Using eddy_correct: consider using eddy (and topup)\n');
[pth,nam,ext] = fileparts(imgNam);
eccNam = fullfile(pth, [nam, '_ecc']); %will generate image "dti_eddy.nii.gz"
setenv('FSLDIR', fsldir);
setenv('PATH', [getenv('PATH') ':/usr/local/fsl/bin'])
command=sprintf('sh -c ". ${FSLDIR}/etc/fslconf/fsl.sh; ${FSLDIR}/bin/eddy_correct %s %s %d"\n',imgNam,eccNam, refVol);
system(command);
eccNam = fullfile(pth, [nam, '_ecc.nii.gz']); %Eddy corrected data
%end eddySub()
function dtiSub(fsldir,imgNam,vNam, refVol, isEddyCorrect) %compute vectors
%%/usr/local/fsl/bin/dtifit --data=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/dti_eddy.nii.gz --out=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/dti --mask=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/dti_mask.nii.gz --bvecs=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/s004a001.bvec --bvals=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/s003a001.bval
[pth,nam,ext] = fileparts(vNam);
bNam = fullfile(pth, [ nam '.bval'] ); %Eddy corrected data
maskNam = betSub(fsldir,imgNam, refVol);
if ~exist(maskNam, 'file'), error('BET failed to create %s', maskNam); end;
if isEddyCorrect
error('x');
eccNam = eddySub(fsldir,imgNam,refVol);
else
eccNam = imgNam;
end
[pth,nam,ext] = fileparts(vNam);
outNam = fullfile(pth, nam); %Eddy corrected data
setenv('FSLDIR', fsldir);
setenv('PATH', [getenv('PATH') ':/usr/local/fsl/bin'])
command=sprintf('sh -c ". ${FSLDIR}/etc/fslconf/fsl.sh; ${FSLDIR}/bin/dtifit --save_tensor --data=%s --out=%s --mask=%s --bvecs=%s --bvals=%s"\n',eccNam, outNam, maskNam,vNam,bNam);
system(command);
fprintf(command);
delete(fullfile(pth, [nam '_V2.nii.gz']));
delete(fullfile(pth, [nam '_V3.nii.gz']));
delete(fullfile(pth, [nam '_L1.nii.gz']));
delete(fullfile(pth, [nam '_L2.nii.gz']));
delete(fullfile(pth, [nam '_L3.nii.gz']));
delete(fullfile(pth, [nam '_MO.nii.gz']));
%delete(fullfile(pth, [nam '_MD.nii.gz']));
delete(fullfile(pth, [nam '_S0.nii.gz']));
%end dtiSub
function viewSub(fsldir,vNam) %compute vectors
%%/usr/local/fsl/bin/dtifit --data=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/dti_eddy.nii.gz --out=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/dti --mask=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/dti_mask.nii.gz --bvecs=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/s004a001.bvec --bvals=/Users/rorden/Desktop/sliceOrder/dicom2/dtitest/s003a001.bval
[pth,nam,ext] = fileparts(vNam);
faNam = fullfile(pth, [nam '_FA.nii.gz']); %Eddy corrected data
v1Nam = fullfile(pth, [nam '_V1.nii.gz']); %Eddy corrected data
setenv('FSLDIR', fsldir);
setenv('PATH', [getenv('PATH') ':/usr/local/fsl/bin'])
exenam = fullfile(fsldir, 'bin', 'fsleyes');
%>fsleyes tosh_DTI_ortho_20191203103319_7006_FA.nii.gz tosh_DTI_ortho_20191203103319_7006_V1.nii.gz
if ~exist(exenam)
exenam = fullfile(fsldir, 'bin', 'fslview_deprecated');
end
if ~exist(exenam)
exenam = fullfile(fsldir, 'bin', 'fslview');
end
if ~exist(exenam)
error('Unable to find fslview - update script to use fsleyes');
end
command=sprintf('sh -c ". ${FSLDIR}/etc/fslconf/fsl.sh; %s %s %s &"\n',exenam, faNam,v1Nam);
%system(command);
fprintf(command);
%end dtiSub