source: aedes_read_nifti.m

Last change on this file was 80, checked in by tjniskan, 9 years ago
  • Changed the historical "an2_" prefix to "aedes_" in all files. NOTE:

Any script or function relying to Aedes functions will be broken
because of this. Just do a search/replace from "an2_" to "aedes_" in
your files and all should be well...

  • Changed the name of an2_readtab.m to a more informative

aedes_readphasetable.m

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

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