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,197 @@
1
+%  Save NIFTI or ANALYZE dataset that is loaded by "load_untouch_nii.m".
2
+%  The output image format and file extension will be the same as the
3
+%  input one (NIFTI.nii, NIFTI.img or ANALYZE.img). Therefore, any file
4
+%  extension that you specified will be ignored.
5
+%
6
+%  Usage: save_untouch_nii(nii, filename)
7
+%  
8
+%  nii - nii structure that is loaded by "load_untouch_nii.m"
9
+%
10
+%  filename  - 	NIFTI or ANALYZE file name.
11
+%
12
+%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
13
+%
14
+function save_untouch_nii(nii, filename)
15
+   
16
+   if ~exist('nii','var') | isempty(nii) | ~isfield(nii,'hdr') | ...
17
+	~isfield(nii,'img') | ~exist('filename','var') | isempty(filename)
18
+
19
+      error('Usage: save_untouch_nii(nii, filename)');
20
+   end
21
+
22
+   if ~isfield(nii,'untouch') | nii.untouch == 0
23
+      error('Usage: please use ''save_nii.m'' for the modified structure.');
24
+   end
25
+
26
+   if isfield(nii.hdr.hist,'magic') & strcmp(nii.hdr.hist.magic(1:3),'ni1')
27
+      filetype = 1;
28
+   elseif isfield(nii.hdr.hist,'magic') & strcmp(nii.hdr.hist.magic(1:3),'n+1')
29
+      filetype = 2;
30
+   else
31
+      filetype = 0;
32
+   end
33
+
34
+   [p,f] = fileparts(filename);
35
+   fileprefix = fullfile(p, f);
36
+
37
+   write_nii(nii, filetype, fileprefix);
38
+
39
+%   %  So earlier versions of SPM can also open it with correct originator
40
+ %  %
41
+  % if filetype == 0
42
+   %   M=[[diag(nii.hdr.dime.pixdim(2:4)) -[nii.hdr.hist.originator(1:3).*nii.hdr.dime.pixdim(2:4)]'];[0 0 0 1]];
43
+    %  save(fileprefix, 'M');
44
+%   elseif filetype == 1
45
+ %     M=[];
46
+  %    save(fileprefix, 'M');
47
+   %end
48
+   
49
+   return					% save_untouch_nii
50
+
51
+
52
+%-----------------------------------------------------------------------------------
53
+function write_nii(nii, filetype, fileprefix)
54
+
55
+   hdr = nii.hdr;
56
+
57
+   if isfield(nii,'ext') & ~isempty(nii.ext)
58
+      ext = nii.ext;
59
+      [ext, esize_total] = verify_nii_ext(ext);
60
+   else
61
+      ext = [];
62
+   end
63
+
64
+   switch double(hdr.dime.datatype),
65
+   case   1,
66
+      hdr.dime.bitpix = int16(1 ); precision = 'ubit1';
67
+   case   2,
68
+      hdr.dime.bitpix = int16(8 ); precision = 'uint8';
69
+   case   4,
70
+      hdr.dime.bitpix = int16(16); precision = 'int16';
71
+   case   8,
72
+      hdr.dime.bitpix = int16(32); precision = 'int32';
73
+   case  16,
74
+      hdr.dime.bitpix = int16(32); precision = 'float32';
75
+   case  32,
76
+      hdr.dime.bitpix = int16(64); precision = 'float32';
77
+   case  64,
78
+      hdr.dime.bitpix = int16(64); precision = 'float64';
79
+   case 128,
80
+      hdr.dime.bitpix = int16(24); precision = 'uint8';
81
+   case 256 
82
+      hdr.dime.bitpix = int16(8 ); precision = 'int8';
83
+   case 512 
84
+      hdr.dime.bitpix = int16(16); precision = 'uint16';
85
+   case 768 
86
+      hdr.dime.bitpix = int16(32); precision = 'uint32';
87
+   case 1024
88
+      hdr.dime.bitpix = int16(64); precision = 'int64';
89
+   case 1280
90
+      hdr.dime.bitpix = int16(64); precision = 'uint64';
91
+   case 1792,
92
+      hdr.dime.bitpix = int16(128); precision = 'float64';
93
+   otherwise
94
+      error('This datatype is not supported');
95
+   end
96
+   
97
+%   hdr.dime.glmax = round(double(max(nii.img(:))));
98
+ %  hdr.dime.glmin = round(double(min(nii.img(:))));
99
+   
100
+   if filetype == 2
101
+      fid = fopen(sprintf('%s.nii',fileprefix),'w');
102
+      
103
+      if fid < 0,
104
+         msg = sprintf('Cannot open file %s.nii.',fileprefix);
105
+         error(msg);
106
+      end
107
+      
108
+      hdr.dime.vox_offset = 352;
109
+
110
+      if ~isempty(ext)
111
+         hdr.dime.vox_offset = hdr.dime.vox_offset + esize_total;
112
+      end
113
+
114
+      hdr.hist.magic = 'n+1';
115
+      save_untouch_nii_hdr(hdr, fid);
116
+
117
+      if ~isempty(ext)
118
+         save_nii_ext(ext, fid);
119
+      end
120
+   elseif filetype == 1
121
+      fid = fopen(sprintf('%s.hdr',fileprefix),'w');
122
+      
123
+      if fid < 0,
124
+         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
125
+         error(msg);
126
+      end
127
+      
128
+      hdr.dime.vox_offset = 0;
129
+      hdr.hist.magic = 'ni1';
130
+      save_untouch_nii_hdr(hdr, fid);
131
+
132
+      if ~isempty(ext)
133
+         save_nii_ext(ext, fid);
134
+      end
135
+      
136
+      fclose(fid);
137
+      fid = fopen(sprintf('%s.img',fileprefix),'w');
138
+   else
139
+      fid = fopen(sprintf('%s.hdr',fileprefix),'w');
140
+      
141
+      if fid < 0,
142
+         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
143
+         error(msg);
144
+      end
145
+      
146
+      save_untouch0_nii_hdr(hdr, fid);
147
+      
148
+      fclose(fid);
149
+      fid = fopen(sprintf('%s.img',fileprefix),'w');
150
+   end
151
+
152
+   ScanDim = double(hdr.dime.dim(5));		% t
153
+   SliceDim = double(hdr.dime.dim(4));		% z
154
+   RowDim   = double(hdr.dime.dim(3));		% y
155
+   PixelDim = double(hdr.dime.dim(2));		% x
156
+   SliceSz  = double(hdr.dime.pixdim(4));
157
+   RowSz    = double(hdr.dime.pixdim(3));
158
+   PixelSz  = double(hdr.dime.pixdim(2));
159
+   
160
+   x = 1:PixelDim;
161
+   
162
+   if filetype == 2 & isempty(ext)
163
+      skip_bytes = double(hdr.dime.vox_offset) - 348;
164
+   else
165
+      skip_bytes = 0;
166
+   end
167
+
168
+   if double(hdr.dime.datatype) == 128
169
+
170
+      %  RGB planes are expected to be in the 4th dimension of nii.img
171
+      %
172
+      if(size(nii.img,4)~=3)
173
+         error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
174
+      end
175
+
176
+      nii.img = permute(nii.img, [4 1 2 3 5]);
177
+   end
178
+
179
+   %  For complex float32 or complex float64, voxel values
180
+   %  include [real, imag]
181
+   %
182
+   if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
183
+      real_img = real(nii.img(:))';
184
+      nii.img = imag(nii.img(:))';
185
+      nii.img = [real_img; nii.img];
186
+   end
187
+
188
+   if skip_bytes
189
+      fwrite(fid, ones(1,skip_bytes), 'uint8');
190
+   end
191
+
192
+   fwrite(fid, nii.img, precision);
193
+%   fwrite(fid, nii.img, precision, skip_bytes);        % error using skip
194
+   fclose(fid);
195
+
196
+   return;					% write_nii
197
+