source: aedes_write_nifti.m

Last change on this file was 211, checked in by tjniskan, 4 years ago
  • Added an option for inputting custom rotation matrix to aedes_write_nifti.m

M aedes_write_nifti.m
M aedes_revision.m

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

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