source: an2_write_nifti.m @ 53

Last change on this file since 53 was 53, checked in by tjniskan, 11 years ago
  • Fixed a bug in an2_readfid.m that in some cases didn't calculate

the array length in procpar correctly. Also changed the handling of
array parameter in an2_readprocpar.m.

  • Fixed a bug in an2_rot3d.m
  • Fixed bugs here and there
  • Added some incomplete code for upcoming new features.

M misclib/nifti4dto3d.m
M misclib/fmri_corr.m
M an2_rot3d.m
M an2_revision.m
M an2_readprocpar.m
M an2_data_read.m
M aedes.m
M an2_res2table.m
M an2_copy_roi.m
M an2_roi_copy_gui.m
M an2_readfid.m
M an2_saveres.m
M an2_write_nifti.m

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

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