1 function [vol, M, mr_parms, volsz] = load_mgh(fname,slices,frames,headeronly)
   2 % [vol, M, mr_parms, Mdc, volsz] = load_mgh(fname,<slices>,<frames>,<headeronly>)
   3 %
   4 % fname - path of the mgh file
   5 % 
   6 % slices - list of one-based slice numbers to load. All
   7 %   slices are loaded if slices is not specified, or
   8 %   if slices is empty, or if slices(1) <= 0.
   9 %
  10 % frames - list of one-based frame numbers to load. All
  11 %   frames are loaded if frames is not specified, or
  12 %   if frames is empty, or if frames(1) <= 0.
  13 %
  14 % M is the 4x4 vox2ras transform such that
  15 % y(i1,i2,i3), xyz1 = M*[i1 i2 i3 1] where the
  16 % indices are 0-based. If the input has multiple frames,
  17 % only the first frame is read.
  18 %
  19 % mr_parms = [tr flipangle te ti fov]
  20 %
  21 % volsz = size(vol). Helpful when using headeronly as vol is [].
  22 %
  23 % See also: save_mgh, vox2ras_0to1
  24 %
  25 % $Id: load_mgh.m,v 1.14 2006/01/18 20:53:03 fischl Exp $
  26 
  27 vol = [];
  28 M = [];
  29 mr_parms = [];
  30 volsz = [];
  31 
  32 if(nargin < 1 | nargin > 4)
  33   msg = 'USAGE: [vol M] = load_mgh(fname,<slices>,<frames>,<headeronly>)';
  34   fprintf('%s',msg);
  35   return;
  36 end
  37 
  38 % unzip if it is compressed 
  39 if (strcmpi(fname((strlen(fname)-3):strlen(fname)), '.MGZ') | ...
  40 		strcmpi(fname((strlen(fname)-3):strlen(fname)), '.GZ'))
  41 	gzipped =  round(rand(1)*10000000);
  42 	ind = findstr(fname, '.');
  43 	new_fname = sprintf('/tmp/tmp%d.mgh', gzipped);
  44 	unix(sprintf('zcat %s > %s', fname, new_fname)) ;
  45 	fname = new_fname ;
  46 else
  47 	gzipped = -1 ;
  48 end
  49 
  50 
  51 if(exist('slices')~=1) slices = []; end
  52 if(isempty(slices)) slices = 0; end
  53 if(slices(1) <= 0) slices = 0; end
  54 
  55 if(exist('frames')~=1) frames = []; end
  56 if(isempty(frames)) frames = 0; end
  57 if(frames(1) <= 0) frames = 0; end
  58 
  59 if(exist('headeronly')~=1) headeronly = 0; end
  60 
  61 fid    = fopen(fname, 'rb', 'b') ;
  62 if(fid == -1)
  63   fprintf('ERROR: could not open %s for reading\n',fname);
  64   return;
  65 end
  66 v       = fread(fid, 1, 'int') ; 
  67 ndim1   = fread(fid, 1, 'int') ; 
  68 ndim2   = fread(fid, 1, 'int') ; 
  69 ndim3   = fread(fid, 1, 'int') ; 
  70 nframes = fread(fid, 1, 'int') ;
  71 type    = fread(fid, 1, 'int') ; 
  72 dof     = fread(fid, 1, 'int') ; 
  73 
  74 if(slices(1) > 0)
  75   ind = find(slices > ndim3);
  76   if(~isempty(ind))
  77     fprintf('ERROR: load_mgh: some slices exceed nslices\n');
  78     return;
  79   end
  80 end
  81 
  82 if(frames(1) > 0)
  83   ind = find(frames > nframes);
  84   if(~isempty(ind))
  85     fprintf('ERROR: load_mgh: some frames exceed nframes\n');
  86     return;
  87   end
  88 end
  89 
  90 UNUSED_SPACE_SIZE= 256;
  91 USED_SPACE_SIZE = (3*4+4*3*4);  % space for ras transform
  92 
  93 unused_space_size = UNUSED_SPACE_SIZE-2 ;
  94 ras_good_flag = fread(fid, 1, 'short') ; 
  95 if (ras_good_flag)
  96   delta  = fread(fid, 3, 'float32') ; 
  97   Mdc    = fread(fid, 9, 'float32') ; 
  98   Mdc    = reshape(Mdc,[3 3]);
  99   Pxyz_c = fread(fid, 3, 'float32') ; 
 100 
 101   D = diag(delta);
 102 
 103   Pcrs_c = [ndim1/2 ndim2/2 ndim3/2]'; % Should this be kept?
 104 
 105   Pxyz_0 = Pxyz_c - Mdc*D*Pcrs_c;
 106 
 107   M = [Mdc*D Pxyz_0;  ...
 108 	0 0 0 1];
 109   ras_xform = [Mdc Pxyz_c; ...
 110 	0 0 0 1];
 111   unused_space_size = unused_space_size - USED_SPACE_SIZE ;
 112 end
 113 
 114 fseek(fid, unused_space_size, 'cof') ;
 115 nv = ndim1 * ndim2 * ndim3 * nframes;  
 116 volsz = [ndim1 ndim2 ndim3 nframes];
 117 
 118 MRI_UCHAR =  0 ;
 119 MRI_INT =    1 ;
 120 MRI_LONG =   2 ;
 121 MRI_FLOAT =  3 ;
 122 MRI_SHORT =  4 ;
 123 MRI_BITMAP = 5 ;
 124 
 125 % Determine number of bytes per voxel
 126 switch type
 127  case MRI_FLOAT,
 128   nbytespervox = 4;
 129  case MRI_UCHAR,
 130   nbytespervox = 1;
 131  case MRI_SHORT,
 132   nbytespervox = 2;
 133  case MRI_INT,
 134   nbytespervox = 4;
 135 end
 136 
 137 if(headeronly)
 138   fseek(fid,nv*nbytespervox,'cof');
 139   if(~feof(fid))
 140     [mr_parms count] = fread(fid,4,'float32');
 141     if(count ~= 4) 
 142       fprintf('WARNING: error reading MR params\n');
 143     end
 144   end
 145   fclose(fid);
 146   if(gzipped >=0)  unix(sprintf('rm %s', fname));  end
 147   return;
 148 end
 149 
 150 
 151 %------------------ Read in the entire volume ----------------%
 152 if(slices(1) <= 0 & frames(1) <= 0)
 153   switch type
 154    case MRI_FLOAT,
 155     vol = fread(fid, nv, 'float32') ; 
 156    case MRI_UCHAR,
 157     vol = fread(fid, nv, 'uchar') ; 
 158    case MRI_SHORT,
 159     vol = fread(fid, nv, 'short') ; 
 160    case MRI_INT,
 161     vol = fread(fid, nv, 'int') ; 
 162   end
 163 
 164   if(~feof(fid))
 165     [mr_parms count] = fread(fid,4,'float32');
 166     if(count ~= 4) 
 167       fprintf('WARNING: error reading MR params\n');
 168     end
 169   end
 170   fclose(fid) ;
 171   if(gzipped >=0)  unix(sprintf('rm %s', fname));  end
 172   
 173   nread = prod(size(vol));
 174   if(nread ~= nv)
 175     fprintf('ERROR: tried to read %d, actually read %d\n',nv,nread);
 176     vol = [];
 177     return;
 178   end
 179   vol = reshape(vol,[ndim1 ndim2 ndim3 nframes]);
 180 
 181   return;
 182 end
 183 
 184 %----- only gets here if a subest of slices/frames are to be loaded ---------%
 185 
 186 
 187 if(frames(1) <= 0) frames = [1:nframes]; end
 188 if(slices(1) <= 0) slices = [1:ndim3]; end
 189 
 190 nvslice = ndim1 * ndim2;
 191 nvvol   = ndim1 * ndim2 * ndim3;
 192 filepos0 = ftell(fid);
 193 vol = zeros(ndim1,ndim2,length(slices),length(frames));
 194 nthframe = 1;
 195 for frame = frames
 196 
 197   nthslice = 1;
 198   for slice = slices
 199     filepos = ((frame-1)*nvvol + (slice-1)*nvslice)*nbytespervox + filepos0;
 200     fseek(fid,filepos,'bof');
 201     
 202     switch type
 203      case MRI_FLOAT,
 204       [tmpslice nread]  = fread(fid, nvslice, 'float32') ; 
 205      case MRI_UCHAR,
 206       [tmpslice nread]  = fread(fid, nvslice, 'uchar') ; 
 207      case MRI_SHORT,
 208       [tmpslice nread]  = fread(fid, nvslice, 'short') ; 
 209      case MRI_INT,
 210       [tmpslice nread]  = fread(fid, nvslice, 'int') ; 
 211     end
 212 
 213     if(nread ~= nvslice)
 214       fprintf('ERROR: load_mgh: reading slice %d, frame %d\n',slice,frame);
 215       fprintf('  tried to read %d, actually read %d\n',nvslice,nread);
 216       fclose(fid);
 217       if(gzipped >=0) unix(sprintf('rm %s', fname)); end
 218       return;
 219     end
 220 
 221     vol(:,:,nthslice,nthframe) = reshape(tmpslice,[ndim1 ndim2]);
 222     nthslice = nthslice + 1;
 223   end
 224 
 225   nthframe = nthframe + 1;
 226 end
 227 
 228 % seek to just beyond the last slice/frame %
 229 filepos = (nframes*nvvol)*nbytespervox + filepos0;
 230 fseek(fid,filepos,'bof');
 231 
 232 if(~feof(fid))
 233   [mr_parms count] = fread(fid,5,'float32');
 234   if(count < 4) 
 235     fprintf('WARNING: error reading MR params\n');
 236   end
 237 end
 238 
 239 fclose(fid) ;
 240 if(gzipped >=0) unix(sprintf('rm %s', fname)); end
 241 
 242 return;