Browse code

searchlight ready. missing nifti-image-write support. added timing to FBS

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@167 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on 30/03/2009 17:54:25
Showing 40 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,554 @@
1
+%  Using 2D or 3D affine matrix to rotate, translate, scale, reflect and
2
+%  shear a 2D image or 3D volume. 2D image is represented by a 2D matrix,
3
+%  3D volume is represented by a 3D matrix, and data type can be real 
4
+%  integer or floating-point.
5
+%
6
+%  You may notice that MATLAB has a function called 'imtransform.m' for
7
+%  2D spatial transformation. However, keep in mind that 'imtransform.m'
8
+%  assumes y for the 1st dimension, and x for the 2nd dimension. They are
9
+%  equivalent otherwise.
10
+%
11
+%  In addition, if you adjust the 'new_elem_size' parameter, this 'affine.m'
12
+%  is equivalent to 'interp2.m' for 2D image, and equivalent to 'interp3.m'
13
+%  for 3D volume.
14
+%
15
+%  Usage: [new_img new_M] = ...
16
+%	affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);
17
+%
18
+%  old_img  -	original 2D image or 3D volume. We assume x for the 1st
19
+%		dimension, y for the 2nd dimension, and z for the 3rd
20
+%		dimension.
21
+%
22
+%  old_M  -	a 3x3 2D affine matrix for 2D image, or a 4x4 3D affine
23
+%		matrix for 3D volume. We assume x for the 1st dimension,
24
+%		y for the 2nd dimension, and z for the 3rd dimension.
25
+%
26
+%  new_elem_size (optional)  -  size of voxel along x y z direction for 
27
+%		a transformed 3D volume, or size of pixel along x y for
28
+%		a transformed 2D image. We assume x for the 1st dimension
29
+%		y for the 2nd dimension, and z for the 3rd dimension.
30
+%		'new_elem_size' is 1 if it is default or empty.
31
+%
32
+%		You can increase its value to decrease the resampling rate,
33
+%		and make the 2D image or 3D volume more coarse. It works
34
+%		just like 'interp3'.
35
+%
36
+%  verbose (optional) - 1, 0
37
+%		1:  show transforming progress in percentage
38
+%		2:  progress will not be displayed
39
+%		'verbose' is 1 if it is default or empty.
40
+%
41
+%  bg (optional)  -	background voxel intensity in any extra corner that
42
+%		is caused by the interpolation. 0 in most cases. If it is
43
+%		default or empty, 'bg' will be the average of two corner
44
+%		voxel intensities in original data.
45
+%
46
+%  method (optional)  -	1, 2, or 3
47
+%		1:  for Trilinear interpolation
48
+%		2:  for Nearest Neighbor interpolation
49
+%		3:  for Fischer's Bresenham interpolation
50
+%		'method' is 1 if it is default or empty.
51
+%
52
+%  new_img  -	transformed 2D image or 3D volume
53
+%
54
+%  new_M  -	transformed affine matrix
55
+%
56
+%  Example 1 (3D rotation):
57
+%	load mri.mat;   old_img = double(squeeze(D));
58
+%	old_M = [0.88 0.5 3 -90; -0.5 0.88 3 -126; 0 0 2 -72; 0 0 0 1];
59
+%	new_img = affine(old_img, old_M, 2);
60
+%	[x y z] = meshgrid(1:128,1:128,1:27);
61
+%	sz = size(new_img);
62
+%	[x1 y1 z1] = meshgrid(1:sz(2),1:sz(1),1:sz(3));
63
+%	figure; slice(x, y, z, old_img, 64, 64, 13.5);
64
+%	shading flat; colormap(map); view(-66, 66);
65
+%	figure; slice(x1, y1, z1, new_img, sz(1)/2, sz(2)/2, sz(3)/2);
66
+%	shading flat; colormap(map); view(-66, 66);
67
+%
68
+%  Example 2 (2D interpolation):
69
+%	load mri.mat;   old_img=D(:,:,1,13)';
70
+%	old_M = [1 0 0; 0 1 0; 0 0 1];
71
+%	new_img = affine(old_img, old_M, [.2 .4]);
72
+%	figure; image(old_img); colormap(map);
73
+%	figure; image(new_img); colormap(map);
74
+%
75
+%  This program is inspired by:
76
+%  SPM5 Software from Wellcome Trust Centre for Neuroimaging
77
+%	http://www.fil.ion.ucl.ac.uk/spm/software
78
+%  Fischer, J., A. del Rio (2004). A Fast Method for Applying Rigid
79
+%	Transformations to Volume Data, WSCG2004 Conference.
80
+%	http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
81
+%  
82
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
83
+%
84
+function [new_img, new_M] = affine(old_img, old_M, new_elem_size, verbose, bg, method)
85
+
86
+   if ~exist('old_img','var') | ~exist('old_M','var')
87
+      error('Usage: [new_img new_M] = affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);');
88
+   end
89
+
90
+   if ndims(old_img) == 3
91
+      if ~isequal(size(old_M),[4 4])
92
+         error('old_M should be a 4x4 affine matrix for 3D volume.');
93
+      end
94
+   elseif ndims(old_img) == 2
95
+      if ~isequal(size(old_M),[3 3])
96
+         error('old_M should be a 3x3 affine matrix for 2D image.');
97
+      end
98
+   else
99
+      error('old_img should be either 2D image or 3D volume.');
100
+   end
101
+
102
+   if ~exist('new_elem_size','var') | isempty(new_elem_size)
103
+      new_elem_size = [1 1 1];
104
+   elseif length(new_elem_size) < 2
105
+      new_elem_size = new_elem_size(1)*ones(1,3);
106
+   elseif length(new_elem_size) < 3
107
+      new_elem_size = [new_elem_size(:); 1]';
108
+   end
109
+
110
+   if ~exist('method','var') | isempty(method)
111
+      method = 1;
112
+   elseif ~exist('bresenham_line3d.m','file') & method == 3
113
+      error([char(10) char(10) 'Please download 3D Bresenham''s line generation program from:' char(10) char(10) 'http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21057' char(10) char(10) 'to test Fischer''s Bresenham interpolation method.' char(10) char(10)]);
114
+   end
115
+
116
+   %  Make compatible to MATLAB earlier than version 7 (R14), which
117
+   %  can only perform arithmetic on double data type
118
+   %
119
+   old_img = double(old_img);
120
+   old_dim = size(old_img);
121
+
122
+   if ~exist('bg','var') | isempty(bg)
123
+      bg = mean([old_img(1) old_img(end)]);
124
+   end
125
+
126
+   if ~exist('verbose','var') | isempty(verbose)
127
+      verbose = 1;
128
+   end
129
+
130
+   if ndims(old_img) == 2
131
+      old_dim(3) = 1;
132
+      old_M = old_M(:, [1 2 3 3]);
133
+      old_M = old_M([1 2 3 3], :);
134
+      old_M(3,:) = [0 0 1 0];
135
+      old_M(:,3) = [0 0 1 0]';
136
+   end
137
+
138
+   %  Vertices of img in voxel
139
+   %
140
+   XYZvox = [	1		1		1
141
+		1		1		old_dim(3)
142
+		1		old_dim(2)	1
143
+		1		old_dim(2)	old_dim(3)
144
+		old_dim(1)	1		1
145
+		old_dim(1)	1		old_dim(3)
146
+		old_dim(1)	old_dim(2)	1
147
+		old_dim(1)	old_dim(2)	old_dim(3)   ]';
148
+
149
+   old_R = old_M(1:3,1:3);
150
+   old_T = old_M(1:3,4);
151
+
152
+   %  Vertices of img in millimeter
153
+   %
154
+   XYZmm = old_R*(XYZvox-1) + repmat(old_T, [1, 8]);
155
+
156
+   %  Make scale of new_M according to new_elem_size
157
+   %
158
+   new_M = diag([new_elem_size 1]);
159
+
160
+   %  Make translation so minimum vertex is moved to [1,1,1]
161
+   %
162
+   new_M(1:3,4) = round( min(XYZmm,[],2) );
163
+
164
+   %  New dimensions will be the maximum vertices in XYZ direction (dim_vox)
165
+   %  i.e. compute   dim_vox   via   dim_mm = R*(dim_vox-1)+T
166
+   %  where, dim_mm = round(max(XYZmm,[],2));
167
+   %
168
+   new_dim = ceil(new_M(1:3,1:3) \ ( round(max(XYZmm,[],2))-new_M(1:3,4) )+1)';
169
+
170
+   %  Initialize new_img with new_dim
171
+   %
172
+   new_img = zeros(new_dim(1:3));
173
+
174
+   %  Mask out any changes from Z axis of transformed volume, since we
175
+   %  will traverse it voxel by voxel below. We will only apply unit
176
+   %  increment of mask_Z(3,4) to simulate the cursor movement
177
+   %
178
+   %  i.e. we will use   mask_Z * new_XYZvox   to replace   new_XYZvox
179
+   %
180
+   mask_Z = diag(ones(1,4));
181
+   mask_Z(3,3) = 0;
182
+
183
+   %  It will be easier to do the interpolation if we invert the process
184
+   %  by not traversing the original volume. Instead, we traverse the
185
+   %  transformed volume, and backproject each voxel in the transformed 
186
+   %  volume back into the original volume. If the backprojected voxel
187
+   %  in original volume is within its boundary, the intensity of that
188
+   %  voxel can be used by the cursor location in the transformed volume.
189
+   %
190
+   %  First, we traverse along Z axis of transformed volume voxel by voxel
191
+   %
192
+   for z = 1:new_dim(3)
193
+
194
+      if verbose & ~mod(z,10)
195
+         fprintf('%.2f percent is done.\n', 100*z/new_dim(3));
196
+      end
197
+
198
+      %  We need to find out the mapping from voxel in the transformed
199
+      %  volume (new_XYZvox) to voxel in the original volume (old_XYZvox)
200
+      %
201
+      %  The following equation works, because they all equal to XYZmm:
202
+      %  new_R*(new_XYZvox-1) + new_T  ==  old_R*(old_XYZvox-1) + old_T
203
+      %
204
+      %  We can use modified new_M1 & old_M1 to substitute new_M & old_M
205
+      %      new_M1 * new_XYZvox       ==       old_M1 * old_XYZvox
206
+      %
207
+      %  where: M1 = M;   M1(:,4) = M(:,4) - sum(M(:,1:3),2);
208
+      %  and:             M(:,4) == [T; 1] == sum(M1,2)
209
+      %
210
+      %  Therefore:   old_XYZvox = old_M1 \ new_M1 * new_XYZvox;
211
+      %
212
+      %  Since we are traverse Z axis, and   new_XYZvox   is replaced
213
+      %  by   mask_Z * new_XYZvox, the above formula can be rewritten
214
+      %  as:    old_XYZvox = old_M1 \ new_M1 * mask_Z * new_XYZvox;
215
+      %
216
+      %  i.e. we find the mapping from new_XYZvox to old_XYZvox:
217
+      %  M = old_M1 \ new_M1 * mask_Z;
218
+      %
219
+      %  First, compute modified old_M1 & new_M1
220
+      %
221
+      old_M1 = old_M;   old_M1(:,4) = old_M(:,4) - sum(old_M(:,1:3),2);
222
+      new_M1 = new_M;   new_M1(:,4) = new_M(:,4) - sum(new_M(:,1:3),2);
223
+
224
+      %  Then, apply unit increment of mask_Z(3,4) to simulate the
225
+      %  cursor movement
226
+      %
227
+      mask_Z(3,4) = z;
228
+
229
+      %  Here is the mapping from new_XYZvox to old_XYZvox
230
+      %
231
+      M = old_M1 \ new_M1 * mask_Z;
232
+
233
+      switch method
234
+      case 1
235
+         new_img(:,:,z) = trilinear(old_img, new_dim, old_dim, M, bg);
236
+      case 2
237
+         new_img(:,:,z) = nearest_neighbor(old_img, new_dim, old_dim, M, bg);
238
+      case 3
239
+         new_img(:,:,z) = bresenham(old_img, new_dim, old_dim, M, bg);
240
+      end
241
+
242
+   end;			% for z
243
+
244
+   if ndims(old_img) == 2
245
+      new_M(3,:) = [];
246
+      new_M(:,3) = [];
247
+   end
248
+
249
+   return;					% affine
250
+
251
+
252
+%--------------------------------------------------------------------
253
+function img_slice = trilinear(img, dim1, dim2, M, bg)
254
+
255
+   img_slice = zeros(dim1(1:2));
256
+   TINY = 5e-2;					% tolerance
257
+
258
+   %  Dimension of transformed 3D volume
259
+   %
260
+   xdim1 = dim1(1);
261
+   ydim1 = dim1(2);
262
+
263
+   %  Dimension of original 3D volume
264
+   %
265
+   xdim2 = dim2(1);
266
+   ydim2 = dim2(2);
267
+   zdim2 = dim2(3);
268
+
269
+   %  initialize new_Y accumulation
270
+   %
271
+   Y2X = 0;
272
+   Y2Y = 0;
273
+   Y2Z = 0;
274
+
275
+   for y = 1:ydim1
276
+
277
+      %  increment of new_Y accumulation
278
+      %
279
+      Y2X = Y2X + M(1,2);		% new_Y to old_X
280
+      Y2Y = Y2Y + M(2,2);		% new_Y to old_Y
281
+      Y2Z = Y2Z + M(3,2);		% new_Y to old_Z
282
+
283
+      %  backproject new_Y accumulation and translation to old_XYZ
284
+      %
285
+      old_X = Y2X + M(1,4);
286
+      old_Y = Y2Y + M(2,4);
287
+      old_Z = Y2Z + M(3,4);
288
+
289
+      for x = 1:xdim1
290
+
291
+         %  accumulate the increment of new_X, and apply it
292
+         %  to the backprojected old_XYZ
293
+         %
294
+         old_X = M(1,1) + old_X  ;
295
+         old_Y = M(2,1) + old_Y  ;
296
+         old_Z = M(3,1) + old_Z  ;
297
+
298
+         %  within boundary of original image
299
+         %
300
+         if (	old_X > 1-TINY & old_X < xdim2+TINY & ...
301
+		old_Y > 1-TINY & old_Y < ydim2+TINY & ...
302
+		old_Z > 1-TINY & old_Z < zdim2+TINY	)
303
+
304
+            %  Calculate distance of old_XYZ to its neighbors for
305
+            %  weighted intensity average
306
+            %
307
+            dx = old_X - floor(old_X);
308
+            dy = old_Y - floor(old_Y);
309
+            dz = old_Z - floor(old_Z);
310
+
311
+            x000 = floor(old_X);
312
+            x100 = x000 + 1;
313
+
314
+            if floor(old_X) < 1
315
+               x000 = 1;
316
+               x100 = x000;
317
+            elseif floor(old_X) > xdim2-1
318
+               x000 = xdim2;
319
+               x100 = x000;
320
+            end
321
+
322
+            x010 = x000;
323
+            x001 = x000;
324
+            x011 = x000;
325
+
326
+            x110 = x100;
327
+            x101 = x100;
328
+            x111 = x100;
329
+
330
+            y000 = floor(old_Y);
331
+            y010 = y000 + 1;
332
+
333
+            if floor(old_Y) < 1
334
+               y000 = 1;
335
+               y100 = y000;
336
+            elseif floor(old_Y) > ydim2-1
337
+               y000 = ydim2;
338
+               y010 = y000;
339
+            end
340
+
341
+            y100 = y000;
342
+            y001 = y000;
343
+            y101 = y000;
344
+
345
+            y110 = y010;
346
+            y011 = y010;
347
+            y111 = y010;
348
+
349
+            z000 = floor(old_Z);
350
+            z001 = z000 + 1;
351
+
352
+            if floor(old_Z) < 1
353
+               z000 = 1;
354
+               z001 = z000;
355
+            elseif floor(old_Z) > zdim2-1
356
+               z000 = zdim2;
357
+               z001 = z000;
358
+            end
359
+
360
+            z100 = z000;
361
+            z010 = z000;
362
+            z110 = z000;
363
+
364
+            z101 = z001;
365
+            z011 = z001;
366
+            z111 = z001;
367
+
368
+            x010 = x000;
369
+            x001 = x000;
370
+            x011 = x000;
371
+
372
+            x110 = x100;
373
+            x101 = x100;
374
+            x111 = x100;
375
+
376
+            v000 = double(img(x000, y000, z000));
377
+            v010 = double(img(x010, y010, z010));
378
+            v001 = double(img(x001, y001, z001));
379
+            v011 = double(img(x011, y011, z011));
380
+
381
+            v100 = double(img(x100, y100, z100));
382
+            v110 = double(img(x110, y110, z110));
383
+            v101 = double(img(x101, y101, z101));
384
+            v111 = double(img(x111, y111, z111));
385
+
386
+            img_slice(x,y) = v000*(1-dx)*(1-dy)*(1-dz) + ...
387
+               v010*(1-dx)*dy*(1-dz) + ...
388
+               v001*(1-dx)*(1-dy)*dz + ...
389
+               v011*(1-dx)*dy*dz + ...
390
+               v100*dx*(1-dy)*(1-dz) + ...
391
+               v110*dx*dy*(1-dz) + ...
392
+               v101*dx*(1-dy)*dz + ...
393
+               v111*dx*dy*dz;
394
+
395
+         else
396
+            img_slice(x,y) = bg;
397
+
398
+         end	% if boundary
399
+
400
+      end	% for x
401
+   end		% for y
402
+
403
+   return;					% trilinear
404
+
405
+
406
+%--------------------------------------------------------------------
407
+function img_slice = nearest_neighbor(img, dim1, dim2, M, bg)
408
+
409
+   img_slice = zeros(dim1(1:2));
410
+
411
+   %  Dimension of transformed 3D volume
412
+   %
413
+   xdim1 = dim1(1);
414
+   ydim1 = dim1(2);
415
+
416
+   %  Dimension of original 3D volume
417
+   %
418
+   xdim2 = dim2(1);
419
+   ydim2 = dim2(2);
420
+   zdim2 = dim2(3);
421
+
422
+   %  initialize new_Y accumulation
423
+   %
424
+   Y2X = 0;
425
+   Y2Y = 0;
426
+   Y2Z = 0;
427
+
428
+   for y = 1:ydim1
429
+
430
+      %  increment of new_Y accumulation
431
+      %
432
+      Y2X = Y2X + M(1,2);		% new_Y to old_X
433
+      Y2Y = Y2Y + M(2,2);		% new_Y to old_Y
434
+      Y2Z = Y2Z + M(3,2);		% new_Y to old_Z
435
+
436
+      %  backproject new_Y accumulation and translation to old_XYZ
437
+      %
438
+      old_X = Y2X + M(1,4);
439
+      old_Y = Y2Y + M(2,4);
440
+      old_Z = Y2Z + M(3,4);
441
+
442
+      for x = 1:xdim1
443
+
444
+         %  accumulate the increment of new_X and apply it
445
+         %  to the backprojected old_XYZ
446
+         %
447
+         old_X = M(1,1) + old_X  ;
448
+         old_Y = M(2,1) + old_Y  ;
449
+         old_Z = M(3,1) + old_Z  ;
450
+
451
+         xi = round(old_X);
452
+         yi = round(old_Y);
453
+         zi = round(old_Z);
454
+
455
+         %  within boundary of original image
456
+         %
457
+         if (	xi >= 1 & xi <= xdim2 & ...
458
+		yi >= 1 & yi <= ydim2 & ...
459
+		zi >= 1 & zi <= zdim2	)
460
+
461
+            img_slice(x,y) = img(xi,yi,zi);
462
+
463
+         else
464
+            img_slice(x,y) = bg;
465
+
466
+         end	% if boundary
467
+
468
+      end	% for x
469
+   end		% for y
470
+
471
+   return;					% nearest_neighbor
472
+
473
+
474
+%--------------------------------------------------------------------
475
+function img_slice = bresenham(img, dim1, dim2, M, bg)
476
+
477
+   img_slice = zeros(dim1(1:2));
478
+
479
+   %  Dimension of transformed 3D volume
480
+   %
481
+   xdim1 = dim1(1);
482
+   ydim1 = dim1(2);
483
+
484
+   %  Dimension of original 3D volume
485
+   %
486
+   xdim2 = dim2(1);
487
+   ydim2 = dim2(2);
488
+   zdim2 = dim2(3);
489
+
490
+   for y = 1:ydim1
491
+
492
+      start_old_XYZ = round(M*[0     y 0 1]');
493
+      end_old_XYZ   = round(M*[xdim1 y 0 1]');
494
+
495
+      [X Y Z] = bresenham_line3d(start_old_XYZ, end_old_XYZ);
496
+
497
+      %  line error correction
498
+      %
499
+%      del = end_old_XYZ - start_old_XYZ;
500
+ %     del_dom = max(del);
501
+  %    idx_dom = find(del==del_dom);
502
+   %   idx_dom = idx_dom(1);
503
+    %  idx_other = [1 2 3];
504
+     % idx_other(idx_dom) = [];
505
+      %del_x1 = del(idx_other(1));
506
+%      del_x2 = del(idx_other(2));
507
+ %     line_slope = sqrt((del_x1/del_dom)^2 + (del_x2/del_dom)^2 + 1);
508
+  %    line_error = line_slope - 1;
509
+% line error correction removed because it is too slow
510
+
511
+      for x = 1:xdim1
512
+
513
+         %  rescale ratio
514
+         %
515
+         i = round(x * length(X) / xdim1);
516
+
517
+         if i < 1
518
+            i = 1;
519
+         elseif i > length(X)
520
+            i = length(X);
521
+         end
522
+
523
+         xi = X(i);
524
+         yi = Y(i);
525
+         zi = Z(i);
526
+
527
+         %  within boundary of the old XYZ space
528
+         %
529
+         if (	xi >= 1 & xi <= xdim2 & ...
530
+		yi >= 1 & yi <= ydim2 & ...
531
+		zi >= 1 & zi <= zdim2	)
532
+
533
+            img_slice(x,y) = img(xi,yi,zi);
534
+
535
+%            if line_error > 1
536
+ %              x = x + 1;
537
+
538
+%               if x <= xdim1
539
+ %                 img_slice(x,y) = img(xi,yi,zi);
540
+  %                line_error = line_slope - 1;
541
+   %            end
542
+    %        end		% if line_error
543
+% line error correction removed because it is too slow
544
+
545
+         else
546
+            img_slice(x,y) = bg;
547
+
548
+         end	% if boundary
549
+
550
+      end	% for x
551
+   end		% for y
552
+
553
+   return;					% bresenham
554
+
0 555
new file mode 100644
... ...
@@ -0,0 +1,94 @@
1
+%BIPOLAR returns an M-by-3 matrix containing a blue-red colormap, in
2
+%	in which red stands for positive, blue stands for negative, 
3
+%	and white stands for 0.
4
+%
5
+%  Usage: cmap = bipolar(M, lo, hi, contrast);  or  cmap = bipolar;
6
+%
7
+%  cmap:  output M-by-3 matrix for BIPOLAR colormap.
8
+%  M:	  number of shades in the colormap. By default, it is the
9
+%	  same length as the current colormap.
10
+%  lo:	  the lowest value to represent.
11
+%  hi:	  the highest value to represent.
12
+%
13
+%  Inspired from the LORETA PASCAL program:
14
+%	http://www.unizh.ch/keyinst/NewLORETA
15
+%
16
+%  jimmy@rotman-baycrest.on.ca
17
+%
18
+%----------------------------------------------------------------
19
+function cmap = bipolar(M, lo, hi, contrast)
20
+
21
+   if ~exist('contrast','var')
22
+      contrast = 128;
23
+   end
24
+
25
+   if ~exist('lo','var')
26
+      lo = -1;
27
+   end
28
+
29
+   if ~exist('hi','var')
30
+      hi = 1;
31
+   end
32
+
33
+   if ~exist('M','var')
34
+      cmap = colormap;
35
+      M = size(cmap,1);
36
+   end
37
+
38
+   steepness = 10 ^ (1 - (contrast-1)/127);
39
+   pos_infs = 1e-99;
40
+   neg_infs = -1e-99;
41
+
42
+   doubleredc = [];
43
+   doublebluec = [];
44
+
45
+   if lo >= 0		% all positive
46
+
47
+      if lo == 0
48
+         lo = pos_infs;
49
+      end
50
+
51
+      for i=linspace(hi/M, hi, M)
52
+         t = exp(log(i/hi)*steepness);
53
+         doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
54
+      end
55
+
56
+      cmap = doubleredc;
57
+
58
+   elseif hi <= 0	% all negative
59
+
60
+      if hi == 0
61
+         hi = neg_infs;
62
+      end
63
+
64
+      for i=linspace(abs(lo)/M, abs(lo), M)
65
+         t = exp(log(i/abs(lo))*steepness);
66
+         doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
67
+      end
68
+
69
+      cmap = flipud(doublebluec);
70
+
71
+   else
72
+
73
+      if hi > abs(lo)
74
+         maxc = hi;
75
+      else
76
+         maxc = abs(lo);
77
+      end
78
+
79
+      for i=linspace(maxc/M, hi, round(M*hi/(hi-lo)))
80
+         t = exp(log(i/maxc)*steepness);
81
+         doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
82
+      end      
83
+
84
+      for i=linspace(maxc/M, abs(lo), round(M*abs(lo)/(hi-lo)))
85
+         t = exp(log(i/maxc)*steepness);
86
+         doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
87
+      end
88
+
89
+      cmap = [flipud(doublebluec); doubleredc];
90
+
91
+   end
92
+
93
+   return;					% bipolar
94
+
0 95
new file mode 100644
... ...
@@ -0,0 +1,189 @@
1
+%  Generate X Y Z coordinates of a 3D Bresenham's line between
2
+%  two given points.
3
+%
4
+%  A very useful application of this algorithm can be found in the
5
+%  implementation of Fischer's Bresenham interpolation method in my
6
+%  another program that can rotate three dimensional image volume
7
+%  with an affine matrix:
8
+%  http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21080
9
+%
10
+%  Usage: [X Y Z] = bresenham_line3d(P1, P2, [precision]);
11
+%
12
+%  P1	- vector for Point1, where P1 = [x1 y1 z1]
13
+%
14
+%  P2	- vector for Point2, where P2 = [x2 y2 z2]
15
+%
16
+%  precision (optional) - Although according to Bresenham's line
17
+%	algorithm, point coordinates x1 y1 z1 and x2 y2 z2 should
18
+%	be integer numbers, this program extends its limit to all
19
+%	real numbers. If any of them are floating numbers, you
20
+%	should specify how many digits of decimal that you would
21
+%	like to preserve. Be aware that the length of output X Y
22
+%	Z coordinates will increase in 10 times for each decimal
23
+%	digit that you want to preserve. By default, the precision
24
+%	is 0, which means that they will be rounded to the nearest
25
+%	integer.
26
+%
27
+%  X	- a set of x coordinates on Bresenham's line
28
+%
29
+%  Y	- a set of y coordinates on Bresenham's line
30
+%
31
+%  Z	- a set of z coordinates on Bresenham's line
32
+%
33
+%  Therefore, all points in XYZ set (i.e. P(i) = [X(i) Y(i) Z(i)])
34
+%  will constitute the Bresenham's line between P1 and P1.
35
+%
36
+%  Example:
37
+%	P1 = [12 37 6];     P2 = [46 3 35];
38
+%	[X Y Z] = bresenham_line3d(P1, P2);
39
+%	figure; plot3(X,Y,Z,'s','markerface','b');
40
+%
41
+%  This program is ported to MATLAB from:
42
+%
43
+%  B.Pendleton.  line3d - 3D Bresenham's (a 3D line drawing algorithm)
44
+%  ftp://ftp.isc.org/pub/usenet/comp.sources.unix/volume26/line3d, 1992
45
+%
46
+%  Which is also referenced by:
47
+%
48
+%  Fischer, J., A. del Rio (2004).  A Fast Method for Applying Rigid
49
+%  Transformations to Volume Data, WSCG2004 Conference.
50
+%  http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
51
+%
52
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
53
+%
54
+function [X,Y,Z] = bresenham_line3d(P1, P2, precision)
55
+
56
+   if ~exist('precision','var') | isempty(precision) | round(precision) == 0
57
+      precision = 0;
58
+      P1 = round(P1);
59
+      P2 = round(P2);
60
+   else
61
+      precision = round(precision);
62
+      P1 = round(P1*(10^precision));
63
+      P2 = round(P2*(10^precision));
64
+   end
65
+
66
+   d = max(abs(P2-P1)+1);
67
+   X = zeros(1, d);
68
+   Y = zeros(1, d);
69
+   Z = zeros(1, d);
70
+
71
+   x1 = P1(1);
72
+   y1 = P1(2);
73
+   z1 = P1(3);
74
+
75
+   x2 = P2(1);
76
+   y2 = P2(2);
77
+   z2 = P2(3);
78
+
79
+   dx = x2 - x1;
80
+   dy = y2 - y1;
81
+   dz = z2 - z1;
82
+
83
+   ax = abs(dx)*2;
84
+   ay = abs(dy)*2;
85
+   az = abs(dz)*2;
86
+
87
+   sx = sign(dx);
88
+   sy = sign(dy);
89
+   sz = sign(dz);
90
+
91
+   x = x1;
92
+   y = y1;
93
+   z = z1;
94
+   idx = 1;
95
+
96
+   if(ax>=max(ay,az))			% x dominant
97
+      yd = ay - ax/2;
98
+      zd = az - ax/2;
99
+
100
+      while(1)
101
+         X(idx) = x;
102
+         Y(idx) = y;
103
+         Z(idx) = z;
104
+         idx = idx + 1;
105
+
106
+         if(x == x2)		% end
107
+            break;
108
+         end
109
+
110
+         if(yd >= 0)		% move along y
111
+            y = y + sy;
112
+            yd = yd - ax;
113
+         end
114
+
115
+         if(zd >= 0)		% move along z
116
+            z = z + sz;
117
+            zd = zd - ax;
118
+         end
119
+
120
+         x  = x  + sx;		% move along x
121
+         yd = yd + ay;
122
+         zd = zd + az;
123
+      end
124
+   elseif(ay>=max(ax,az))		% y dominant
125
+      xd = ax - ay/2;
126
+      zd = az - ay/2;
127
+
128
+      while(1)
129
+         X(idx) = x;
130
+         Y(idx) = y;
131
+         Z(idx) = z;
132
+         idx = idx + 1;
133
+
134
+         if(y == y2)		% end
135
+            break;
136
+         end
137
+
138
+         if(xd >= 0)		% move along x
139
+            x = x + sx;
140
+            xd = xd - ay;
141
+         end
142
+
143
+         if(zd >= 0)		% move along z
144
+            z = z + sz;
145
+            zd = zd - ay;
146
+         end
147
+
148
+         y  = y  + sy;		% move along y
149
+         xd = xd + ax;
150
+         zd = zd + az;
151
+      end
152
+   elseif(az>=max(ax,ay))		% z dominant
153
+      xd = ax - az/2;
154
+      yd = ay - az/2;
155
+
156
+      while(1)
157
+         X(idx) = x;
158
+         Y(idx) = y;
159
+         Z(idx) = z;
160
+         idx = idx + 1;
161
+
162
+         if(z == z2)		% end
163
+            break;
164
+         end
165
+
166
+         if(xd >= 0)		% move along x
167
+            x = x + sx;
168
+            xd = xd - az;
169
+         end
170
+
171
+         if(yd >= 0)		% move along y
172
+            y = y + sy;
173
+            yd = yd - az;
174
+         end
175
+
176
+         z  = z  + sz;		% move along z
177
+         xd = xd + ax;
178
+         yd = yd + ay;
179
+      end
180
+   end
181
+
182
+   if precision ~= 0
183
+      X = X/(10^precision);
184
+      Y = Y/(10^precision);
185
+      Z = Z/(10^precision);
186
+   end
187
+
188
+   return;					% bresenham_line3d
189
+
0 190
new file mode 100644
... ...
@@ -0,0 +1,202 @@
1
+%  Collapse multiple single-scan NIFTI files into a multiple-scan NIFTI file
2
+%
3
+%  Usage: collapse_nii_scan(scan_file_pattern, [collapsed_filename], [scan_file_folder])
4
+%
5
+%  Here, scan_file_pattern should look like: 'myscan_0*.img'
6
+%  If collapsed_filename is omit, 'multi_scan.nii' will be used
7
+%  If scan_file_folder is omit, current file folder will be used
8
+%
9
+%  The order of volumes in the collapsed file will be the order of 
10
+%  corresponding filenames for those selected scan files.
11
+%
12
+%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
13
+%
14
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
15
+%
16
+function collapse_nii_scan(scan_pattern, fileprefix, scan_path)
17
+
18
+   if ~exist('fileprefix','var'), fileprefix = 'multi_scan.nii'; end
19
+   if ~exist('scan_path','var'), scan_path = pwd; end
20
+
21
+   filetype = 1;
22
+
23
+   %  Note: fileprefix is actually the filename you want to save
24
+   %   
25
+   if findstr('.nii',fileprefix)
26
+      filetype = 2;
27
+      fileprefix = strrep(fileprefix,'.nii','');
28
+   end
29
+   
30
+   if findstr('.hdr',fileprefix)
31
+      fileprefix = strrep(fileprefix,'.hdr','');
32
+   end
33
+   
34
+   if findstr('.img',fileprefix)
35
+      fileprefix = strrep(fileprefix,'.img','');
36
+   end
37
+
38
+   pnfn = fullfile(scan_path, scan_pattern);
39
+
40
+   file_lst = dir(pnfn);
41
+   flist = {file_lst.name};
42
+   flist = flist(:);
43
+
44
+   nii = load_nii(flist{1});
45
+   nii.hdr.dime.dim(5) = length(flist);
46
+   hdr = nii.hdr;
47
+
48
+   if isfield(nii,'ext') & ~isempty(nii.ext)
49
+      ext = nii.ext;
50
+      [ext, esize_total] = verify_nii_ext(ext);
51
+   else
52
+      ext = [];
53
+   end
54
+
55
+   switch double(hdr.dime.datatype),
56
+   case   1,
57
+      hdr.dime.bitpix = int16(1 ); precision = 'ubit1';
58
+   case   2,
59
+      hdr.dime.bitpix = int16(8 ); precision = 'uint8';
60
+   case   4,
61
+      hdr.dime.bitpix = int16(16); precision = 'int16';
62
+   case   8,
63
+      hdr.dime.bitpix = int16(32); precision = 'int32';
64
+   case  16,
65
+      hdr.dime.bitpix = int16(32); precision = 'float32';
66
+   case  32,
67
+      hdr.dime.bitpix = int16(64); precision = 'float32';
68
+   case  64,
69
+      hdr.dime.bitpix = int16(64); precision = 'float64';
70
+   case 128,
71
+      hdr.dime.bitpix = int16(24); precision = 'uint8';
72
+   case 256 
73
+      hdr.dime.bitpix = int16(8 ); precision = 'int8';
74
+   case 512 
75
+      hdr.dime.bitpix = int16(16); precision = 'uint16';
76
+   case 768 
77
+      hdr.dime.bitpix = int16(32); precision = 'uint32';
78
+   case 1024
79
+      hdr.dime.bitpix = int16(64); precision = 'int64';
80
+   case 1280
81
+      hdr.dime.bitpix = int16(64); precision = 'uint64';
82
+   case 1792,
83
+      hdr.dime.bitpix = int16(128); precision = 'float64';
84
+   otherwise
85
+      error('This datatype is not supported');
86
+   end
87
+
88
+   if filetype == 2
89
+      fid = fopen(sprintf('%s.nii',fileprefix),'w');
90
+      
91
+      if fid < 0,
92
+         msg = sprintf('Cannot open file %s.nii.',fileprefix);
93
+         error(msg);
94
+      end
95
+      
96
+      hdr.dime.vox_offset = 352;
97
+
98
+      if ~isempty(ext)
99
+         hdr.dime.vox_offset = hdr.dime.vox_offset + esize_total;
100
+      end
101
+
102
+      hdr.hist.magic = 'n+1';
103
+      save_nii_hdr(hdr, fid);
104
+
105
+      if ~isempty(ext)
106
+         save_nii_ext(ext, fid);
107
+      end
108
+   else
109
+      fid = fopen(sprintf('%s.hdr',fileprefix),'w');
110
+      
111
+      if fid < 0,
112
+         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
113
+         error(msg);
114
+      end
115
+      
116
+      hdr.dime.vox_offset = 0;
117
+      hdr.hist.magic = 'ni1';
118
+      save_nii_hdr(hdr, fid);
119
+
120
+      if ~isempty(ext)
121
+         save_nii_ext(ext, fid);
122
+      end
123
+      
124
+      fclose(fid);
125
+      fid = fopen(sprintf('%s.img',fileprefix),'w');
126
+   end
127
+
128
+   if filetype == 2 & isempty(ext)
129
+      skip_bytes = double(hdr.dime.vox_offset) - 348;
130
+   else
131
+      skip_bytes = 0;
132
+   end
133
+
134
+   glmax = -inf;
135
+   glmin = inf;
136
+
137
+   for i = 1:length(flist)
138
+      nii = load_nii(flist{i});
139
+
140
+      if double(hdr.dime.datatype) == 128
141
+
142
+         %  RGB planes are expected to be in the 4th dimension of nii.img
143
+         %
144
+         if(size(nii.img,4)~=3)
145
+            error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
146
+         end
147
+
148
+         if old_RGB
149
+            nii.img = permute(nii.img, [1 2 4 3 5 6 7 8]);
150
+         else
151
+            nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]);
152
+         end
153
+      end
154
+
155
+      %  For complex float32 or complex float64, voxel values
156
+      %  include [real, imag]
157
+      %
158
+      if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
159
+         real_img = real(nii.img(:))';
160
+         nii.img = imag(nii.img(:))';
161
+         nii.img = [real_img; nii.img];
162
+      end
163
+
164
+      if nii.hdr.dime.glmax > glmax
165
+         glmax = nii.hdr.dime.glmax;
166
+      end
167
+
168
+      if nii.hdr.dime.glmin < glmin
169
+         glmin = nii.hdr.dime.glmin;
170
+      end
171
+
172
+      fwrite(fid, nii.img, precision);
173
+   end
174
+
175
+   hdr.dime.glmax = round(glmax);
176
+   hdr.dime.glmin = round(glmin);
177
+
178
+   if filetype == 2
179
+      fseek(fid, 140, 'bof');
180
+      fwrite(fid, hdr.dime.glmax, 'int32');
181
+      fwrite(fid, hdr.dime.glmin, 'int32');
182
+   else
183
+      fid2 = fopen(sprintf('%s.hdr',fileprefix),'w');
184
+
185
+      if fid2 < 0,
186
+         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
187
+         error(msg);
188
+      end
189
+
190
+      save_nii_hdr(hdr, fid2);
191
+
192
+      if ~isempty(ext)
193
+         save_nii_ext(ext, fid2);
194
+      end
195
+
196
+      fclose(fid2);
197
+   end
198
+
199
+   fclose(fid);
200
+
201
+   return;					% collapse_nii_scan
202
+
0 203
new file mode 100644
... ...
@@ -0,0 +1,120 @@
1
+
2
+- Examples to load, make and save a nii struct:
3
+
4
+	To load Analyze data or NIFTI data to a structure:
5
+
6
+		nii = load_nii(NIFTI_file_name, [img_idx], [old_RGB24]);
7
+
8
+	img_idx is a numerical array of image indices along the temporal
9
+	axis, which is only available in NIFTI data. After you specify 
10
+	img_idx, only those images indexed by img_idx will be loaded. If
11
+	there is no img_idx or img_idx is empty, all available images 
12
+	will be loaded.
13
+
14
+	For RGB image, most people use RGB triple sequentially for each
15
+	voxel, like [R1 G1 B1 R2 G2 B2 ...]. However, some program like
16
+	Analyze 6.0 developed by AnalyzeDirect uses old RGB24, in a way
17
+	like [R1 R2 ... G1 G2 ... B1 B2 ...] for each slices. In this
18
+	case, you can set old_RGB24 flag to 1 and load data correctly:
19
+
20
+		nii = load_nii(NIFTI_file_name, [], 1);
21
+
22
+	To get a total number of images along the temporal axis:
23
+
24
+		num_scan = get_nii_frame(NIFTI_file_name);
25
+
26
+	You can also load the header extension if it exists:
27
+
28
+		nii.ext = load_nii_ext(NIFTI_file_name);
29
+
30
+	You can just load the Analyze or NIFTI header:
31
+	(header contains: hk, dime, and hist)
32
+
33
+		hdr = load_nii_hdr(NIFTI_file_name);
34
+
35
+	You can also save the structure to a new file:
36
+	(header extension will be saved if there is nii.ext structure)
37
+
38
+		save_nii(nii, NIFTI_file_name);
39
+
40
+	To make the structure from any 3D (or 4D) data:
41
+
42
+		img = rand(91,109,91); or
43
+		img = rand(64,64,21,18);
44
+		nii = make_nii(img [, voxel_size, origin, datatype] );
45
+
46
+	Use "help load_nii", "help save_nii", "help make_nii" etc.
47
+	to get more detail information.
48
+
49
+
50
+- Examples to plot a nii struct:
51
+  (More detail descriptions are available on top of "view_nii.m")
52
+
53
+	Simple way to plot a nii struct:
54
+
55
+		view_nii(nii);
56
+
57
+	The default colormap will use the Gray if all data values are
58
+	non-negative; otherwise, the default colormap will use BiPolar.
59
+	You can choose other colormap, including customized colormap
60
+	from panel.
61
+
62
+	To imbed the plot into your existing figure:
63
+
64
+		h = gcf;
65
+		opt.command = 'init';
66
+		opt.setarea = [0.3 0.1 0.6 0.8];
67
+		view_nii(h, nii, opt);
68
+
69
+	To add a colorbar:
70
+
71
+		opt.usecolorbar = 1;
72
+		view_nii(gcf, opt);
73
+
74
+		Here, opt.command is implicitly set to 'update'.
75
+
76
+	To display in real aspect ratio:
77
+
78
+		opt.usestretch = 0;
79
+		view_nii(gcf, opt);
80
+
81
+	If you want the data value to be directly used as the index
82
+	of colormap, instead of scale to the whole colormap:
83
+
84
+		opt.useimagesc = 0;
85
+		view_nii(gcf, opt);
86
+
87
+	If you modified the data value without changing the dimension,
88
+	voxel_size, and origin, you can update the display by:
89
+
90
+		opt.command = 'updateimg';
91
+		view_nii(gcf, nii.img, opt);
92
+
93
+	If the data is completely different, display can be updated by:
94
+
95
+		opt.command = 'updatenii';
96
+		view_nii(gcf, nii, opt);
97
+
98
+
99
+- Contrast and Brightness are available under Gray and Bipolar colormap:
100
+
101
+
102
+	Increase contrast in Gray colormap will make high end values
103
+	more distinguishable by sacrificing the low end values; The
104
+	minimum contrast (default) will display the whole range.
105
+
106
+	Increase or decrease contrast in BiPolar colormap will shift
107
+	the distinguishable position for both positive and negative
108
+	values.
109
+
110
+	Increase or decrease brightness in Gray colormap will shift
111
+	the distinguishable position.
112
+
113
+	Increase or decrease brightness in BiPolar colormap will make
114
+	both positive and negative values more distinguishable.
115
+
116
+
117
+- Required files:
118
+
119
+	All files in this package.
120
+
0 121
new file mode 100644
... ...
@@ -0,0 +1,24 @@
1
+%  Expand a multiple-scan NIFTI file into multiple single-scan NIFTI files
2
+%
3
+%  Usage: expand_nii_scan(multi_scan_filename, [img_idx], [path_to_save])
4
+%
5
+%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
6
+%
7
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
8
+%
9
+function expand_nii_scan(filename, img_idx, newpath)
10
+
11
+   if ~exist('newpath','var'), newpath = pwd; end
12
+   if ~exist('img_idx','var'), img_idx = 1:get_nii_frame(filename); end
13
+
14
+   for i=img_idx
15
+      nii_i = load_nii(filename, i);
16
+
17
+      fn = [nii_i.fileprefix '_' sprintf('%04d',i)];
18
+      pnfn = fullfile(newpath, fn);
19
+
20
+      save_nii(nii_i, pnfn);
21
+   end
22
+
23
+   return;					% expand_nii_scan
24
+
0 25
new file mode 100644
... ...
@@ -0,0 +1,255 @@
1
+%  Decode extra NIFTI header information into hdr.extra
2
+%
3
+%  Usage: hdr = extra_nii_hdr(hdr)
4
+%
5
+%  hdr can be obtained from load_nii_hdr
6
+%
7
+%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
8
+%
9
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
10
+%
11
+function hdr = extra_nii_hdr(hdr)
12
+
13
+   switch hdr.dime.datatype
14
+   case 1
15
+      extra.NIFTI_DATATYPES = 'DT_BINARY';
16
+   case 2
17
+      extra.NIFTI_DATATYPES = 'DT_UINT8';
18
+   case 4
19
+      extra.NIFTI_DATATYPES = 'DT_INT16';
20
+   case 8
21
+      extra.NIFTI_DATATYPES = 'DT_INT32';
22
+   case 16
23
+      extra.NIFTI_DATATYPES = 'DT_FLOAT32';
24
+   case 32
25
+      extra.NIFTI_DATATYPES = 'DT_COMPLEX64';
26
+   case 64
27
+      extra.NIFTI_DATATYPES = 'DT_FLOAT64';
28
+   case 128
29
+      extra.NIFTI_DATATYPES = 'DT_RGB24';
30
+   case 256
31
+      extra.NIFTI_DATATYPES = 'DT_INT8';
32
+   case 512
33
+      extra.NIFTI_DATATYPES = 'DT_UINT16';
34
+   case 768
35
+      extra.NIFTI_DATATYPES = 'DT_UINT32';
36
+   case 1024
37
+      extra.NIFTI_DATATYPES = 'DT_INT64';
38
+   case 1280
39
+      extra.NIFTI_DATATYPES = 'DT_UINT64';
40
+   case 1536
41
+      extra.NIFTI_DATATYPES = 'DT_FLOAT128';
42
+   case 1792
43
+      extra.NIFTI_DATATYPES = 'DT_COMPLEX128';
44
+   case 2048
45
+      extra.NIFTI_DATATYPES = 'DT_COMPLEX256';
46
+   otherwise
47
+      extra.NIFTI_DATATYPES = 'DT_UNKNOWN';
48
+   end
49
+
50
+   switch hdr.dime.intent_code
51
+   case 2
52
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CORREL';
53
+   case 3
54
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TTEST';
55
+   case 4
56
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_FTEST';
57
+   case 5
58
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_ZSCORE';
59
+   case 6
60
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHISQ';
61
+   case 7
62
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_BETA';
63
+   case 8
64
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_BINOM';
65
+   case 9
66
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_GAMMA';
67
+   case 10
68
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_POISSON';
69
+   case 11
70
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NORMAL';
71
+   case 12
72
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_FTEST_NONC';
73
+   case 13
74
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHISQ_NONC';
75
+   case 14
76
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOGISTIC';
77
+   case 15
78
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LAPLACE';
79
+   case 16
80
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_UNIFORM';
81
+   case 17
82
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TTEST_NONC';
83
+   case 18
84
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_WEIBULL';
85
+   case 19
86
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHI';
87
+   case 20
88
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_INVGAUSS';
89
+   case 21
90
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_EXTVAL';
91
+   case 22
92
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_PVAL';
93
+   case 23
94
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOGPVAL';
95
+   case 24
96
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOG10PVAL';
97
+   case 1001
98
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_ESTIMATE';
99
+   case 1002
100
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LABEL';
101
+   case 1003
102
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NEURONAME';
103
+   case 1004
104
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_GENMATRIX';
105
+   case 1005
106
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_SYMMATRIX';
107
+   case 1006
108
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_DISPVECT';
109
+   case 1007
110
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_VECTOR';
111
+   case 1008
112
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_POINTSET';
113
+   case 1009
114
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TRIANGLE';
115
+   case 1010
116
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_QUATERNION';
117
+   case 1011
118
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_DIMLESS';
119
+   otherwise
120
+      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NONE';
121
+   end
122
+
123
+   extra.NIFTI_INTENT_NAMES = hdr.hist.intent_name;
124
+
125
+   if hdr.hist.sform_code > 0
126
+      switch hdr.hist.sform_code
127
+      case 1
128
+         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_SCANNER_ANAT';
129
+      case 2
130
+         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_ALIGNED_ANAT';
131
+      case 3
132
+         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_TALAIRACH';
133
+      case 4
134
+         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_MNI_152';
135
+      otherwise
136
+         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
137
+      end
138
+
139
+      extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
140
+   elseif hdr.hist.qform_code > 0
141
+      extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
142
+
143
+      switch hdr.hist.qform_code
144
+      case 1
145
+         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_SCANNER_ANAT';
146
+      case 2
147
+         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_ALIGNED_ANAT';
148
+      case 3
149
+         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_TALAIRACH';
150
+      case 4
151
+         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_MNI_152';
152
+      otherwise
153
+         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
154
+      end
155
+   else
156
+      extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
157
+      extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
158
+   end
159
+
160
+   switch bitand(hdr.dime.xyzt_units, 7)	% mask with 0x07
161
+   case 1
162
+      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_METER';
163
+   case 2
164
+      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_MM';	% millimeter
165
+   case 3
166
+      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_MICRO';
167
+   otherwise
168
+      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_UNKNOWN';
169
+   end
170
+
171
+   switch bitand(hdr.dime.xyzt_units, 56)	% mask with 0x38
172
+   case 8
173
+      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_SEC';
174
+   case 16
175
+      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_MSEC';
176
+   case 24
177
+      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_USEC';	% microsecond
178
+   otherwise
179
+      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_UNKNOWN';
180
+   end
181
+
182
+   switch hdr.dime.xyzt_units
183
+   case 32
184
+      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_HZ';
185
+   case 40
186
+      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_PPM';	% part per million
187
+   case 48
188
+      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_RADS';	% radians per second
189
+   otherwise
190
+      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_UNKNOWN';
191
+   end
192
+
193
+   %  MRI-specific spatial and temporal information
194
+   %
195
+   dim_info = hdr.hk.dim_info;
196
+   extra.NIFTI_FREQ_DIM = bitand(dim_info, 3);
197
+   extra.NIFTI_PHASE_DIM = bitand(bitshift(dim_info, -2), 3);
198
+   extra.NIFTI_SLICE_DIM = bitand(bitshift(dim_info, -4), 3);
199
+
200
+   %  Check slice code
201
+   %
202
+   switch hdr.dime.slice_code
203
+   case 1
204
+      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_SEQ_INC';	% sequential increasing
205
+   case 2
206
+      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_SEQ_DEC';	% sequential decreasing
207
+   case 3
208
+      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_INC';	% alternating increasing
209
+   case 4
210
+      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_DEC';	% alternating decreasing