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 on30/03/2009 17:54:25
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,385 @@
1
+%  internal function
2
+
3
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
4
+
5
+function [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)
6
+
7
+   if ~exist('hdr','var') | ~exist('filetype','var') | ~exist('fileprefix','var') | ~exist('machine','var')
8
+      error('Usage: [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,[img_idx],[dim5_idx],[dim6_idx],[dim7_idx],[old_RGB]);');
9
+   end
10
+
11
+   if ~exist('img_idx','var') | isempty(img_idx) | hdr.dime.dim(5)<1
12
+      img_idx = [];
13
+   end
14
+
15
+   if ~exist('dim5_idx','var') | isempty(dim5_idx) | hdr.dime.dim(6)<1
16
+      dim5_idx = [];
17
+   end
18
+
19
+   if ~exist('dim6_idx','var') | isempty(dim6_idx) | hdr.dime.dim(7)<1
20
+      dim6_idx = [];
21
+   end
22
+
23
+   if ~exist('dim7_idx','var') | isempty(dim7_idx) | hdr.dime.dim(8)<1
24
+      dim7_idx = [];
25
+   end
26
+
27
+   if ~exist('old_RGB','var') | isempty(old_RGB)
28
+      old_RGB = 0;
29
+   end
30
+
31
+   %  check img_idx
32
+   %
33
+   if ~isempty(img_idx) & ~isnumeric(img_idx)
34
+      error('"img_idx" should be a numerical array.');
35
+   end
36
+
37
+   if length(unique(img_idx)) ~= length(img_idx)
38
+      error('Duplicate image index in "img_idx"');
39
+   end
40
+
41
+   if ~isempty(img_idx) & (min(img_idx) < 1 | max(img_idx) > hdr.dime.dim(5))
42
+      max_range = hdr.dime.dim(5);
43
+
44
+      if max_range == 1
45
+         error(['"img_idx" should be 1.']);
46
+      else
47
+         range = ['1 ' num2str(max_range)];
48
+         error(['"img_idx" should be an integer within the range of [' range '].']);
49
+      end
50
+   end
51
+
52
+   %  check dim5_idx
53
+   %
54
+   if ~isempty(dim5_idx) & ~isnumeric(dim5_idx)
55
+      error('"dim5_idx" should be a numerical array.');
56
+   end
57
+
58
+   if length(unique(dim5_idx)) ~= length(dim5_idx)
59
+      error('Duplicate image index in "dim5_idx"');
60
+   end
61
+
62
+   if ~isempty(dim5_idx) & (min(dim5_idx) < 1 | max(dim5_idx) > hdr.dime.dim(6))
63
+      max_range = hdr.dime.dim(6);
64
+
65
+      if max_range == 1
66
+         error(['"dim5_idx" should be 1.']);
67
+      else
68
+         range = ['1 ' num2str(max_range)];
69
+         error(['"dim5_idx" should be an integer within the range of [' range '].']);
70
+      end
71
+   end
72
+
73
+   %  check dim6_idx
74
+   %
75
+   if ~isempty(dim6_idx) & ~isnumeric(dim6_idx)
76
+      error('"dim6_idx" should be a numerical array.');
77
+   end
78
+
79
+   if length(unique(dim6_idx)) ~= length(dim6_idx)
80
+      error('Duplicate image index in "dim6_idx"');
81
+   end
82
+
83
+   if ~isempty(dim6_idx) & (min(dim6_idx) < 1 | max(dim6_idx) > hdr.dime.dim(7))
84
+      max_range = hdr.dime.dim(7);
85
+
86
+      if max_range == 1
87
+         error(['"dim6_idx" should be 1.']);
88
+      else
89
+         range = ['1 ' num2str(max_range)];
90
+         error(['"dim6_idx" should be an integer within the range of [' range '].']);
91
+      end
92
+   end
93
+
94
+   %  check dim7_idx
95
+   %
96
+   if ~isempty(dim7_idx) & ~isnumeric(dim7_idx)
97
+      error('"dim7_idx" should be a numerical array.');
98
+   end
99
+
100
+   if length(unique(dim7_idx)) ~= length(dim7_idx)
101
+      error('Duplicate image index in "dim7_idx"');
102
+   end
103
+
104
+   if ~isempty(dim7_idx) & (min(dim7_idx) < 1 | max(dim7_idx) > hdr.dime.dim(8))
105
+      max_range = hdr.dime.dim(8);
106
+
107
+      if max_range == 1
108
+         error(['"dim7_idx" should be 1.']);
109
+      else
110
+         range = ['1 ' num2str(max_range)];
111
+         error(['"dim7_idx" should be an integer within the range of [' range '].']);
112
+      end
113
+   end
114
+
115
+   [img,hdr] = read_image(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB);
116
+
117
+   return					% load_nii_img
118
+
119
+
120
+%---------------------------------------------------------------------
121
+function [img,hdr] = read_image(hdr, filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)
122
+
123
+   switch filetype
124
+   case {0, 1}
125
+      fn = [fileprefix '.img'];
126
+   case 2
127
+      fn = [fileprefix '.nii'];
128
+   end
129
+
130
+   fid = fopen(fn,'r',machine);
131
+
132
+   if fid < 0,
133
+      msg = sprintf('Cannot open file %s.',fn);
134
+      error(msg);
135
+   end
136
+
137
+   %  Set bitpix according to datatype
138
+   %
139
+   %  /*Acceptable values for datatype are*/ 
140
+   %
141
+   %     0 None                     (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN 
142
+   %     1 Binary                         (ubit1, bitpix=1) % DT_BINARY 
143
+   %     2 Unsigned char         (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8 
144
+   %     4 Signed short                  (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16 
145
+   %     8 Signed integer                (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32 
146
+   %    16 Floating point    (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32 
147
+   %    32 Complex, 2 float32      (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
148
+   %    64 Double precision  (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64 
149
+   %   128 uint8 RGB                 (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24 
150
+   %   256 Signed char            (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8 
151
+   %   511 Single RGB              (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96
152
+   %   512 Unsigned short               (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16 
153
+   %   768 Unsigned integer             (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32 
154
+   %  1024 Signed long long              (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64
155
+   %  1280 Unsigned long long           (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64 
156
+   %  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128 
157
+   %  1792 Complex128, 2 float64  (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128 
158
+   %  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128 
159
+   %
160
+   switch hdr.dime.datatype
161
+   case   1,
162
+      hdr.dime.bitpix = 1;  precision = 'ubit1';
163
+   case   2,
164
+      hdr.dime.bitpix = 8;  precision = 'uint8';
165
+   case   4,
166
+      hdr.dime.bitpix = 16; precision = 'int16';
167
+   case   8,
168
+      hdr.dime.bitpix = 32; precision = 'int32';
169
+   case  16,
170
+      hdr.dime.bitpix = 32; precision = 'float32';
171
+   case  32,
172
+      hdr.dime.bitpix = 64; precision = 'float32';
173
+   case  64,
174
+      hdr.dime.bitpix = 64; precision = 'float64';
175
+   case 128,
176
+      hdr.dime.bitpix = 24; precision = 'uint8';
177
+   case 256 
178
+      hdr.dime.bitpix = 8;  precision = 'int8';
179
+   case 511 
180
+      hdr.dime.bitpix = 96; precision = 'float32';
181
+   case 512 
182
+      hdr.dime.bitpix = 16; precision = 'uint16';
183
+   case 768 
184
+      hdr.dime.bitpix = 32; precision = 'uint32';
185
+   case 1024
186
+      hdr.dime.bitpix = 64; precision = 'int64';
187
+   case 1280
188
+      hdr.dime.bitpix = 64; precision = 'uint64';
189
+   case 1792,
190
+      hdr.dime.bitpix = 128; precision = 'float64';
191
+   otherwise
192
+      error('This datatype is not supported'); 
193
+   end
194
+
195
+   hdr.dime.dim(find(hdr.dime.dim < 1)) = 1;
196
+
197
+   %  move pointer to the start of image block
198
+   %
199
+   switch filetype
200
+   case {0, 1}
201
+      fseek(fid, 0, 'bof');
202
+   case 2
203
+      fseek(fid, hdr.dime.vox_offset, 'bof');
204
+   end
205
+
206
+   %  Load whole image block for old Analyze format or binary image;
207
+   %  otherwise, load images that are specified in img_idx, dim5_idx,
208
+   %  dim6_idx, and dim7_idx
209
+   %
210
+   %  For binary image, we have to read all because pos can not be
211
+   %  seeked in bit and can not be calculated the way below.
212
+   %
213
+   if filetype == 0 | hdr.dime.datatype == 1 | isequal(hdr.dime.dim(5:8),ones(1,4)) | ...
214
+	(isempty(img_idx) & isempty(dim5_idx) & isempty(dim6_idx) & isempty(dim7_idx))
215
+
216
+      %  For each frame, precision of value will be read 
217
+      %  in img_siz times, where img_siz is only the 
218
+      %  dimension size of an image, not the byte storage
219
+      %  size of an image.
220
+      %
221
+      img_siz = prod(hdr.dime.dim(2:8));
222
+
223
+      %  For complex float32 or complex float64, voxel values
224
+      %  include [real, imag]
225
+      %
226
+      if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
227
+         img_siz = img_siz * 2;
228
+      end
229
+	 
230
+      %MPH: For RGB24, voxel values include 3 separate color planes
231
+      %
232
+      if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
233
+	 img_siz = img_siz * 3;
234
+      end
235
+
236
+      img = fread(fid, img_siz, sprintf('*%s',precision));
237
+
238
+      d1 = hdr.dime.dim(2);
239
+      d2 = hdr.dime.dim(3);
240
+      d3 = hdr.dime.dim(4);
241
+      d4 = hdr.dime.dim(5);
242
+      d5 = hdr.dime.dim(6);
243
+      d6 = hdr.dime.dim(7);
244
+      d7 = hdr.dime.dim(8);
245
+
246
+      if isempty(img_idx)
247
+         img_idx = 1:d4;
248
+      end
249
+
250
+      if isempty(dim5_idx)
251
+         dim5_idx = 1:d5;
252
+      end
253
+
254
+      if isempty(dim6_idx)
255
+         dim6_idx = 1:d6;
256
+      end
257
+
258
+      if isempty(dim7_idx)
259
+         dim7_idx = 1:d7;
260
+      end
261
+   else
262
+
263
+      img = [];
264
+
265
+      %  For each frame, precision of value will be read 
266
+      %  in img_siz times, where img_siz is only the 
267
+      %  dimension size of an image, not the byte storage
268
+      %  size of an image.
269
+      %
270
+      img_siz = prod(hdr.dime.dim(2:4));
271
+
272
+      d1 = hdr.dime.dim(2);
273
+      d2 = hdr.dime.dim(3);
274
+      d3 = hdr.dime.dim(4);
275
+      d4 = hdr.dime.dim(5);
276
+      d5 = hdr.dime.dim(6);
277
+      d6 = hdr.dime.dim(7);
278
+      d7 = hdr.dime.dim(8);
279
+
280
+      if isempty(img_idx)
281
+         img_idx = 1:d4;
282
+      end
283
+
284
+      if isempty(dim5_idx)
285
+         dim5_idx = 1:d5;
286
+      end
287
+
288
+      if isempty(dim6_idx)
289
+         dim6_idx = 1:d6;
290
+      end
291
+
292
+      if isempty(dim7_idx)
293
+         dim7_idx = 1:d7;
294
+      end
295
+
296
+      for i7=1:length(dim7_idx)
297
+         for i6=1:length(dim6_idx)
298
+            for i5=1:length(dim5_idx)
299
+               for t=1:length(img_idx)
300
+
301
+                  %  Position is seeked in bytes. To convert dimension size
302
+                  %  to byte storage size, hdr.dime.bitpix/8 will be
303
+                  %  applied.
304
+                  %
305
+                  pos = sub2ind([d1 d2 d3 d4 d5 d6 d7], 1, 1, 1, ...
306
+			img_idx(t), dim5_idx(i5),dim6_idx(i6),dim7_idx(i7)) -1;
307
+                  pos = pos * hdr.dime.bitpix/8;
308
+
309
+                  %  For complex float32 or complex float64, voxel values
310
+                  %  include [real, imag]
311
+                  %
312
+                  if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
313
+                     img_siz = img_siz * 2;
314
+                  end
315
+
316
+                  %MPH: For RGB24, voxel values include 3 separate color planes
317
+                  %
318
+                  if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
319
+	             img_siz = img_siz * 3;
320
+                  end
321
+         
322
+                  if filetype == 2
323
+                     fseek(fid, pos + hdr.dime.vox_offset, 'bof');
324
+                  else
325
+                     fseek(fid, pos, 'bof');
326
+                  end
327
+
328
+                  %  For each frame, fread will read precision of value
329
+                  %  in img_siz times
330
+                  %
331
+                  img = [img fread(fid, img_siz, sprintf('*%s',precision))];
332
+               end
333
+            end
334
+         end
335
+      end
336
+   end
337
+
338
+   %  For complex float32 or complex float64, voxel values
339
+   %  include [real, imag]
340
+   %
341
+   if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
342
+      img = reshape(img, [2, length(img)/2]);
343
+      img = complex(img(1,:)', img(2,:)');
344
+   end
345
+
346
+   fclose(fid);
347
+
348
+   %  Update the global min and max values 
349
+   %
350
+   hdr.dime.glmax = max(double(img(:)));
351
+   hdr.dime.glmin = min(double(img(:)));
352
+
353
+   if old_RGB & hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
354
+      img = squeeze(reshape(img, [hdr.dime.dim(2:3) 3 hdr.dime.dim(4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
355
+      img = permute(img, [1 2 4 3 5 6 7 8]);
356
+   elseif hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
357
+      img = squeeze(reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
358
+      img = permute(img, [2 3 4 1 5 6 7 8]);
359
+   elseif hdr.dime.datatype == 511 & hdr.dime.bitpix == 96
360
+      img = double(img);
361
+      img = (img - min(img))/(max(img) - min(img));
362
+      img = squeeze(reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
363
+      img = permute(img, [2 3 4 1 5 6 7 8]);
364
+   else
365
+      img = squeeze(reshape(img, [hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
366
+   end
367
+
368
+   if ~isempty(img_idx)
369
+      hdr.dime.dim(5) = length(img_idx);
370
+   end
371
+
372
+   if ~isempty(dim5_idx)
373
+      hdr.dime.dim(6) = length(dim5_idx);
374
+   end
375
+
376
+   if ~isempty(dim6_idx)
377
+      hdr.dime.dim(7) = length(dim6_idx);
378
+   end
379
+
380
+   if ~isempty(dim7_idx)
381
+      hdr.dime.dim(8) = length(dim7_idx);
382
+   end
383
+
384
+   return						% read_image
385
+