utilities_ESRF/edfread.m

   1 ## Copyright (C) 2007 P. Cloetens
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, write to the Free Software
  15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  16
  17 ## edfread
  18 ## function im=edfread(filename,varargin)
  19 ##      reads an image in edf format into the matrix im
  20 ##      uses pmedf_read by Petr Mikulik if only filename is given
  21 ##      if both header and image are needed, use [hd,im]=pmedf_read(filename)
  22 ##      does not transpose the image, 1st dimension remains 1st dimension of original file
  23 ##
  24 ##      arguments:
  25 ##      argument 1: filename
  26 ##      argument 2 (optional): vector of adjacent rows to be read
  27 ##      argument 3 (optional): vector of adjacent colums to be read
  28 ##      argument 4 (optional): vector of layers to be read
  29 ##
  30 ##      examples:
  31 ##      im = edfread('dark.edf');
  32 ##          reads complete image
  33 ##      im = edfread('dark.edf',10:200,30:50);
  34 ##          reads subimage corresponding to row (Dim_1) 10 until 200 and column (Dim_2) 30 until 50
  35 ##      im = edfread('dark.edf',10:200,30:50,[1 4 5]);
  36 ##          reads layers 1, 4 and 5 of a 3-dimensional image (see Dim_3 and W. Ludwig)
  37 ##
  38 ##      See also: edfwrite, pmedf_read, pmedf_write
  39
  40 ## Author: P. Cloetens <cloetens@esrf.fr>
  41 ##
  42 ## 2007-03-03 P. Cloetens <cloetens@esrf.fr>
  43 ## * Initial revision
  44
  45 function im=edfread(filename,varargin)
  46
  47 usepm = 0;
  48
  49 switch nargin
  50     case 1
  51         usepm = isoctave&exist('pmedf_read');
  52     case 3
  53         rows = varargin{1};
  54         columns = varargin{2};
  55     case 4
  56         rows = varargin{1};
  57         columns = varargin{2};
  58         layers =varargin{3};
  59 end
  60
  61 if usepm
  62     [hd,im] = pmedf_read(filename);
  63 else
  64     fid=fopen(filename,'r');
  65
  66     if fid==-1
  67             im=fid;
  68             disp(sprintf('!!! Error opening file %s !!!',filename))
  69     else
  70             hd=readedfheader(fid);
  71
  72         byteorder=findheader(hd,'ByteOrder','string');
  73         fidpos=ftell(fid); % store present position in file
  74         fclose(fid);
  75         if strcmp(byteorder,'HighByteFirst')
  76             byteorder='b';
  77         else
  78             byteorder='l';
  79         end
  80         fid=fopen(filename,'r',byteorder); % re-open with good format
  81         fseek(fid,fidpos,0); % re-position at end of header
  82
  83         xsize=findheader(hd,'Dim_1','integer');
  84             ysize=findheader(hd,'Dim_2','integer');
  85         zsize=findheader(hd,'Dim_3','integer');
  86         if isempty(zsize)
  87             zsize=1;
  88         end
  89
  90         datatype=findheader(hd,'DataType','string');
  91             switch datatype
  92                     case 'UnsignedByte',
  93                             datatype='uint8';
  94                             nbytes=1;
  95                     case 'UnsignedShort',
  96                             datatype='uint16';
  97                             nbytes=2;
  98                     case {'UnsignedInteger','UnsignedLong'}
  99                             datatype='uint32';
 100                             nbytes=4;
 101                     case {'Float','FloatValue','FLOATVALUE','Real'}
 102                             datatype='float32';
 103                             nbytes=4;
 104                     case 'DoubleValue'
 105                             datatype='float64';
 106                             nbytes=8;
 107     %etcetera
 108             end
 109
 110             if isempty(who('rows'))
 111                     rows=1:xsize;
 112             end
 113             if isempty(who('columns'))
 114                     columns=1:ysize;
 115         end
 116         if isempty(who('layers'))
 117                     layers=1:zsize;
 118         end
 119
 120         if zsize==1
 121             fseek(fid,nbytes*(rows(1)-1+(columns(1)-1)*xsize),0);
 122                 im=fread(fid,[length(rows),length(columns)],sprintf('%d*%s',length(rows),datatype),nbytes*(xsize-length(rows)));
 123         else
 124             j=1;
 125             for i=layers
 126                 fseek(fid,1024+nbytes*(rows(1)-1+(columns(1)-1)*xsize+xsize*ysize*(i-1)),-1);
 127                 im(:,:,j)=fread(fid,[length(rows),length(columns)],sprintf('%d*%s',length(rows),datatype),nbytes*(xsize-length(rows)));
 128                 j=j+1;
 129             end
 130         end
 131         st=fclose(fid);
 132
 133     end
 134 end % if usepm