source: an2_write_nifti.m @ 41

Last change on this file since 41 was 39, checked in by tjniskan, 11 years ago
  • Added "Recent Files" quick access menu in File-menu
  • Fixed a bug in plugins/copy_data_to_workspace.m that made writing a

new name for a variable difficult

  • Fixed a bug in the file prefix editbox in an2_export_gui.m
  • Changed references to Analyze 2.0 in all license notices to Aedes.

The name change should now be complete...

M an2_export_gui.m
M an2_cellsprintf.m
M an2_calc_wait.m
M an2_check_file_exist.m
M an2_iseven.m
M an2_cellwrite.m
M an2_wbar.m
M an2_rot3d.m
M an2_readfdf.m
M an2_revision.m
M an2_viewprocpar.m
M an2_checkcthdr.m
M an2_readprocpar.m
M an2_fitmaps.m
M an2_read_nifti.m
M an2_data_read.m
M an2_resviewer.m
M an2_maptool.m
M aedes.m
M an2_res2table.m
M an2_copy_roi.m
M plugins/save_roi_as_mask.m
M plugins/write_difference_images.m
M plugins/plot_profile.m
M plugins/calculate_t2_map.m
M plugins/calculate_t1r_map.m
M plugins/view_kspace.m
M plugins/copy_data_to_workspace.m
M plugins/take_snapshot.m
M an2_inputdlg.m
M an2_roi_copy_gui.m
M an2_readctdata.m
M an2_readfid.m
M an2_readfidprefs.m
M an2_readtab.m
M an2_check_updates.m
M an2_killfigs.m
M an2_roi_stats.m
M an2_saveres.m
M an2_rotateflip.m
M an2_juigetfiles.m
M an2_gui_defaults.m
M an2_editstack.m
M an2_update.m
M an2_write_nifti.m

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

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