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,320 @@
1
+%  Load NIFTI dataset header. Support both *.nii and *.hdr/*.img file
2
+%  extension. If file extension is not provided, *.hdr/*.img will be 
3
+%  used as default.
4
+%  
5
+%  Usage: [hdr, filetype, fileprefix, machine] = load_nii_hdr(filename)
6
+%  
7
+%  filename - NIFTI file name.
8
+%  
9
+%  Returned values:
10
+%  
11
+%  hdr - struct with NIFTI header fields.
12
+%  
13
+%  filetype	- 0 for Analyze format (*.hdr/*.img);
14
+%		  1 for NIFTI format in 2 files (*.hdr/*.img);
15
+%		  2 for NIFTI format in 1 file (*.nii).
16
+%  
17
+%  fileprefix - NIFTI file name without extension.
18
+%  
19
+%  machine    - a string, see below for details. The default here is 'ieee-le'.
20
+%
21
+%    'native'      or 'n' - local machine format - the default
22
+%    'ieee-le'     or 'l' - IEEE floating point with little-endian
23
+%                           byte ordering
24
+%    'ieee-be'     or 'b' - IEEE floating point with big-endian
25
+%                           byte ordering
26
+%    'vaxd'        or 'd' - VAX D floating point and VAX ordering
27
+%    'vaxg'        or 'g' - VAX G floating point and VAX ordering
28
+%    'cray'        or 'c' - Cray floating point with big-endian
29
+%                           byte ordering
30
+%    'ieee-le.l64' or 'a' - IEEE floating point with little-endian
31
+%                           byte ordering and 64 bit long data type
32
+%    'ieee-be.l64' or 's' - IEEE floating point with big-endian byte
33
+%                           ordering and 64 bit long data type.
34
+%  
35
+%  Number of scanned images in the file can be obtained by:
36
+%  num_scan = hdr.dime.dim(5)
37
+%
38
+%  Part of this file is copied and modified under GNU license from
39
+%  MRI_TOOLBOX developed by CNSP in Flinders University, Australia
40
+%
41
+%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
42
+%
43
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
44
+%
45
+function [hdr, filetype, fileprefix, machine] = load_nii_hdr(fileprefix)
46
+
47
+   if ~exist('fileprefix','var'),
48
+      error('Usage: [hdr, filetype, fileprefix, machine] = load_nii_hdr(filename)');
49
+   end
50
+
51
+   machine = 'ieee-le';
52
+   new_ext = 0;
53
+
54
+   if findstr('.nii',fileprefix)
55
+      new_ext = 1;
56
+      fileprefix = strrep(fileprefix,'.nii','');
57
+   end
58
+
59
+   if findstr('.hdr',fileprefix)
60
+      fileprefix = strrep(fileprefix,'.hdr','');
61
+   end
62
+
63
+   if findstr('.img',fileprefix)
64
+      fileprefix = strrep(fileprefix,'.img','');
65
+   end
66
+
67
+   if new_ext
68
+      fn = sprintf('%s.nii',fileprefix);
69
+
70
+      if ~exist(fn)
71
+         msg = sprintf('Cannot find file "%s.nii".', fileprefix);
72
+         error(msg);
73
+      end
74
+   else
75
+      fn = sprintf('%s.hdr',fileprefix);
76
+
77
+      if ~exist(fn)
78
+         msg = sprintf('Cannot find file "%s.hdr".', fileprefix);
79
+         error(msg);
80
+      end
81
+   end
82
+
83
+   fid = fopen(fn,'r',machine);
84
+    
85
+   if fid < 0,
86
+      msg = sprintf('Cannot open file %s.',fn);
87
+      error(msg);
88
+   else
89
+      fseek(fid,0,'bof');
90
+
91
+      if fread(fid,1,'int32') == 348
92
+         hdr = read_header(fid);
93
+         fclose(fid);
94
+      else
95
+         fclose(fid);
96
+
97
+         %  first try reading the opposite endian to 'machine'
98
+         %
99
+         switch machine,
100
+         case 'ieee-le', machine = 'ieee-be';
101
+         case 'ieee-be', machine = 'ieee-le';
102
+         end
103
+
104
+         fid = fopen(fn,'r',machine);
105
+
106
+         if fid < 0,
107
+            msg = sprintf('Cannot open file %s.',fn);
108
+            error(msg);
109
+         else
110
+            fseek(fid,0,'bof');
111
+
112
+            if fread(fid,1,'int32') ~= 348
113
+
114
+               %  Now throw an error
115
+               %
116
+               msg = sprintf('File "%s" is corrupted.',fn);
117
+               error(msg);
118
+            end
119
+
120
+            hdr = read_header(fid);
121
+            fclose(fid);
122
+         end
123
+      end
124
+   end
125
+
126
+   if strcmp(hdr.hist.magic, 'n+1')
127
+      filetype = 2;
128
+   elseif strcmp(hdr.hist.magic, 'ni1')
129
+      filetype = 1;
130
+   else
131
+      filetype = 0;
132
+   end
133
+
134
+   return					% load_nii_hdr
135
+
136
+
137
+%---------------------------------------------------------------------
138
+function [ dsr ] = read_header(fid)
139
+
140
+        %  Original header structures
141
+	%  struct dsr
142
+	%       { 
143
+	%       struct header_key hk;            /*   0 +  40       */
144
+	%       struct image_dimension dime;     /*  40 + 108       */
145
+	%       struct data_history hist;        /* 148 + 200       */
146
+	%       };                               /* total= 348 bytes*/
147
+
148
+    dsr.hk   = header_key(fid);
149
+    dsr.dime = image_dimension(fid);
150
+    dsr.hist = data_history(fid);
151
+
152
+    %  For Analyze data format
153
+    %
154
+    if ~strcmp(dsr.hist.magic, 'n+1') & ~strcmp(dsr.hist.magic, 'ni1')
155
+        dsr.hist.qform_code = 0;
156
+        dsr.hist.sform_code = 0;
157
+    end
158
+
159
+    return					% read_header
160
+
161
+
162
+%---------------------------------------------------------------------
163
+function [ hk ] = header_key(fid)
164
+
165
+    fseek(fid,0,'bof');
166
+    
167
+	%  Original header structures	
168
+	%  struct header_key                     /* header key      */ 
169
+	%       {                                /* off + size      */
170
+	%       int sizeof_hdr                   /*  0 +  4         */
171
+	%       char data_type[10];              /*  4 + 10         */
172
+	%       char db_name[18];                /* 14 + 18         */
173
+	%       int extents;                     /* 32 +  4         */
174
+	%       short int session_error;         /* 36 +  2         */
175
+	%       char regular;                    /* 38 +  1         */
176
+	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
177
+	%       };                               /* total=40 bytes  */
178
+	%
179
+	% int sizeof_header   Should be 348.
180
+	% char regular        Must be 'r' to indicate that all images and 
181
+	%                     volumes are the same size. 
182
+
183
+    v6 = version;
184
+    if str2num(v6(1))<6
185
+       directchar = '*char';
186
+    else
187
+       directchar = 'uchar=>char';
188
+    end
189
+
190
+    hk.sizeof_hdr    = fread(fid, 1,'int32')';	% should be 348!
191
+    hk.data_type     = deblank(fread(fid,10,directchar)');
192
+    hk.db_name       = deblank(fread(fid,18,directchar)');
193
+    hk.extents       = fread(fid, 1,'int32')';
194
+    hk.session_error = fread(fid, 1,'int16')';
195
+    hk.regular       = fread(fid, 1,directchar)';
196
+    hk.dim_info      = fread(fid, 1,'uchar')';
197
+    
198
+    return					% header_key
199
+
200
+
201
+%---------------------------------------------------------------------
202
+function [ dime ] = image_dimension(fid)
203
+
204
+	%  Original header structures    
205
+	%  struct image_dimension
206
+	%       {                                /* off + size      */
207
+	%       short int dim[8];                /* 0 + 16          */
208
+        %       /*
209
+        %           dim[0]      Number of dimensions in database; usually 4. 
210
+        %           dim[1]      Image X dimension;  number of *pixels* in an image row. 
211
+        %           dim[2]      Image Y dimension;  number of *pixel rows* in slice. 
212
+        %           dim[3]      Volume Z dimension; number of *slices* in a volume. 
213
+        %           dim[4]      Time points; number of volumes in database
214
+        %       */
215
+	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
216
+	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
217
+	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
218
+	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
219
+	%       short int datatype;              /* 30 + 2          */
220
+	%       short int bitpix;                /* 32 + 2          */
221
+	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
222
+	%       float pixdim[8];                 /* 36 + 32         */
223
+	%	/*
224
+	%		pixdim[] specifies the voxel dimensions:
225
+	%		pixdim[1] - voxel width, mm
226
+	%		pixdim[2] - voxel height, mm
227
+	%		pixdim[3] - slice thickness, mm
228
+	%		pixdim[4] - volume timing, in msec
229
+	%					..etc
230
+	%	*/
231
+	%       float vox_offset;                /* 68 + 4          */
232
+	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
233
+	%       float scl_inter;   % float funused1;      /* 76 + 4 */
234
+	%       short slice_end;   % float funused2;      /* 80 + 2 */
235
+	%       char slice_code;   % float funused2;      /* 82 + 1 */
236
+	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
237
+	%       float cal_max;                   /* 84 + 4          */
238
+	%       float cal_min;                   /* 88 + 4          */
239
+	%       float slice_duration;   % int compressed; /* 92 + 4 */
240
+	%       float toffset;   % int verified;          /* 96 + 4 */
241
+	%       int glmax;                       /* 100 + 4         */
242
+	%       int glmin;                       /* 104 + 4         */
243
+	%       };                               /* total=108 bytes */
244
+	
245
+    dime.dim        = fread(fid,8,'int16')';
246
+    dime.intent_p1  = fread(fid,1,'float32')';
247
+    dime.intent_p2  = fread(fid,1,'float32')';
248
+    dime.intent_p3  = fread(fid,1,'float32')';
249
+    dime.intent_code = fread(fid,1,'int16')';
250
+    dime.datatype   = fread(fid,1,'int16')';
251
+    dime.bitpix     = fread(fid,1,'int16')';
252
+    dime.slice_start = fread(fid,1,'int16')';
253
+    dime.pixdim     = abs(fread(fid,8,'float32')');
254
+    dime.vox_offset = fread(fid,1,'float32')';
255
+    dime.scl_slope  = fread(fid,1,'float32')';
256
+    dime.scl_inter  = fread(fid,1,'float32')';
257
+    dime.slice_end  = fread(fid,1,'int16')';
258
+    dime.slice_code = fread(fid,1,'uchar')';
259
+    dime.xyzt_units = fread(fid,1,'uchar')';
260
+    dime.cal_max    = fread(fid,1,'float32')';
261
+    dime.cal_min    = fread(fid,1,'float32')';
262
+    dime.slice_duration = fread(fid,1,'float32')';
263
+    dime.toffset    = fread(fid,1,'float32')';
264
+    dime.glmax      = fread(fid,1,'int32')';
265
+    dime.glmin      = fread(fid,1,'int32')';
266
+        
267
+    return					% image_dimension
268
+
269
+
270
+%---------------------------------------------------------------------
271
+function [ hist ] = data_history(fid)
272
+        
273
+	%  Original header structures
274
+	%  struct data_history       
275
+	%       {                                /* off + size      */
276
+	%       char descrip[80];                /* 0 + 80          */
277
+	%       char aux_file[24];               /* 80 + 24         */
278
+	%       short int qform_code;            /* 104 + 2         */
279
+	%       short int sform_code;            /* 106 + 2         */
280
+	%       float quatern_b;                 /* 108 + 4         */
281
+	%       float quatern_c;                 /* 112 + 4         */
282
+	%       float quatern_d;                 /* 116 + 4         */
283
+	%       float qoffset_x;                 /* 120 + 4         */
284
+	%       float qoffset_y;                 /* 124 + 4         */
285
+	%       float qoffset_z;                 /* 128 + 4         */
286
+	%       float srow_x[4];                 /* 132 + 16        */
287
+	%       float srow_y[4];                 /* 148 + 16        */
288
+	%       float srow_z[4];                 /* 164 + 16        */
289
+	%       char intent_name[16];            /* 180 + 16        */
290
+	%       char magic[4];   % int smin;     /* 196 + 4         */
291
+	%       };                               /* total=200 bytes */
292
+
293
+    v6 = version;
294
+    if str2num(v6(1))<6
295
+       directchar = '*char';
296
+    else
297
+       directchar = 'uchar=>char';
298
+    end
299
+
300
+    hist.descrip     = deblank(fread(fid,80,directchar)');
301
+    hist.aux_file    = deblank(fread(fid,24,directchar)');
302
+    hist.qform_code  = fread(fid,1,'int16')';
303
+    hist.sform_code  = fread(fid,1,'int16')';
304
+    hist.quatern_b   = fread(fid,1,'float32')';
305
+    hist.quatern_c   = fread(fid,1,'float32')';
306
+    hist.quatern_d   = fread(fid,1,'float32')';
307
+    hist.qoffset_x   = fread(fid,1,'float32')';
308
+    hist.qoffset_y   = fread(fid,1,'float32')';
309
+    hist.qoffset_z   = fread(fid,1,'float32')';
310
+    hist.srow_x      = fread(fid,4,'float32')';
311
+    hist.srow_y      = fread(fid,4,'float32')';
312
+    hist.srow_z      = fread(fid,4,'float32')';
313
+    hist.intent_name = deblank(fread(fid,16,directchar)');
314
+    hist.magic       = deblank(fread(fid,4,directchar)');
315
+
316
+    fseek(fid,253,'bof');
317
+    hist.originator  = fread(fid, 5,'int16')';
318
+    
319
+    return					% data_history
320
+