-
Notifications
You must be signed in to change notification settings - Fork 27
/
nii_makeDTI.m
177 lines (170 loc) · 4.76 KB
/
nii_makeDTI.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
function nii_makeDTI (RASorder)
%Generates a simple NIfTI format DTI image
% RASorder : (optional) if false then columns:rows:slices not in right:anterior:superior order
%Examples
% nii_makeDTI; %order: RAS
% nii_makeDTI(1) %order: ARS
% nii_makeDTI(2); %order: LAS
% nii_makeDTI(3); %order: SRA
fnm = 'test';
pixDimMM = 3; %distance between voxel centers
%bvec = makeBVec7; %6 diffusion bvectors small file but somewhat uneven
bvec = makeBVec33; %32 diffusion bvectors large file
nVol = size(bvec,1);
dim = [9, 9, 9, nVol]; %image resolution in columns, rows, slices, volumes
%create matrix
mat = eye(4);
mat(1:3, 1:3) = mat(1:3, 1:3) * pixDimMM;
mat(1,4) = -pixDimMM * ((dim(1)+1)/2); %center, indexed from 1 not 0
mat(2,4) = -pixDimMM * ((dim(2)+1)/2);
mat(3,4) = -pixDimMM * ((dim(3)+1)/2);
%SPM's matrices are indexed from 1
% hdr = spm_vol('test.nii,1');
% mm = (hdr.mat * [0 0 1 1; 0 0 9 1]') %mm of 1st and 9th slice (z=-12..12)
%if you use fslhd you will see matrix indexed from 0
% >fslhd test.nii
% mat0 = [3 0 0 -12; 0 3 0 -12; 0 0 3 -12; 0 0 0 1]
% mm = (mat0 * [0 0 0 1; 0 0 8 1]') %mm of 1st and 9th slice (z=-12..12)
if exist('RASorder','var') && (RASorder > 0) && (RASorder < 4)
%permute rows so no longer RAS
if RASorder == 1
mat = mat(:,[2 1 3 4]); %order: ARS
fnm = [fnm 'ARS'];
elseif RASorder == 2
mat(1,:) = -mat(1,:);
fnm = [fnm 'LAS'];
elseif RASorder == 3
mat = mat(:,[3 1 2 4]); %order: SRA
fnm = [fnm 'SRA'];
end;
else
fnm = [fnm 'RAS'];
end
%create BVec/BVal files
writeBvecBval (bvec, fnm);
%create NIfTI image
dtype = 32; %precision of data
ofile = [fnm, '.nii'];%spm_file(parfile,'path',opts.outdir,'ext',opts.ext);
scale = 1;
inter = 0;
switch dtype
case 8
dtype = spm_type('int8'); % uint8?
case 16
dtype = spm_type('int16'); % uint16?
case 32
dtype = spm_type('float32');
case 64
dtype = spm_type('float64');
otherwise
error('Unknown data type.');
end
dato = file_array(ofile,dim,[dtype 0],0,scale,inter);
N = nifti;
N.dat = dato;
N.mat = mat;
N.mat0 = mat;
N.mat_intent = 'Scanner';
N.mat0_intent = 'Scanner';
N.descrip = 'made with Matlab';
%N.timing = struct('toffset',[],'tspace',[]); % store TR
create(N);
for i=1:nVol
N.dat(:,:,:,i) = makeDtiSub(dim(1:3), bvec(i,:));
end
%end nii_mk()
function writeBvecBval (bvec, fnm)
dlmwrite([fnm '.bvec'],bvec', 'delimiter','\t');
bval = zeros(size(bvec,1),1);
for i = 1 : size(bvec,1) %B=0 images have zero-length vectors, otherwise B-weighted
if norm(bvec(i,:)) >= 0.01
bval(i) = 1000;
end
end
dlmwrite([fnm '.bval'],bval', 'delimiter','\t');
dlmwrite([fnm '.trakvis.bvec'],bvec(bval > 0,:), 'delimiter','\t');
%end writeBvecBval()
function bvec = makeBVec33
%Returns a set of bvectors: 0,0,0 for B=0
% http://www.emmanuelcaruyer.com/q-space-sampling.php
bvec=[0 0 0;
0.049 -0.919 -0.391
0.726 0.301 -0.618
-0.683 0.255 -0.684
0.845 -0.502 -0.186
-0.73 -0.619 -0.288
-0.051 0.039 0.998
-0.018 0.871 -0.491
-0.444 0.494 0.747
-0.989 -0.086 -0.116
-0.47 -0.855 0.221
0.412 0.4 0.819
-0.552 0.79 -0.267
-0.123 -0.477 0.871
-0.848 0.141 0.51
-0.341 -0.788 -0.512
0.361 -0.529 0.768
-0.472 0.85 0.234
-0.856 -0.481 0.189
0.797 0.162 0.582
0.467 -0.009 -0.884
0.013 0.998 -0.056
0.882 -0.387 0.267
0.017 -0.536 -0.844
-0.442 -0.651 0.617
0.365 -0.058 0.929
0.977 -0.004 -0.213
-0.406 -0.902 -0.145
-0.627 0.614 0.479
-0.354 0.772 -0.528
-0.658 -0.472 -0.586
0.423 0.322 -0.847
0.212 -0.754 -0.622];
%end makeBVec33()
function bvec = makeBVec7
%Returns a set of bvectors: 0,0,0 for B=0
% http://www.emmanuelcaruyer.com/q-space-sampling.php
bvec=[0 0 0;
0.049 -0.919 -0.391;
0.726 0.301 -0.618;
-0.683 0.255 -0.684;
0.845 -0.502 -0.186;
-0.730 -0.619 -0.288;
-0.051 0.039 0.998;
-0.018 0.871 -0.491];
%end makeBVec7()
function img = makeDtiSub (dim, bvec)
isBZero = norm(bvec) < 0.01; %if B=0 then unweighted image
img = zeros(dim);
center = (dim+1)/2; %e.g. center of 9 voxels is the 5th
sz = norm(center);
for z = 1 : dim(3)
vec = z/dim(3);
if vec < 0.5
vec = [0 0 1]; %head/foot
else
vec = [1 0 0]; %left/right
end
dwi = norm(cross(bvec,vec)); %is this orthogonal
if isBZero
dwi = 1;
end
%fprintf('%g\n', dwi);
for y = 1 : dim(2)
yRamp = (y/dim(2)) * 2; %make low rows dimmer
%if yRamp > 1, yRamp = 1; end;
for x = 1 : dim(1)
dx = (sz -norm(center-[x,y,z]))/sz;
if (~isBZero) && (dx < 0.5)
dx = 0;
end
if yRamp < 1.0 %no FA in dark rows
img(x,y,z)= dx * yRamp;
else
img(x,y,z)= dwi * dx;
end;
end;
end;
end;
img = img * 4096; %TrackVis likes large integer values
%end makeDtiSub()