function [vol meta] = readmeta(filename)
% READMETA - Reads MetaImage format file from disk
%
% USAGE:
% vol = dzip('file.mha')
% [vol meta] = dzip('file.mha');
%
% VARIABLES:
% filename = filename to open
% vol = image data
% meta = structure containg image metadata such as spacing and orientation

% Created by Casey Goodlett 

fid = fopen(filename,'rb');
  if fid == -1
      error('Could not open MetaImage');
  end
  
  dflag = 1;
  
  binarydata = true;
  compresseddata = false;
  byteorder = 'native';
  
  while dflag
    l = fgetl(fid);
    key = sscanf(l,'%s',1);
    sscanf(l,'%*s %s',1);
    value = sscanf(l,'%*s %*s %s',1);

    if strcmp(key,'ObjectType')
      if ~strcmp(value,'Image')
          error('Only Image ObjectType allowed');
      end
    elseif strcmp(key,'NDims')
      dims = uint32(str2num(value));
    elseif strcmp(key,'BinaryDataByteOrderMSB')
      if strcmp(value,'False')
        byteorder = 'ieee-le';
      else
        byteorder = 'ieee-be';
      end
    elseif strcmp(key,'BinaryData')
      if strcmp(value,'True')
          binarydata = true;
      end
    elseif strcmp(key,'CompressedData')
        if strcmp(value,'True')
            compresseddata = true;
        end
    elseif strcmp(key,'DimSize')
      isize = uint32(sscanf( l(11:end) ,'%d'));
    elseif strcmp(key,'ElementSpacing')
      spacing = sscanf( l(18:end) ,'%f');
    elseif strcmp(key,'Offset')
      origin = sscanf( l(10:end), '%f');
    elseif strcmp(key,'ElementType')
      if strcmp(value,'MET_USHORT')
        etype = 'uint16';
      elseif strcmp(value,'MET_SHORT')
        etype = 'int16';
      elseif strcmp(value,'MET_CHAR')
        etype = 'int8';
      elseif strcmp(value,'MET_UCHAR')
        etype = 'uint8';
      elseif strcmp(value,'MET_INT')
        etype = 'int32';
      elseif strcmp(value,'MET_UINT')
        etype = 'uint32';
      elseif strcmp(value,'MET_LONG')
        etype = 'int64';
      elseif strcmp(value,'MET_ULONG')
        etype = 'uint64';
      elseif strcmp(value,'MET_FLOAT')
        etype = 'single';
      elseif strcmp(value,'MET_DOUBLE')
        etype = 'double';
      else
        error('Unknown pixel type');
      end      
    elseif strcmp(key,'ElementDataFile')
        if ~strcmp(value,'LOCAL')
            fclose(fid);
            fid = fopen(value,'rb');
        end
        dflag = 0;
    end
    
  end
  
  if length(isize) ~= dims
      error('Number sizes != dims');
  end
  nelem = prod(double(isize));

  if binarydata
    if compresseddata
        vol = fread(fid,'uint8');
        vol = deflate(vol);
        vol = typecast(vol, etype);
        [fn pm mf] = fopen(1);
        if ~strcmp(byteorder, mf(1:7))
            vol = byteswap(vol);
        end
    else
        vol = fread(fid,nelem,etype,0,byteorder);
    end
  else
    if strcmp(etype,'float') || strcmp(etype,'double')
      vol = fscanf(fid,'%f',nelem);
    else
      vol = fscanf(fid,'%d',nelem);
    end
  end

  vol = cast(vol,etype);
  vol = reshape(vol,isize');
  
  if nargout > 1
    meta = struct('spacing',{0});
    if exist('spacing')
        %setfield(meta,'spacing', spacing);
        meta.spacing = spacing;
    end
    if exist('origin')
        meta.origin = origin;
    end
  end

function q = deflate(Z);

import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
a=java.io.ByteArrayInputStream(Z);
b=java.util.zip.InflaterInputStream(a);
isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
c = java.io.ByteArrayOutputStream;
isc.copyStream(b,c);

q = typecast(c.toByteArray, 'uint8');