source: an2_read_nifti.m @ 44

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

new name for a variable difficult

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

The name change should now be complete...

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

File size: 23.4 KB
Line 
1function [DATA,msg] = an2_read_nifti(filename,varargin)
2% AN2_READ_NIFTI - Read NIfTI (*.nii) and Analyze 7.5 (*.hdr,*.img) files
3%   
4%
5% Synopsis:
6%       [DATA,msg]=an2_read_nifti(filename,opt)
7%       [DATA,msg]=an2_read_nifti(DATA.HDR)
8%
9% Description:
10%       The function reads the NIfTI and Analyze 7.5 file formats and
11%       returns data in DATA.FTDATA and header in DATA.HDR (for more
12%       information about the structure of DATA, see the help of AN2_READFID
13%       function, i.e. type "help an2_readfid"). If error has occurred during
14%       reading, the DATA structure is empty and the second output
15%       argument msg contains the error message. The function is capable
16%       of reading NIfTI and Analyze 7.5 files in the two file
17%       format (*.img and *.hdr) and the NIfTI on file format (*.nii).
18%
19%       First input argument can be either the full path to the NIfTI or
20%       Analyze 7.5 file or the NIfTI/Analyze7.5 header structure. If the
21%       second input argument is a string 'header', only the header of
22%       the file is read. If the function is called without input
23%       arguments, a file dialog will appear to browse for a file.
24%
25%
26%       Acknowledgment: This function is modified under GNU license from
27%       MRI_TOOLBOX developed by CNSP in Flinders University,
28%       Australia. Some parts are also originally written by Jimmy Shen.
29%       See the following links:
30%       http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8797&objectType=file
31%       http://eeg.sourceforge.net/
32
33%       NIfTI data format can be found on: http://nifti.nimh.nih.gov
34%
35% Examples:
36%       [DATA,msg]=an2_read_nifti(filename,'header'); % Read only header
37%       [DATA,msg]=an2_read_nifti(DATA.HDR);          % Read data
38%       or
39%       [DATA,msg]=an2_read_nifti;                    % Browse for a file and
40%                                                 % read both header and
41%                                                 % data
42%       or
43%       [DATA,msg]=an2_read_nifti(filename)           % Read both header and
44%                                                 % data
45%
46% See also:
47%       AN2_READFID, AEDES, AN2_WRITE_NIFTI
48
49% This function is a part of Aedes - A graphical tool for analyzing
50% medical images
51%
52% Copyright (C) 2006 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uku.fi>
53% Copyright (C) 2006 Jimmy Shen
54%
55% Department of Physics, (or Department of Neurobiology)
56% University of Kuopio, FINLAND
57%
58% This program may be used under the terms of the GNU General Public
59% License version 2.0 as published by the Free Software Foundation
60% and appearing in the file LICENSE.TXT included in the packaging of
61% this program.
62%
63% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
64% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
65
66
67DATA=[];
68data=[];
69msg='';
70
71ReadHdr = true;
72ReadData = true;
73machine=[];
74
75%% Parse input arguments
76switch nargin
77 
78  %% Ask for a file name, if not provided in the input argument
79 case 0,
80  [f_name, f_path, f_index] = uigetfile( ...
81       {'*.nii;*.NII;*.nii.gz;*.hdr;*.HDR',...
82        'NIfTI and Analyze 7.5 files (*.nii,*.nii.gz,*.hdr)'; ...
83        '*.*', 'All Files (*.*)'}, ...
84        'Select NIfTI or Analyze 7.5 files');
85  if ( all(f_name==0) | all(f_path==0) ) % Cancel is pressed
86    DATA=[];
87    msg_out='Canceled';
88    return
89  end
90  filename=[f_path,f_name];
91 
92 case 1,
93  if isstruct(filename)
94    hdr=filename;
95    filename=[hdr.fpath,hdr.fname];
96        machine = hdr.byteorder;
97    ReadHdr = false;
98    ReadData = true;
99  else
100    ReadHdr = true;
101    ReadData = true;
102  end
103 case 2,
104  % Read only header
105  if strcmpi(varargin{1},'header')
106    ReadHdr = true;
107    ReadData = false;
108  else
109    msg=['Unknown option "' varargin{1} '".'];
110    return
111  end
112 otherwise
113  msg='Too many input arguments';
114  return
115end
116
117%% Parse filename
118[f_path,f_name,f_ext] = fileparts(filename);
119if isempty(f_path)
120  % If there is no path information in "filename", do some extra stuff
121  % here...
122  filename = which(filename);
123  [f_path,f_name,f_ext] = fileparts(filename);
124end
125
126%% Check if the file is gzipped NIfTI (*.nii.gz)
127gz_filename = '';
128if strcmpi(f_ext,'.gz') && length(f_name)>3 && ...
129        strcmpi(f_name(end-3:end),'.nii')
130 
131  % Try to gunzip the file to a temporary folder and read it from there
132  try
133        gunzip(filename,tempdir);
134  catch
135        error('Could not gunzip "%s" to "%s"',filename,tempdir)
136  end
137  gz_filename = fullfile(tempdir,f_name);
138  [f_path,f_name,f_ext]=fileparts(gz_filename);
139end
140
141%% Read header
142if ReadHdr
143  [hdr,tmp_msg,machine]=l_ReadHdr([f_path,filesep,f_name],f_ext,machine);
144  if isempty(hdr)
145    DATA=[];
146    if iscell(tmp_msg)
147      msg = {'Error while reading data header.',...
148             tmp_msg{:}};
149    else
150      msg = {'Error while reading data header.',...
151             tmp_msg};
152    end
153    return
154  end
155end
156
157%% Set dataformat string
158if strcmp(hdr.FileHeader.hist.magic,'n+1')
159  DATA.DataFormat = 'NIfTI(1)'; % NIfTI format (*.nii), one file
160  filetype=2;
161elseif strcmp(hdr.FileHeader.hist.magic,'ni1')
162  DATA.DataFormat = 'NIfTI(2)'; % NIfTI format (*.hdr,*.img), two files
163  filetype=1;
164else
165  DATA.DataFormat = 'Analyze75'; % Analyze 7.5 format (*.hdr,*.img), two
166                                % files
167  filetype=0;
168end
169
170%% Read data
171if ReadData
172  [data,hdr,tmp_msg]=l_ReadData(hdr,machine,[f_path,filesep,f_name,f_ext],filetype);
173  if isempty(data)
174    DATA=[];
175    if iscell(tmp_msg)
176      msg = {'Error while reading data.',...
177             tmp_msg{:}};
178    else
179      msg = {'Error while reading data.',...
180             tmp_msg};
181    end
182    return
183  end
184end
185
186%% Set fields to DATA
187if isempty(gz_filename)
188  DATA.HDR.FileHeader = hdr.FileHeader;
189  DATA.HDR.fname = [f_name,f_ext];
190  DATA.HDR.fpath = [f_path,filesep];
191  DATA.HDR.byteorder = machine;
192else
193  % In .nii.gz files, use the original filename
194  [f_path,f_name,f_ext] = fileparts(filename);
195 
196  % Clean up -> remove the temporary gunzipped file
197  try
198        delete(gz_filename)
199  catch
200        warning(['Could not remove temporary file "%s". ',...
201          'Maybe you should see what''s wrong...'],gz_filename)
202  end
203 
204  DATA.HDR.FileHeader = hdr.FileHeader;
205  DATA.HDR.fname = [f_name,f_ext];
206  DATA.HDR.fpath = [f_path,filesep];
207  DATA.HDR.byteorder = machine;
208end
209
210
211% Set xyz and time units
212tmp=dec2bin(DATA.HDR.FileHeader.dime.xyzt_units,8);
213
214xyzunits = tmp(6:8);
215timeunits = tmp(3:5);
216
217% Units of spatial and temporal dimensions:
218% 0 = Unknown, 1 = meter, 2 = millimeter, 3 = micrometer, 8 = seconds
219% 16 = milliseconds, 24 = microseconds, 32 = Hertz, 40 = ppm,
220% 48 = rad/sec.
221% Bits 0..2 of xyzt_units specify the units of pixdim[1..3]. Bits
222% 3..5 of xyzt_units specify the units of pixdim[4] and are multiples
223% of 8.   
224if bin2dec(xyzunits)==0
225  DATA.HDR.xyzunits = 'Unknown';
226elseif bin2dec(xyzunits)==1
227  DATA.HDR.xyzunits = 'meter';
228elseif bin2dec(xyzunits)==2
229  DATA.HDR.xyzunits = 'mm';
230elseif bin2dec(xyzunits)==3
231  DATA.HDR.xyzunits = 'micron';
232end
233if bin2dec(timeunits)==0
234  DATA.HDR.timeunits = 'Unknown';
235elseif bin2dec(timeunits)==1
236  DATA.HDR.timeunits = 'sec';
237elseif bin2dec(timeunits)==2
238  DATA.HDR.timeunits = 'msec';
239elseif bin2dec(timeunits)==3
240  DATA.HDR.timeunits = 'usec';
241elseif bin2dec(timeunits)==4
242  DATA.HDR.timeunits = 'Hz';
243elseif bin2dec(timeunits)==5
244  DATA.HDR.timeunits = 'ppm';
245elseif bin2dec(timeunits)==6
246  DATA.HDR.timeunits = 'rad/sec';
247end
248if ReadData
249  DATA.FTDATA = data;
250  DATA.KSPACE=[];
251end
252
253
254%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255%%
256%% Read NIfTI/Analyze75 header
257%%
258%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
259function [hdr,msg,machine]=l_ReadHdr(filename,ext,machine)
260
261hdr=[];
262msg='';
263
264%% Try first little-endian byte ordering
265if isempty(machine)
266  machine = 'ieee-le';
267end
268
269new_ext = 0;
270
271if strcmpi(ext,'.nii')
272  new_ext = 1;
273else
274  new_ext = 0;
275end
276
277if any(strcmpi(ext,{'.hdr','.nii'}))
278  fname = [filename,ext];
279elseif strcmpi(ext,'.img')
280  fname = [filename,'.hdr'];
281end
282
283
284fid = fopen(fname,'r',machine);
285if fid < 0,
286  hdr=[];
287  msg = sprintf('Cannot open file %s.',fname);
288  return
289end
290
291%% Try to determine if the header file is Analyze75 or NIfTI
292fseek(fid,344,'bof');
293tmp=fread(fid,4,'*char');
294tmp=deblank(tmp');
295if isempty(tmp)
296  filetype=0; % Analyze 7.5
297elseif strcmpi(tmp,'ni1')
298  filetype=1; % Two file NIfTI
299elseif strcmpi(tmp,'n+1')
300  filetype=2; % One file NIfTI
301else
302  % Somethings wrong -> throw an error
303  fclose(fid);
304  hdr=[];
305  msg = sprintf(['The file "%s" is not a valid NIfTI/Analyze 7.5 header ' ...
306                 'file.'],fname);
307  return
308end
309
310% Rewind file
311fseek(fid,0,'bof');
312
313% Read header
314hdr.FileHeader = read_header(fid,filetype);
315fclose(fid);
316
317
318
319
320
321% Check if machine format was correct
322if hdr.FileHeader.hk.sizeof_hdr ~= 348
323  % first try reading the opposite endian to 'machine'
324  switch machine,
325   case 'ieee-le', machine = 'ieee-be';
326   case 'ieee-be', machine = 'ieee-le';
327  end
328 
329  fid = fopen(fname,'r',machine);
330  if fid < 0,
331    hdr=[];
332    msg = sprintf('Cannot open file %s.',fname);
333  else
334    hdr.FileHeader = read_header(fid,filetype);
335    fclose(fid);
336  end
337end
338
339
340if hdr.FileHeader.hk.sizeof_hdr ~= 348
341  % Now throw an error
342  hdr=[];
343  msg = sprintf('File "%s" is probably corrupted.',fname);
344end
345
346% $$$ if strcmp(hdr.hist.magic, 'n+1')
347% $$$   filetype = 2;
348% $$$ elseif strcmp(hdr.hist.magic, 'ni1')
349% $$$   filetype = 1;
350% $$$ else
351% $$$   filetype = 0;
352% $$$ end
353
354return                                  % load_nii_hdr
355
356
357%---------------------------------------------------------------------
358function [ dsr ] = read_header(fid,filetype)
359
360%  Original header structures
361%  struct dsr
362%       {
363%       struct header_key hk;            /*   0 +  40       */
364%       struct image_dimension dime;     /*  40 + 108       */
365%       struct data_history hist;        /* 148 + 200       */
366%       };                               /* total= 348 bytes*/
367
368dsr.hk   = header_key(fid);
369dsr.dime = image_dimension(fid);
370dsr.hist = data_history(fid,filetype);
371
372%  For Analyze data format
373%
374% $$$ if ~strcmp(dsr.hist.magic, 'n+1') & ~strcmp(dsr.hist.magic, 'ni1')
375% $$$   dsr.dime.scl_slope = 0;
376% $$$   dsr.hist.qform_code = 0;
377% $$$   dsr.hist.sform_code = 0;
378% $$$ end
379
380return                                  % read_header
381
382
383%---------------------------------------------------------------------
384function [ hk ] = header_key(fid)
385
386fseek(fid,0,'bof');
387
388%  Original header structures   
389%  struct header_key                     /* header key      */
390%       {                                /* off + size      */
391%       int sizeof_hdr                   /*  0 +  4         */
392%       char data_type[10];              /*  4 + 10         */
393%       char db_name[18];                /* 14 + 18         */
394%       int extents;                     /* 32 +  4         */
395%       short int session_error;         /* 36 +  2         */
396%       char regular;                    /* 38 +  1         */
397%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
398%       };                               /* total=40 bytes  */
399%
400% int sizeof_header   Should be 348.
401% char regular        Must be 'r' to indicate that all images and
402%                     volumes are the same size.
403
404hk.sizeof_hdr    = fread(fid, 1,'int32')';      % should be 348!
405hk.data_type     = deblank(fread(fid,10,'*char')');
406hk.db_name       = deblank(fread(fid,18,'*char')');
407hk.extents       = fread(fid, 1,'int32')';
408hk.session_error = fread(fid, 1,'int16')';
409hk.regular       = fread(fid, 1,'*char')';
410hk.dim_info      = fread(fid, 1,'char')';
411
412return                                  % header_key
413
414
415%---------------------------------------------------------------------
416function [ dime ] = image_dimension(fid)
417
418%  Original header structures   
419%  struct image_dimension
420%       {                                /* off + size      */
421%       short int dim[8];                /* 0 + 16          */
422%       /*
423%           dim[0]      Number of dimensions in database; usually 4.
424%           dim[1]      Image X dimension;  number of *pixels* in an image row.
425%           dim[2]      Image Y dimension;  number of *pixel rows* in slice.
426%           dim[3]      Volume Z dimension; number of *slices* in a volume.
427%           dim[4]      Time points; number of volumes in database
428%       */
429%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
430%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
431%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
432%       short int intent_code;   % short int unused1;   /* 28 + 2 */
433%       short int datatype;              /* 30 + 2          */
434%       short int bitpix;                /* 32 + 2          */
435%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
436%       float pixdim[8];                 /* 36 + 32         */
437%       /*
438%               pixdim[] specifies the voxel dimensions:
439%               pixdim[1] - voxel width, mm
440%               pixdim[2] - voxel height, mm
441%               pixdim[3] - slice thickness, mm
442%               pixdim[4] - volume timing, in msec
443%                                       ..etc
444%       */
445%       float vox_offset;                /* 68 + 4          */
446%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
447%       float scl_inter;   % float funused1;      /* 76 + 4 */
448%       short slice_end;   % float funused2;      /* 80 + 2 */
449%       char slice_code;   % float funused2;      /* 82 + 1 */
450%       char xyzt_units;   % float funused2;      /* 83 + 1 */
451%       float cal_max;                   /* 84 + 4          */
452%       float cal_min;                   /* 88 + 4          */
453%       float slice_duration;   % int compressed; /* 92 + 4 */
454%       float toffset;   % int verified;          /* 96 + 4 */
455%       int glmax;                       /* 100 + 4         */
456%       int glmin;                       /* 104 + 4         */
457%       };                               /* total=108 bytes */
458
459dime.dim        = fread(fid,8,'int16')';
460dime.intent_p1  = fread(fid,1,'float32')';
461dime.intent_p2  = fread(fid,1,'float32')';
462dime.intent_p3  = fread(fid,1,'float32')';
463dime.intent_code = fread(fid,1,'int16')';
464dime.datatype   = fread(fid,1,'int16')';
465dime.bitpix     = fread(fid,1,'int16')';
466dime.slice_start = fread(fid,1,'int16')';
467dime.pixdim     = fread(fid,8,'float32')';
468dime.vox_offset = fread(fid,1,'float32')';
469dime.scl_slope  = fread(fid,1,'float32')';
470dime.scl_inter  = fread(fid,1,'float32')';
471dime.slice_end  = fread(fid,1,'int16')';
472dime.slice_code = fread(fid,1,'char')';
473dime.xyzt_units = fread(fid,1,'char')';
474dime.cal_max    = fread(fid,1,'float32')';
475dime.cal_min    = fread(fid,1,'float32')';
476dime.slice_duration = fread(fid,1,'float32')';
477dime.toffset    = fread(fid,1,'float32')';
478dime.glmax      = fread(fid,1,'int32')';
479dime.glmin      = fread(fid,1,'int32')';
480
481return                                  % image_dimension
482
483
484%---------------------------------------------------------------------
485function [ hist ] = data_history(fid,filetype)
486
487%  Original header structures /* (NIfTI format) */ 
488%  struct data_history 
489%       {                                /* off + size      */
490%       char descrip[80];                /* 0 + 80          */
491%       char aux_file[24];               /* 80 + 24         */
492%       short int qform_code;            /* 104 + 2         */
493%       short int sform_code;            /* 106 + 2         */
494%       float quatern_b;                 /* 108 + 4         */
495%       float quatern_c;                 /* 112 + 4         */
496%       float quatern_d;                 /* 116 + 4         */
497%       float qoffset_x;                 /* 120 + 4         */
498%       float qoffset_y;                 /* 124 + 4         */
499%       float qoffset_z;                 /* 128 + 4         */
500%       float srow_x[4];                 /* 132 + 16        */
501%       float srow_y[4];                 /* 148 + 16        */
502%       float srow_z[4];                 /* 164 + 16        */
503%       char intent_name[16];            /* 180 + 16        */
504%       char magic[4];   % int smin;     /* 196 + 4         */
505%       };                               /* total=200 bytes */
506
507hist.descrip     = deblank(fread(fid,80,'*char')');
508hist.aux_file    = deblank(fread(fid,24,'*char')');
509
510%% Read NIfTI fields
511if any(filetype==[1 2])
512  hist.qform_code  = fread(fid,1,'int16')';
513  hist.sform_code  = fread(fid,1,'int16')';
514  hist.quatern_b   = fread(fid,1,'float32')';
515  hist.quatern_c   = fread(fid,1,'float32')';
516  hist.quatern_d   = fread(fid,1,'float32')';
517  hist.qoffset_x   = fread(fid,1,'float32')';
518  hist.qoffset_y   = fread(fid,1,'float32')';
519  hist.qoffset_z   = fread(fid,1,'float32')';
520  hist.srow_x      = fread(fid,4,'float32')';
521  hist.srow_y      = fread(fid,4,'float32')';
522  hist.srow_z      = fread(fid,4,'float32')';
523  hist.intent_name = deblank(fread(fid,16,'*char')');
524  hist.magic       = deblank(fread(fid,4,'*char')');
525else % Read Analyze 7.5 fields
526  hist.orient      = fread(fid, 1,'*char');
527  hist.originator  = fread(fid,5,'*int16')';
528  hist.generated   = fread(fid,10,'*char')';
529  hist.scannum     = fread(fid,10,'*char')';
530  hist.patient_id  = fread(fid,10,'*char')';
531  hist.exp_date    = fread(fid,10,'*char')';
532  hist.exp_time    = fread(fid,10,'*char')';
533  hist.hist_un0    = fread(fid, 3,'*char')';
534  hist.views       = fread(fid, 1,'*int32');
535  hist.vols_added  = fread(fid, 1,'*int32');
536  hist.start_field = fread(fid, 1,'*int32');
537  hist.field_skip  = fread(fid, 1,'*int32');
538  hist.omax        = fread(fid, 1,'*int32');
539  hist.omin        = fread(fid, 1,'*int32');
540  hist.smax        = fread(fid, 1,'*int32');
541  hist.smin        = fread(fid, 1,'*int32');
542 
543  % Read also the magic field
544  fseek(fid,344,'bof');
545  hist.magic       = deblank(fread(fid,4,'*char')');
546end
547
548fseek(fid,253,'bof');
549hist.originator  = fread(fid, 5,'int16')';
550
551return                                  % data_history
552
553
554
555%%%%%%%%%%%%%%%%%%%%%%%%%%%%
556%%
557%% Read NIfTI data
558%%
559%%%%%%%%%%%%%%%%%%%%%%%%%%%%
560function [data,hdr,msg]=l_ReadData(hdr,machine,filename,filetype)
561
562img_idx = [];
563old_RGB = 0;
564msg='';
565
566%% Use little endian byte ordering by default
567if isempty(machine)
568  machine='ieee-le';
569end
570
571[fpath,fname,fext]=fileparts(filename);
572if any(filetype==[0 1])
573  filename = [fpath,filesep,fname,'.img'];
574 
575  %% Make sure that the *.img file is found
576  if exist(filename,'file')~=2
577    data=[];
578    hdr=[];
579    msg = {'Could not find data file',...
580           ['"' filename '"']};
581    return
582  end
583else
584  filename = [fpath,filesep,fname,'.nii'];
585end
586 
587try
588  [data,hdr.FileHeader] = read_image(hdr.FileHeader,filetype,filename,machine,img_idx, ...
589                          old_RGB);
590catch
591  data=[];
592  hdr=[];
593  msg = lasterr;
594end
595
596return
597
598
599%---------------------------------------------------------------------
600function [img,hdr] = read_image(hdr, filetype,filename,machine,img_idx,old_RGB)
601
602%switch filetype
603% case {0, 1}
604%  fn = [fileprefix '.img'];
605% case 2
606%  fn = [fileprefix '.nii'];
607%end
608
609fid = fopen(filename,'r',machine);
610
611if fid < 0,
612  msg = sprintf('Cannot open file %s.',filename);
613  error(msg);
614end
615
616%  Set bitpix according to datatype
617%
618%  /*Acceptable values for datatype are*/
619%
620%     0 None                   (Unknown bit per voxel)   % DT_NONE, DT_UNKNOWN
621%     1 Binary                       (ubit1, bitpix=1)   % DT_BINARY
622%     2 Unsigned char       (uchar or uint8, bitpix=8)   % DT_UINT8, NIFTI_TYPE_UINT8
623%     4 Signed short                 (int16, bitpix=16)  % DT_INT16, NIFTI_TYPE_INT16
624%     8 Signed integer               (int32, bitpix=32)  % DT_INT32, NIFTI_TYPE_INT32
625%    16 Floating point   (single or float32, bitpix=32)  % DT_FLOAT32, NIFTI_TYPE_FLOAT32
626%    32 Complex, 2 float32     (Unsupported, bitpix=64)  % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
627%    64 Double precision (double or float64, bitpix=64)  % DT_FLOAT64, NIFTI_TYPE_FLOAT64
628%   128 uint8 RGB                (Use uint8, bitpix=24)  % DT_RGB24, NIFTI_TYPE_RGB24
629%   256 Signed char          (schar or int8, bitpix=8)   % DT_INT8, NIFTI_TYPE_INT8
630%   511 Single RGB             (Use float32, bitpix=96)  % DT_RGB96, NIFTI_TYPE_RGB96
631%   512 Unsigned short              (uint16, bitpix=16)  % DT_UNINT16, NIFTI_TYPE_UNINT16
632%   768 Unsigned integer            (uint32, bitpix=32)  % DT_UNINT32, NIFTI_TYPE_UNINT32
633%  1024 Signed long long             (int64, bitpix=64)  % DT_INT64, NIFTI_TYPE_INT64
634%  1280 Unsigned long long          (uint64, bitpix=64)  % DT_UINT64, NIFTI_TYPE_UINT64
635%  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
636%  1792 Complex128, 2 float64  (Unsupported, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
637%  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
638%
639switch hdr.dime.datatype
640 case   1,
641  hdr.dime.bitpix = 1;  precision = 'ubit1';
642 case   2,
643  hdr.dime.bitpix = 8;  precision = 'uint8';
644 case   4,
645  hdr.dime.bitpix = 16; precision = 'int16';
646 case   8,
647  hdr.dime.bitpix = 32; precision = 'int32';
648 case  16,
649  hdr.dime.bitpix = 32; precision = 'float32';
650 case  64,
651  hdr.dime.bitpix = 64; precision = 'float64';
652 case 128,
653  hdr.dime.bitpix = 24; precision = 'uint8';
654 case 256
655  hdr.dime.bitpix = 8;  precision = 'int8';
656 case 511
657  hdr.dime.bitpix = 96; precision = 'float32';
658 case 512
659  hdr.dime.bitpix = 16; precision = 'uint16';
660 case 768
661  hdr.dime.bitpix = 32; precision = 'uint32';
662 case 1024
663  hdr.dime.bitpix = 64; precision = 'int64';
664 case 1280
665  hdr.dime.bitpix = 64; precision = 'uint64';
666 otherwise
667  error('This datatype is not supported');
668end
669
670%  move pointer to the start of image block
671%
672switch filetype
673 case {0, 1}
674  fseek(fid, 0, 'bof');
675 case 2
676  fseek(fid, hdr.dime.vox_offset, 'bof');
677end
678
679%  Load whole image block for old Analyze format, or binary image,
680%  or img_idx is empty; otherwise, load images that are specified
681%  in img_idx
682%
683%  For binary image, we have to read all because pos can not be
684%  seeked in bit and can not be calculated the way below.
685%
686if filetype == 0 | hdr.dime.datatype == 1 | isempty(img_idx)
687  img = fread(fid,inf,sprintf('*%s',precision));
688else
689  img = [];
690 
691  for i=1:length(img_idx)
692   
693    %  Precision of value will be read in img_siz times, where
694    %  img_siz is only the dimension size of an image, not the
695    %  byte storage size of an image.
696    %
697    img_siz = prod(hdr.dime.dim(2:4));
698
699    %  Position is seeked in bytes. To convert dimension size
700    %  to byte storage size, hdr.dime.bitpix/8 will be
701    %  applied.
702    %
703    pos = (img_idx(i) - 1) * img_siz * hdr.dime.bitpix/8;
704   
705    if filetype == 2
706      fseek(fid, pos + hdr.dime.vox_offset, 'bof');
707    else
708      fseek(fid, pos, 'bof');
709    end
710
711    %  fread will read precision of value in img_siz times
712    %
713    img = [img fread(fid,img_siz,sprintf('*%s',precision))];
714  end
715end
716
717fclose(fid);
718
719%  Update the global min and max values
720%hdr.dime.glmax = max(double(img(:)));
721%hdr.dime.glmin = min(double(img(:)));
722
723
724if isempty(img_idx)
725  img_idx = 1:hdr.dime.dim(5);
726end
727
728if old_RGB & hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
729  img = squeeze(reshape(img, [hdr.dime.dim(2:3) 3 hdr.dime.dim(4) length(img_idx)]));
730  img = permute(img, [1 2 4 3 5]);
731elseif hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
732  img = squeeze(reshape(img, [3 hdr.dime.dim(2:4) length(img_idx)]));
733  img = permute(img, [2 3 4 1 5]);
734elseif hdr.dime.datatype == 511 & hdr.dime.bitpix == 96
735  img = double(img);
736  img = (img - min(img))/(max(img) - min(img));
737  img = squeeze(reshape(img, [3 hdr.dime.dim(2:4) length(img_idx)]));
738  img = permute(img, [2 3 4 1 5]);
739else
740  img = squeeze(reshape(img, [hdr.dime.dim(2:4) length(img_idx)]));
741  img = flipdim(permute(img,[2 1 3 4]),1);
742end
743
744if ~isempty(img_idx)
745  hdr.dime.dim(5) = length(img_idx);
746end
747
748%% Scale data if requested
749if isfield(hdr.dime,'scl_slope') && hdr.dime.scl_slope~=0
750  if isfield(hdr.dime,'scl_inter')
751        if hdr.dime.scl_slope~=1 && hdr.dime.scl_inter~=0
752          img = double(img);
753          img = hdr.dime.scl_slope*img+hdr.dime.scl_inter;
754        end
755  else
756        if hdr.dime.scl_slope~=1
757          img = double(img);
758          img = hdr.dime.scl_slope*img;
759        end
760  end
761end
762
763return                                          % read_image
Note: See TracBrowser for help on using the repository browser.

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