source: an2_write_nifti.m @ 44

Last change on this file since 44 was 42, checked in by tjniskan, 11 years ago
  • Rewrote the overlay feature of Aedes to work with absolute rather

than relative values.

  • The export tool can now also print overlays and ROIs with

transparencies. In these cases, however, it has to use OpenGL renderer
which is more unreliable in terms of print quality than other renderers.

  • Added a license notification that prints into the command window

every time Aedes is started. This has to do with the coming
"internationalization" of Aedes. The notification can be suppressed y
typing setpref('Aedes','ShowLicenseAtStartUp?',false) into the command
window.

  • Added some further error checking in an2_fitmaps.m
  • Fixed a minor bug in an2_write_nifti.m

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

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