source: aedes_write_nifti.m @ 151

Last change on this file since 151 was 101, checked in by tjniskan, 9 years ago
  • aedes_write_nifti now accepts logical data

M aedes_write_nifti.m
M aedes_revision.m

File size: 33.9 KB
Line 
1function [done,msg] = aedes_write_nifti(DATA,filename,varargin)
2% AEDES_WRITE_NIFTI - Write data in NIfTI or Analyze 7.5 formats
3%   
4%
5% Synopsis:
6%       [done,msg] = aedes_write_nifti(DATA,filename,prop1,value1,prop2,value2,...);
7%
8% Description:
9%       The function writes "DATA" to a file "filename" in NIfTI
10%       one file (default), NIfTI two file, or Analyze 7.5
11%       format. Function returns a done=1, if the writing was
12%       successful. Otherwise, done=0 and the second output argument msg
13%       contains the error message.
14%
15%       DATA can be a valid Aedes DATA-structure or a 2-D, 3-D, or 4-D
16%       matrix. Filename is the full path to the output file. If the
17%       filename is given without path, the data file is written into the
18%       working directory (pwd). If only one input argument is given or
19%       the filename is an empty string, the output file name is prompted.
20%
21%       The possible property/value pairs are the following:
22%
23%       Property        Value ({ }=default)       Desc.
24%       --------        --------                  --------
25%
26%       'FileType'      [ 0 | 1 | {2} ]           % 0 = Analyze 7.5
27%                                                 % 1 = NIfTI (two file)
28%                                                 % 2 = NIfTI (one file)
29%                                                 % (default)
30%
31%       'DataType'      [ {0}    | 'uint8' |      % 0 = Determine
32%                        'int16' | 'uint16'|      % datatype from data
33%                        'int32' | 'uint32'|      % (default)
34%                        'int64' | 'uint64'|      % 'str' = specify
35%                        'single'| 'double' ]     % datatype
36%
37%       'VoxelSize'     vox_size                  % A vector
38%                                                 % corresponding to the
39%                                                 % voxel width in
40%                                                 % dimension i:
41%                                                 % vox_size(1)= x width
42%                                                 % vox_size(2)= y width
43%                                                 % vox_size(3)= z width
44%                                                 % vox_size(4)= time width
45%                                                 
46%       'XYZUnits'      [ 'meter' | 'mm' |        % Units of the spatial
47%                        'micron' | {'Unknown'} ] % x, y, and z dimensions
48%                                               
49%       'TimeUnits'     [ 'sec' | 'msec' |        % Units of the temporal
50%                        'usec' | 'hz' |          % dimensions
51%                         'ppm' | 'rad/sec' |
52%                         {'Unknown'} ]
53%
54%       'machine'       [ {'ieee-le'}|'ieee-be' ] % little-endian
55%                                                 % (default) or
56%                                                 % big-endian Byte
57%                                                 % ordering
58%                       
59%       'HeaderOnly'    [ {0} | 1 ]               % 0 = write header and
60%                                                 % data (default)
61%                                                 % 1 = write only header
62%
63%       'Clim'          display min/max           % The displayed min and
64%                       (1x2 vector, [min,max])   % max intensity values.
65%                                                 % The NIfTI viewer may
66%                                                 % or may not use these
67%                                                 % values. Clim=[0 0] is
68%                                                 % the default and means
69%                                                 % that the whole range of
70%                                                 % the data is shown.
71%
72%       'description'   (string)                  % A string to be written
73%                                                 % in the description
74%                                                 % field in NIfTI header.
75%
76%
77% Examples:
78%       DATA = aedes_readfid('testi.fid');
79%
80%       % Write one file NIfTI
81%       aedes_write_nifti(DATA,'testi.nii');
82%
83%       % Write two file NIfTI
84%       aedes_write_nifti(DATA,'testi.img','FileType',1)
85%
86%       % Write Analyze 7.5 *.hdr and *.img files
87%       aedes_write_nifti(DATA,'testi.img','FileType',0)
88%       
89% See also:
90%       AEDES, AEDES_READ_NIFTI
91
92%       Acknowledgment: This function is modified under GNU license from
93%       MRI_TOOLBOX developed by CNSP in Flinders University,
94%       Australia. Some parts are also originally written by Jimmy Shen.
95%       See the following links:
96%       http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8797&objectType=file
97%       http://eeg.sourceforge.net/
98
99%       NIfTI data format specifications can be found here:
100%       http://nifti.nimh.nih.gov
101
102% This function is a part of Aedes - A graphical tool for analyzing
103% medical images
104%
105% Copyright (C) 2006 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uku.fi>
106% Copyright (C) 2006 Jimmy Shen
107%
108% Department of Physics, Department of Neurobiology
109% University of Kuopio, FINLAND
110%
111% This program may be used under the terms of the GNU General Public
112% License version 2.0 as published by the Free Software Foundation
113% and appearing in the file LICENSE.TXT included in the packaging of
114% this program.
115%
116% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
117% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
118
119
120%% Default values
121Dat.filetype=2;
122Dat.machine='ieee-le';
123Dat.datatype=0;
124Dat.xyzunits = '';
125Dat.timeunits = '';
126Dat.voxelsize = [];
127Dat.Clim=[0 0];
128Dat.Descrip = '';
129Dat.FlipOrient = 0;
130done=false;
131writeOnlyHeader=0;
132
133
134%% Parse input arguments
135switch nargin
136 case 0 % ---------------------------------
137   done=0;
138   msg=['Too few input arguments.'];
139   return
140 
141 case {1, 2} % --------------------------------
142   if isnumeric(DATA) || islogical(DATA)
143     data = DATA;
144     if islogical(data)
145       data=uint8(data);
146     end
147     DATA=[];
148     DATA.FTDATA = data;
149     DATA.DataFormat = '';
150   elseif isstruct(DATA)
151     if ~isfield(DATA,'FTDATA')
152       error('The field FTDATA not found from the inputted structure!')
153     end
154     if ~isfield(DATA,'DataFormat')
155       DATA.DataFormat = '';
156     end
157   else
158     done=0;
159     msg=['First input argument must be structure or numeric type!'];
160     return
161   end
162   
163   %% Prompt for file name
164   if ~exist('filename','var') || isempty(filename)
165     [fname,fpath,findex]=uiputfile({'*.nii;*.NII',...
166                         'NIfTI Files - One File Format (*.nii)';...
167                   '*.hdr;*.HDR','NIfTI Files - Two File Format (*.hdr)';...
168                   '*.hdr;*.HDR','Analyze 7.5 Files (*.hdr)';...
169                   '*.*','All Files (*.*)'},'Save Data As');
170     if isequal(fname,0) || isequal(fpath,0)
171       % Cancelled
172       done = false;
173       msg = 'Cancelled';
174       return
175     else
176       if findex==1
177         Dat.filetype=2;
178       elseif findex==2
179         Dat.filetype=1;
180       elseif findex==3
181         Dat.filetype=0;
182       else
183         Dat.filetype=2;
184       end
185       filename = [fpath,fname];
186     end
187   end
188  otherwise
189   
190    if isnumeric(DATA)
191      data = DATA;
192      DATA=[];
193      DATA.FTDATA = data;
194      DATA.DataFormat = '';
195    elseif isstruct(DATA)
196      if ~isfield(DATA,'FTDATA')
197        error('The field FTDATA not found from the inputted structure!')
198      end
199      if ~isfield(DATA,'DataFormat')
200        DATA.DataFormat = '';
201      end
202    else
203      done=0;
204      msg=['First input argument must be structure or numeric type!'];
205      return
206    end
207   
208  %% Parse parameter/value pairs
209  for ii=1:2:length(varargin)
210    switch lower(varargin{ii})
211      case 'machine'
212        Dat.machine=varargin{ii+1};
213       
214      case 'datatype'
215        Dat.datatype=varargin{ii+1};
216       
217      case 'filetype'
218        Dat.filetype=varargin{ii+1};
219       
220      case 'headeronly'
221        writeOnlyHeader = varargin{ii+1};
222       
223      case 'xyzunits'
224        Dat.xyzunits = varargin{ii+1};
225       
226      case 'timeunits'
227        Dat.timeunits = varargin{ii+1};
228       
229      case 'voxelsize'
230        Dat.voxelsize = varargin{ii+1};
231       
232      case 'clim'
233        Dat.Clim = varargin{ii+1};
234               
235          case 'description'
236                Dat.Descrip = varargin{ii+1};
237               
238          case 'FlipInplaneOrient'
239                Dat.FlipOrient = varargin{ii+1};
240       
241     otherwise
242      msg=['Unknown parameter "' varargin{ii} '".'];
243      return
244    end
245  end
246end
247
248%% Get default datatype from DATA structure if not specified as an input
249if ischar(Dat.datatype)
250  % Datatype given as a string 'uint16', 'int32', 'single',...
251  Dat.datatype=l_GetDataType(Dat.datatype);
252elseif Dat.datatype==0
253  Dat.datatype=l_GetDataType(DATA);
254end
255
256%% Parse filename
257[fp,fn,fe]=fileparts(filename);
258if ~isempty(fp)
259  fp = [fp,filesep];
260else
261  fp = [pwd,filesep];
262end
263 
264if ~any(strcmpi(fe,{'.nii','.hdr','.img'}))
265  fn=[fn,fe];
266end
267
268
269%% Construct header structure
270[done,msg,hdr] = l_ConstructValidHeader(DATA,Dat);
271if ~done
272  return
273end
274
275%% Open file(s) for writing
276if any(Dat.filetype==[0 1]) % Analyze75 and two file NIfTI
277  fid_hdr = fopen([fp,fn,'.hdr'],'w',Dat.machine);
278  if fid_hdr<0
279        done = false;
280    msg={'Could not open file',['"',fp,fn,'.hdr"'],...
281         'for writing'};
282    return
283  end
284  fid_img = fopen([fp,fn,'.img'],'w',Dat.machine);
285  if fid_hdr<0
286        done = false;
287    msg={'Could not open file',['"',fp,fn,'.img"'],...
288         'for writing'};
289    return
290  end
291 
292  %% Write header
293  [done,msg] = l_WriteNiftiHdr(hdr,fid_hdr,Dat);
294  if ~done
295    fclose(fid_hdr);
296    return
297  end
298  fclose(fid_hdr);
299 
300  %% Write data
301  [done,msg] = l_WriteNiftiData(DATA,hdr,fid_img,Dat);
302  if ~done
303    fclose(fid_img);
304    return
305  end
306  fclose(fid_img);
307elseif Dat.filetype==2 % On file NIfTI (default)
308  fid = fopen([fp,fn,'.nii'],'w',Dat.machine);
309  if fid<0
310        done = false;
311    msg={'Could not open file',['"',fp,fn,'.nii"'],...
312         'for writing'};
313    return
314  end
315 
316  %% Write header
317  [done,msg] = l_WriteNiftiHdr(hdr,fid,Dat);
318  if ~done
319    fclose(fid);
320    return
321  end
322 
323  %% Write data
324  [done,msg] = l_WriteNiftiData(DATA,hdr,fid,Dat);
325  if ~done
326    fclose(fid);
327    return
328  end
329  fclose(fid);
330else
331  msg={'Unsupported filetype.',...
332       '0=Analyze75 format (*.hdr,*.img)',...
333       '1=NIfTI (*.hdr, *.img)',...
334       '2=NIfTI (*.nii)'};
335end
336
337
338%% All seems to be well and we can exit normally without errors
339done=true;
340return
341
342
343
344%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
345% Construct header structure
346%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
347function [done,msg,hdr] = l_ConstructValidHeader(DATA,Dat)
348
349hdr = [];
350msg='';
351done=false;
352
353%% If Analyze75 <-> NIfTI
354try
355  if isfield(DATA,'DataFormat') && any(strcmpi(DATA.DataFormat,{'Analyze75','NIfTI(1)','NIfTI(2)'}))
356    hdr = DATA.HDR.FileHeader;
357   
358    %% NIfTI -> Analyze 7.5
359    if Dat.filetype==0 && any(strcmpi(DATA.DataFormat,{'NIfTI(1)','NIfTI(2)'}))
360         
361      hdr.hist.orient      = char(0);
362      hdr.hist.originator=zeros(1,5);
363      hdr.hist.originator(1)=round(hdr.dime.dim(2)/2);
364      hdr.hist.originator(2)=round(hdr.dime.dim(3)/2);
365      hdr.hist.originator(3)=round(hdr.dime.dim(4)/2);
366      hdr.hist.generated   = ['Aedes',char([0 0 0 0 0])];
367      hdr.hist.scannum     = char(zeros(1,10));
368      hdr.hist.patient_id  = char(zeros(1,10));
369      hdr.hist.exp_date    = char(zeros(1,10));
370      hdr.hist.exp_time    = char(zeros(1,10));
371      hdr.hist.hist_un0    = char(zeros(1,3));
372      hdr.hist.views       = 0;
373      hdr.hist.vols_added  = 0;
374      hdr.hist.start_field = 0;
375      hdr.hist.field_skip  = 0;
376      hdr.hist.omax        = 0;
377      hdr.hist.omin        = 0;
378      hdr.hist.smax        = 0;
379      hdr.hist.smin        = 0;
380      hdr.hist.magic = '';
381      hdr.dime.vox_offset=0;
382          hdr.dime.pixdim(1)=0;
383      if any(Dat.Clim~=0)
384        hdr.dime.cal_max = Dat.Clim(2);
385        hdr.dime.cal_min = Dat.Clim(1);
386      end
387     
388    %% Analyze 7.5 -> NIfTI
389    elseif any(Dat.filetype==[1 2]) && strcmpi(DATA.DataFormat,'Analyze75')
390      hdr.hist.qform_code=0;
391      hdr.hist.sform_code=0;
392      hdr.hist.quatern_b=0;
393      hdr.hist.quatern_c=0;
394      hdr.hist.quatern_d=0;
395      hdr.hist.qoffset_x=0;
396      hdr.hist.qoffset_y=0;
397      hdr.hist.qoffset_z=0;
398      hdr.hist.srow_x=zeros(1,4);
399      hdr.hist.srow_y=zeros(1,4);
400      hdr.hist.srow_z=zeros(1,4);
401      hdr.hist.intent_name='';
402      hdr.hist.originator=zeros(1,5);
403      hdr.hist.originator(1)=round(hdr.dime.dim(2)/2);
404      hdr.hist.originator(2)=round(hdr.dime.dim(3)/2);
405      hdr.hist.originator(3)=round(hdr.dime.dim(4)/2);
406      if Dat.filetype==1
407        hdr.hist.magic='ni1';
408        hdr.dime.vox_offset=0;
409      elseif Dat.filetype==2
410        hdr.hist.magic='n+1';
411        hdr.dime.vox_offset=352;
412        hdr.hist.sform_code = 1;
413        hdr.hist.srow_x(1) = hdr.dime.pixdim(2);
414        hdr.hist.srow_y(2) = hdr.dime.pixdim(3);
415        hdr.hist.srow_z(3) = hdr.dime.pixdim(4);
416        hdr.hist.srow_x(4) = (1-hdr.hist.originator(1))*hdr.dime.pixdim(2);
417        hdr.hist.srow_y(4) = (1-hdr.hist.originator(2))*hdr.dime.pixdim(3);
418        hdr.hist.srow_z(4) = (1-hdr.hist.originator(3))*hdr.dime.pixdim(4);
419      end
420      if any(Dat.Clim~=0)
421        hdr.dime.cal_max = Dat.Clim(2);
422        hdr.dime.cal_min = Dat.Clim(1);
423      end
424     
425    %% NIfTI -> NIfTI or Analyze75 -> Analyze75
426    else
427      if any(Dat.Clim~=0)
428        hdr.dime.cal_max = Dat.Clim(2);
429        hdr.dime.cal_min = Dat.Clim(1);
430          end
431          if ~isempty(Dat.voxelsize)
432                hdr.dime.pixdim(2:length(Dat.voxelsize)+1) = Dat.voxelsize;
433          end
434          if ~isempty(Dat.Descrip)
435                hdr.hist.descrip=Dat.Descrip;
436          end
437          data_dim(1)=size(DATA.FTDATA,1);
438          data_dim(2)=size(DATA.FTDATA,2);
439          data_dim(3)=size(DATA.FTDATA,3);
440          data_dim(4)=size(DATA.FTDATA,4);
441         
442          % Swap 1st and 2nd dimesions
443          data_dim = [data_dim(2) data_dim(1) data_dim(3:end)];
444         
445          hdr.dime.dim=[length(size(DATA.FTDATA)) data_dim 1 1 1]; % Data dimensions:
446                                           % dim[0] = nbr of dimensions
447                                           % dim[i] = length of dimension i
448                                           % (i=1..7)
449         
450          % Make sure that datatype is valid
451          % Possible value pairs for DataType and BitPix
452          d_type = [0 1 2 4 8 16 32 64 128 256 511 512 768 1024 1280 1536 1792 2048];
453          b_pix = [0 1 8 16 32 32 64 64 24 8 96 16 32 64 64 128 128 256];
454          hdr.dime.datatype = Dat.datatype;
455          hdr.dime.bitpix = b_pix(d_type==Dat.datatype);
456         
457        end
458% $$$       if hdr.hist.qform_code == 0 & hdr.hist.sform_code == 0
459% $$$         hdr.hist.sform_code = 1;
460% $$$         hdr.hist.srow_x(1) = hdr.dime.pixdim(2);
461% $$$         hdr.hist.srow_y(2) = hdr.dime.pixdim(3);
462% $$$         hdr.hist.srow_z(3) = hdr.dime.pixdim(4);
463% $$$         hdr.hist.srow_x(4) = (1-hdr.hist.originator(1))*hdr.dime.pixdim(2);
464% $$$         hdr.hist.srow_y(4) = (1-hdr.hist.originator(2))*hdr.dime.pixdim(3);
465% $$$         hdr.hist.srow_z(4) = (1-hdr.hist.originator(3))*hdr.dime.pixdim(4);
466% $$$       end
467% $$$       
468% $$$       %% Set the magic field
469% $$$       switch filetype
470% $$$        case 1 % NIfTI one file
471% $$$         hdr.hist.magic = 'ni1';
472% $$$         hdr.dime.vox_offset=0;
473% $$$        case 2 % NIfTI two file
474% $$$         hdr.hist.magic = 'n+1';
475% $$$         hdr.dime.vox_offset=352;
476% $$$        otherwise
477% $$$         msg=['Unexpected filetype encountered while constructing header.'];
478% $$$         return
479% $$$       end
480% $$$     end
481  else
482    %% For arbitrary DATA structure formats the header structure has to be
483    %  constructed from scratch...
484   
485    % Possible value pairs for DataType and BitPix
486    d_type = [0 1 2 4 8 16 32 64 128 256 511 512 768 1024 1280 1536 1792 2048];
487    b_pix = [0 1 8 16 32 32 64 64 24 8 96 16 32 64 64 128 128 256];
488   
489    % Header key
490    hdr.hk.sizeof_hdr=348;
491    hdr.hk.data_type='';
492    hdr.hk.db_name='';
493    hdr.hk.extents=0;
494    hdr.hk.session_error=0;
495    hdr.hk.regular='r';
496    hdr.hk.dim_info=0;
497   
498    % Image dimensions
499    data_dim(1)=size(DATA.FTDATA,1);
500        data_dim(2)=size(DATA.FTDATA,2);
501        data_dim(3)=size(DATA.FTDATA,3);
502        data_dim(4)=size(DATA.FTDATA,4);
503       
504        % Swap 1st and 2nd dimesions
505        data_dim = [data_dim(2) data_dim(1) data_dim(3:end)];
506       
507    hdr.dime.dim=[length(size(DATA.FTDATA)) data_dim 1 1 1]; % Data dimensions:
508                                     % dim[0] = nbr of dimensions
509                                     % dim[i] = length of dimension i (i=1..7)
510    hdr.dime.intent_p1=0;
511    hdr.dime.intent_p2=0;
512    hdr.dime.intent_p3=0;
513    hdr.dime.intent_code=0;
514    hdr.dime.datatype = Dat.datatype;
515    hdr.dime.bitpix = b_pix(d_type==Dat.datatype);
516    hdr.dime.slice_start = 0;
517    hdr.dime.pixdim = zeros(1,8); % pixdim[i] = voxel width along
518                                  % dimension i
519    if ~isempty(Dat.voxelsize)
520      %hdr.dime.pixdim(1)=length(Dat.voxelsize);
521      hdr.dime.pixdim(2:length(Dat.voxelsize)+1) = Dat.voxelsize;
522    else
523      hdr.dime.pixdim(2:4)=1;
524    end
525    hdr.dime.vox_offset=0;
526    hdr.dime.scl_slope=0;
527    hdr.dime.scl_inter=0;
528    hdr.dime.slice_end=0;
529    hdr.dime.slice_code=0;
530   
531
532        xyz_bits = '000'; % Unknown
533        time_bits = '000'; % Unknown
534    hdr.dime.xyzt_units=bin2dec(['00',time_bits,xyz_bits]);
535   
536    hdr.dime.cal_max=Dat.Clim(2);
537    hdr.dime.cal_min=Dat.Clim(1);
538    hdr.dime.slice_duration=0;
539    hdr.dime.toffset=0;
540    hdr.dime.glmax=round(max(DATA.FTDATA(:)));
541    hdr.dime.glmin=round(min(DATA.FTDATA(:)));
542   
543    % Data history
544        if ~isempty(Dat.Descrip)
545          hdr.hist.descrip=Dat.Descrip;
546        else
547          hdr.hist.descrip=['Converted from ' DATA.DataFormat];
548        end
549    hdr.hist.aux_file='none';
550   
551   
552    %% Analyze 7.5 style data history
553    if Dat.filetype==0
554      hdr.hist.orient      = char(0);
555      hdr.hist.originator=zeros(1,5);
556      hdr.hist.originator(1)=round(hdr.dime.dim(2)/2);
557      hdr.hist.originator(2)=round(hdr.dime.dim(3)/2);
558      hdr.hist.originator(3)=round(hdr.dime.dim(4)/2);
559      hdr.hist.generated   = ['Aedes',char([0 0])];
560      hdr.hist.scannum     = char(zeros(1,10));
561      hdr.hist.patient_id  = char(zeros(1,10));
562      hdr.hist.exp_date    = char(zeros(1,10));
563      hdr.hist.exp_time    = char(zeros(1,10));
564      hdr.hist.hist_un0    = char(zeros(1,3));
565      hdr.hist.views       = 0;
566      hdr.hist.vols_added  = 0;
567      hdr.hist.start_field = 0;
568      hdr.hist.field_skip  = 0;
569      hdr.hist.omax        = 0;
570      hdr.hist.omin        = 0;
571      hdr.hist.smax        = 0;
572      hdr.hist.smin        = 0;
573     
574      % Set the magic field (NIfTI identifier)
575      hdr.hist.magic       = '';
576    elseif any(Dat.filetype==[1 2])
577      hdr.hist.qform_code=0;
578      hdr.hist.sform_code=0;
579      hdr.hist.quatern_b=0;
580      hdr.hist.quatern_c=0;
581      hdr.hist.quatern_d=0;
582      hdr.hist.qoffset_x=0;
583      hdr.hist.qoffset_y=0;
584      hdr.hist.qoffset_z=0;
585      hdr.hist.srow_x=zeros(1,4);
586      hdr.hist.srow_y=zeros(1,4);
587      hdr.hist.srow_z=zeros(1,4);
588      hdr.hist.intent_name='';
589      hdr.hist.originator=zeros(1,5);
590      hdr.hist.originator(1)=round(hdr.dime.dim(2)/2);
591      hdr.hist.originator(2)=round(hdr.dime.dim(3)/2);
592      hdr.hist.originator(3)=round(hdr.dime.dim(4)/2);
593      if Dat.filetype==1
594        hdr.hist.magic='ni1';
595      elseif Dat.filetype==2
596        hdr.hist.magic='n+1';
597        hdr.dime.vox_offset=352;
598        hdr.hist.sform_code = 1;
599        hdr.hist.srow_x(1) = hdr.dime.pixdim(2);
600        hdr.hist.srow_y(2) = hdr.dime.pixdim(3);
601        hdr.hist.srow_z(3) = hdr.dime.pixdim(4);
602        hdr.hist.srow_x(4) = (1-hdr.hist.originator(1))*hdr.dime.pixdim(2);
603        hdr.hist.srow_y(4) = (1-hdr.hist.originator(2))*hdr.dime.pixdim(3);
604        hdr.hist.srow_z(4) = (1-hdr.hist.originator(3))*hdr.dime.pixdim(4);
605      end
606    end
607  end
608 
609  % Units of spatial and temporal dimensions:
610  % 0 = Unknown, 1 = meter, 2 = millimeter, 3 = micrometer, 8 = seconds
611  % 16 = milliseconds, 24 = microseconds, 32 = Hertz, 40 = ppm,
612  % 48 = rad/sec.
613  % Bits 0..2 of xyzt_units specify the units of pixdim[1..3]. Bits
614  % 3..5 of xyzt_units specify the units of pixdim[4] and are multiples
615  % of 8.
616  if not(isempty(Dat.xyzunits))
617        if not(isfield(DATA,'DataFormat') && strcmpi(DATA.DataFormat,'Analyze75'))
618          if isempty(Dat.xyzunits) || strcmpi(Dat.xyzunits,'unknown')
619                xyz_bits = '000';
620          elseif strcmpi(Dat.xyzunits,'meter')
621                xyz_bits = dec2bin(1,3);
622          elseif strcmpi(Dat.xyzunits,'mm')
623                xyz_bits = dec2bin(2,3);
624          elseif strcmpi(Dat.xyzunits,'micron')
625                xyz_bits = dec2bin(3,3);
626          end
627          if isempty(Dat.timeunits) || strcmpi(Dat.timeunits,'unknown')
628                time_bits = '000';
629          elseif strcmpi(Dat.timeunits,'sec')
630                time_bits = dec2bin(1,3);
631          elseif strcmpi(Dat.timeunits,'msec')
632                time_bits = dec2bin(2,3);
633          elseif strcmpi(Dat.timeunits,'usec')
634                time_bits = dec2bin(3,3);
635          elseif strcmpi(Dat.timeunits,'hz')
636                time_bits = dec2bin(4,3);
637          elseif strcmpi(Dat.timeunits,'ppm')
638                time_bits = dec2bin(5,3);
639          elseif strcmpi(Dat.timeunits,'rad/sec')
640                time_bits = dec2bin(6,3);
641          end
642          hdr.dime.xyzt_units=bin2dec(['00',time_bits,xyz_bits]);
643        end
644  end
645  done=true;
646catch
647  msg=lasterr;
648end
649
650
651%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
652% Define datatype
653%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
654function [datatype,bitpix]=l_GetDataType(DATA)
655
656%  Set bitpix according to datatype
657%
658%  /*Acceptable values for datatype are*/
659%
660%     0 None                   (Unknown bit per voxel)   % DT_NONE, DT_UNKNOWN
661%     1 Binary                       (ubit1, bitpix=1)   % DT_BINARY
662%     2 Unsigned char       (uchar or uint8, bitpix=8)   % DT_UINT8, NIFTI_TYPE_UINT8
663%     4 Signed short                 (int16, bitpix=16)  % DT_INT16, NIFTI_TYPE_INT16
664%     8 Signed integer               (int32, bitpix=32)  % DT_INT32, NIFTI_TYPE_INT32
665%    16 Floating point   (single or float32, bitpix=32)  % DT_FLOAT32, NIFTI_TYPE_FLOAT32
666%    32 Complex, 2 float32     (Unsupported, bitpix=64)  % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
667%    64 Double precision (double or float64, bitpix=64)  % DT_FLOAT64, NIFTI_TYPE_FLOAT64
668%   128 uint8 RGB                (Use uint8, bitpix=24)  % DT_RGB24, NIFTI_TYPE_RGB24
669%   256 Signed char          (schar or int8, bitpix=8)   % DT_INT8, NIFTI_TYPE_INT8
670%   511 Single RGB             (Use float32, bitpix=96)  % DT_RGB96, NIFTI_TYPE_RGB96
671%   512 Unsigned short              (uint16, bitpix=16)  % DT_UNINT16, NIFTI_TYPE_UNINT16
672%   768 Unsigned integer            (uint32, bitpix=32)  % DT_UNINT32, NIFTI_TYPE_UNINT32
673%  1024 Signed long long             (int64, bitpix=64)  % DT_INT64, NIFTI_TYPE_INT64
674%  1280 Unsigned long long          (uint64, bitpix=64)  % DT_UINT64, NIFTI_TYPE_UINT64
675%  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
676%  1792 Complex128, 2 float64  (Unsupported, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX12
677%  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128,
678%  NIFTI_TYPE_COMPLEX128
679if ischar(DATA)
680  action = DATA;
681else
682  action = class(DATA.FTDATA);
683end
684
685switch action
686 case {'uint8','logical'}
687  datatype = 2;
688  bitpix = 8;
689 case 'int16'
690  datatype = 4;
691  bitpix = 16;
692 case 'uint16'
693  datatype = 512;
694  bitpix = 16;
695 case 'int32'
696  datatype = 8;
697  bitpix = 32;
698 case 'uint32'
699  datatype = 768;
700  bitpix = 32;
701 case 'int64'
702  datatype = 1024;
703  bitpix = 64;
704 case 'uint64'
705  datatype = 1280;
706  bitpix = 64;
707 case {'single','float'}
708  datatype = 16;
709  bitpix = 32;
710 case 'double'
711  datatype = 64;
712  bitpix = 64;
713 otherwise
714  % Default to float
715  datatype = 16;
716  bitpix = 32;
717end
718
719
720%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
721% Write NIfTI/Analyze75 header
722%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
723function [done,msg] = l_WriteNiftiHdr(hdr,fid,Dat)
724
725done=false;
726msg='';
727
728%  Original header structures
729%  struct dsr                           /* dsr = hdr */
730%       {
731%       struct header_key hk;            /*   0 +  40       */
732%       struct image_dimension dime;     /*  40 + 108       */
733%       struct data_history hist;        /* 148 + 200       */
734%       };                               /* total= 348 bytes*/
735try
736  header_key(fid, hdr.hk);
737  image_dimension(fid, hdr.dime);
738  data_history(fid, hdr.hist,Dat.filetype);
739catch
740  msg=lasterr;
741  return
742end
743
744%  check the file size is 348 bytes
745%
746fbytes = ftell(fid);
747 
748if ~isequal(fbytes,348),
749  msg = sprintf('Header size is not 348 bytes.');
750  return
751end
752
753done=true;
754return;                                 % write_header
755
756
757%---------------------------------------------------------------------
758function header_key(fid, hk)
759
760fseek(fid,0,'bof');
761
762%  Original header structures   
763%  struct header_key                      /* header key      */
764%       {                                /* off + size      */
765%       int sizeof_hdr                   /*  0 +  4         */
766%       char data_type[10];              /*  4 + 10         */
767%       char db_name[18];                /* 14 + 18         */
768%       int extents;                     /* 32 +  4         */
769%       short int session_error;         /* 36 +  2         */
770%       char regular;                    /* 38 +  1         */
771%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
772%       };                               /* total=40 bytes  */
773
774fwrite(fid, hk.sizeof_hdr(1),    'int32');      % must be 348.
775
776% data_type = sprintf('%-10s',hk.data_type);    % ensure it is 10 chars from left
777% fwrite(fid, data_type(1:10), 'uchar');
778pad = zeros(1, 10-length(hk.data_type));
779hk.data_type = [hk.data_type  char(pad)];
780fwrite(fid, hk.data_type(1:10), 'uchar');
781
782% db_name   = sprintf('%-18s', hk.db_name);     % ensure it is 18 chars from left
783% fwrite(fid, db_name(1:18), 'uchar');
784pad = zeros(1, 18-length(hk.db_name));
785hk.db_name = [hk.db_name  char(pad)];
786fwrite(fid, hk.db_name(1:18), 'uchar');
787
788fwrite(fid, hk.extents(1),       'int32');
789fwrite(fid, hk.session_error(1), 'int16');
790fwrite(fid, hk.regular(1),       'uchar');      % might be uint8
791
792% fwrite(fid, hk.hkey_un0(1),    'uchar');
793% fwrite(fid, hk.hkey_un0(1),    'uint8');
794fwrite(fid, hk.dim_info(1),      'uchar');
795return;                                 % header_key
796
797
798%---------------------------------------------------------------------
799function image_dimension(fid, dime)
800
801%  Original header structures       
802%  struct image_dimension
803%       {                                /* off + size      */
804%       short int dim[8];                /* 0 + 16          */
805%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
806%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
807%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
808%       short int intent_code;   % short int unused1;   /* 28 + 2 */
809%       short int datatype;              /* 30 + 2          */
810%       short int bitpix;                /* 32 + 2          */
811%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
812%       float pixdim[8];                 /* 36 + 32         */
813%                       /*
814%                               pixdim[] specifies the voxel dimensions:
815%                               pixdim[1] - voxel width
816%                               pixdim[2] - voxel height
817%                               pixdim[3] - interslice distance
818%                               pixdim[4] - volume timing, in msec
819%                                       ..etc
820%                       */
821%       float vox_offset;                /* 68 + 4          */
822%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
823%       float scl_inter;   % float funused1;      /* 76 + 4 */
824%       short slice_end;   % float funused2;      /* 80 + 2 */
825%       char slice_code;   % float funused2;      /* 82 + 1 */
826%       char xyzt_units;   % float funused2;      /* 83 + 1 */
827%       float cal_max;                   /* 84 + 4          */
828%       float cal_min;                   /* 88 + 4          */
829%       float slice_duration;   % int compressed; /* 92 + 4 */
830%       float toffset;   % int verified;          /* 96 + 4 */
831%       int glmax;                       /* 100 + 4         */
832%       int glmin;                       /* 104 + 4         */
833%       };                               /* total=108 bytes */
834
835fwrite(fid, dime.dim(1:8),        'int16');
836fwrite(fid, dime.intent_p1(1),  'float32');
837fwrite(fid, dime.intent_p2(1),  'float32');
838fwrite(fid, dime.intent_p3(1),  'float32');
839fwrite(fid, dime.intent_code(1),  'int16');
840fwrite(fid, dime.datatype(1),     'int16');
841fwrite(fid, dime.bitpix(1),       'int16');
842fwrite(fid, dime.slice_start(1),  'int16');
843fwrite(fid, dime.pixdim(1:8),   'float32');
844fwrite(fid, dime.vox_offset(1), 'float32');
845fwrite(fid, dime.scl_slope(1),  'float32');
846fwrite(fid, dime.scl_inter(1),  'float32');
847fwrite(fid, dime.slice_end(1),    'int16');
848fwrite(fid, dime.slice_code(1),   'uchar');
849fwrite(fid, dime.xyzt_units(1),   'uchar');
850fwrite(fid, dime.cal_max(1),    'float32');
851fwrite(fid, dime.cal_min(1),    'float32');
852fwrite(fid, dime.slice_duration(1), 'float32');
853fwrite(fid, dime.toffset(1),    'float32');
854fwrite(fid, dime.glmax(1),        'int32');
855fwrite(fid, dime.glmin(1),        'int32');
856return;                                 % image_dimension
857
858
859%---------------------------------------------------------------------
860function data_history(fid, hist, filetype)
861
862% Original header structures
863%struct data_history       
864%       {                                /* off + size      */
865%       char descrip[80];                /* 0 + 80          */
866%       char aux_file[24];               /* 80 + 24         */
867%       short int qform_code;            /* 104 + 2         */
868%       short int sform_code;            /* 106 + 2         */
869%       float quatern_b;                 /* 108 + 4         */
870%       float quatern_c;                 /* 112 + 4         */
871%       float quatern_d;                 /* 116 + 4         */
872%       float qoffset_x;                 /* 120 + 4         */
873%       float qoffset_y;                 /* 124 + 4         */
874%       float qoffset_z;                 /* 128 + 4         */
875%       float srow_x[4];                 /* 132 + 16        */
876%       float srow_y[4];                 /* 148 + 16        */
877%       float srow_z[4];                 /* 164 + 16        */
878%       char intent_name[16];            /* 180 + 16        */
879%       char magic[4];   % int smin;     /* 196 + 4         */
880%       };                               /* total=200 bytes */
881
882% descrip     = sprintf('%-80s', hist.descrip);     % 80 chars from left
883% fwrite(fid, descrip(1:80),    'uchar');
884pad = zeros(1, 80-length(hist.descrip));
885hist.descrip = [hist.descrip  char(pad)];
886fwrite(fid, hist.descrip(1:80), 'uchar');
887
888% aux_file    = sprintf('%-24s', hist.aux_file);    % 24 chars from left
889% fwrite(fid, aux_file(1:24),   'uchar');
890pad = zeros(1, 24-length(hist.aux_file));
891hist.aux_file = [hist.aux_file  char(pad)];
892fwrite(fid, hist.aux_file(1:24), 'uchar');
893
894% Write NIfTI style data history
895if any(filetype==[1 2])
896  fwrite(fid, hist.qform_code,    'int16');
897  fwrite(fid, hist.sform_code,    'int16');
898  fwrite(fid, hist.quatern_b,   'float32');
899  fwrite(fid, hist.quatern_c,   'float32');
900  fwrite(fid, hist.quatern_d,   'float32');
901  fwrite(fid, hist.qoffset_x,   'float32');
902  fwrite(fid, hist.qoffset_y,   'float32');
903  fwrite(fid, hist.qoffset_z,   'float32');
904  fwrite(fid, hist.srow_x(1:4), 'float32');
905  fwrite(fid, hist.srow_y(1:4), 'float32');
906  fwrite(fid, hist.srow_z(1:4), 'float32');
907
908  % intent_name = sprintf('%-16s', hist.intent_name);   % 16 chars from left
909  % fwrite(fid, intent_name(1:16),    'uchar');
910  pad = zeros(1, 16-length(hist.intent_name));
911  hist.intent_name = [hist.intent_name  char(pad)];
912  fwrite(fid, hist.intent_name(1:16), 'uchar');
913 
914  % magic       = sprintf('%-4s', hist.magic);          % 4 chars from left
915  % fwrite(fid, magic(1:4),           'uchar');
916  pad = zeros(1, 4-length(hist.magic));
917  hist.magic = [hist.magic  char(pad)];
918  fwrite(fid, hist.magic(1:4),        'uchar');
919elseif filetype==0 % Write Analyze 7.5 style (old) data history
920  fwrite(fid,hist.orient,'uchar');
921  fwrite(fid,hist.originator,'int16');
922  fwrite(fid,hist.generated,'uchar');
923  fwrite(fid,hist.scannum,'uchar');
924  fwrite(fid,hist.patient_id,'uchar');
925  fwrite(fid,hist.exp_date,'uchar');
926  fwrite(fid,hist.exp_time,'uchar');
927  fwrite(fid,hist.hist_un0,'uchar');
928  fwrite(fid,hist.views,'int32');
929  fwrite(fid,hist.vols_added,'int32');
930  fwrite(fid,hist.start_field,'int32');
931  fwrite(fid,hist.field_skip,'int32');
932  fwrite(fid,hist.omax,'int32');
933  fwrite(fid,hist.omin,'int32');
934  fwrite(fid,hist.smax,'int32');
935  fwrite(fid,hist.smin,'int32');
936 
937  % Write also the magic field NIfTI/Analyze75 identifier
938  pad=zeros(1,4-length(hist.magic));
939  hist.magic=[hist.magic char(pad)];
940  fseek(fid,344,'bof');
941  fwrite(fid,hist.magic(1:4),'uchar');
942end
943return;                                 % data_history
944
945
946%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
947% Write NIfTI/Analyze75 file
948%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
949function [done,msg] = l_WriteNiftiData(DATA,hdr,fid,Dat)
950
951done=false;
952msg='';
953
954switch double(hdr.dime.datatype),
955 case   1,
956  precision = 'ubit1';
957 case   2,
958  precision = 'uint8';
959 case   4,
960  precision = 'int16';
961 case   8,
962  precision = 'int32';
963 case  16,
964  precision = 'float32';
965 case  64,
966  precision = 'float64';
967 case 128,
968  precision = 'uint8';
969 case 256
970  precision = 'int8';
971 case 512
972  precision = 'uint16';
973 case 768
974  precision = 'uint32';
975 case 1024
976  precision = 'int64';
977 case 1280
978  precision = 'uint64';
979 otherwise
980  msg='Unsupported datatype.';
981  return
982end
983
984%% if original data is readed from an Analyze75 or NIfTI file, it is very
985%% likely that the data is already calibrated using scl_slope and
986%% scl_interp. Because of this, the data is "decalibrated".
987% ---- On second thought, don't "decalibrate"
988% if isfield(DATA,'DataFormat') && ...
989%     any(strcmpi(DATA.DataFormat,{'analyze75','nifti(1)','nifti(2)'}))
990%   if DATA.HDR.FileHeader.dime.scl_slope~=0
991%     DATA.FTDATA = (DATA.FTDATA-DATA.HDR.FileHeader.dime.scl_inter)/...
992%         DATA.HDR.FileHeader.dime.scl_slope;
993%   end
994% end
995
996
997try
998  ScanDim = double(hdr.dime.dim(5));            % t
999  SliceDim = double(hdr.dime.dim(4));           % z
1000  RowDim   = double(hdr.dime.dim(3));           % y
1001  PixelDim = double(hdr.dime.dim(2));           % x
1002  SliceSz  = double(hdr.dime.pixdim(4));
1003  RowSz    = double(hdr.dime.pixdim(3));
1004  PixelSz  = double(hdr.dime.pixdim(2));
1005 
1006  x = 1:PixelDim;
1007 
1008  if Dat.filetype == 2
1009    skip_bytes = double(hdr.dime.vox_offset) - 348;
1010  else
1011    skip_bytes = 0;
1012  end
1013 
1014  if double(hdr.dime.datatype) == 128
1015    DATA.FTDATA = permute(DATA.FTDATA, [4 1 2 3 5]);
1016  end
1017 
1018  if skip_bytes
1019    fwrite(fid, ones(1,skip_bytes), 'uint8');
1020  end
1021  %flipdim(permute(img,[2 1 3 4]),1)
1022  fwrite(fid,permute(flipdim(DATA.FTDATA,1),[2 1 3 4]), precision);
1023  %   fwrite(fid, nii.img, precision, skip_bytes);        % error using skip
1024  %fclose(fid);
1025catch
1026  msg=lasterr;
1027end
1028
1029done=true;
1030return;                                 % write_nii
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008