source: an2_write_nifti.m @ 46

Last change on this file since 46 was 46, checked in by tjniskan, 11 years ago
  • Added preliminary support for reading multiple 3D files into Aedes

as 4D data.

  • Added new function an2_getdataformat.m which tries to determine data

format of a file. The code in an2_checkcthdr.m is now part of
an2_getdataformat.m and hence the depricated an2_checkcthdr.m was
removed.

  • Removed depricated function an2_mrdread.m
  • Fixed a small bug in an2_juigetfiles.m that appeared in Linux

version of Matlab.

M an2_getfilefilter.m
D an2_checkcthdr.m
M an2_read_nifti.m
M an2_data_read.m
D an2_mrdread.m
M aedes.m
M an2_juigetfiles.m
M an2_editstack.m
A an2_getdataformat.m
M an2_write_nifti.m

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

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