useful functions for simulations, created by ESRF people mainly (free to use)

[CreaPhase.git] / utilities_ESRF / edfread.m

diff --git a/utilities_ESRF/edfread.m b/utilities_ESRF/edfread.m
new file mode 100644 (file)
index 0000000..75eca2a
--- /dev/null
+++ b/utilities_ESRF/edfread.m
@@ -0,0 +1,134 @@
+## Copyright (C) 2007 P. Cloetens
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+## edfread
+## function im=edfread(filename,varargin)
+##      reads an image in edf format into the matrix im
+##      uses pmedf_read by Petr Mikulik if only filename is given
+##      if both header and image are needed, use [hd,im]=pmedf_read(filename)
+##      does not transpose the image, 1st dimension remains 1st dimension of original file
+##      
+##      arguments:
+##      argument 1: filename
+##      argument 2 (optional): vector of adjacent rows to be read
+##      argument 3 (optional): vector of adjacent colums to be read
+##      argument 4 (optional): vector of layers to be read
+##      
+##      examples:
+##      im = edfread('dark.edf');
+##          reads complete image
+##      im = edfread('dark.edf',10:200,30:50);
+##          reads subimage corresponding to row (Dim_1) 10 until 200 and column (Dim_2) 30 until 50
+##      im = edfread('dark.edf',10:200,30:50,[1 4 5]);
+##          reads layers 1, 4 and 5 of a 3-dimensional image (see Dim_3 and W. Ludwig)
+##
+##      See also: edfwrite, pmedf_read, pmedf_write
+
+## Author: P. Cloetens <cloetens@esrf.fr>
+##
+## 2007-03-03 P. Cloetens <cloetens@esrf.fr>
+## * Initial revision
+
+function im=edfread(filename,varargin)
+
+usepm = 0;
+
+switch nargin
+    case 1
+        usepm = isoctave&exist('pmedf_read');
+    case 3
+        rows = varargin{1};
+        columns = varargin{2};
+    case 4
+        rows = varargin{1};
+        columns = varargin{2};
+        layers =varargin{3};  
+end
+
+if usepm
+    [hd,im] = pmedf_read(filename);
+else
+    fid=fopen(filename,'r');
+
+    if fid==-1
+           im=fid;
+           disp(sprintf('!!! Error opening file %s !!!',filename))
+    else
+           hd=readedfheader(fid);
+
+        byteorder=findheader(hd,'ByteOrder','string');
+        fidpos=ftell(fid); % store present position in file
+        fclose(fid); 
+        if strcmp(byteorder,'HighByteFirst')
+            byteorder='b';
+        else
+            byteorder='l';
+        end
+        fid=fopen(filename,'r',byteorder); % re-open with good format
+        fseek(fid,fidpos,0); % re-position at end of header
+
+        xsize=findheader(hd,'Dim_1','integer');
+           ysize=findheader(hd,'Dim_2','integer');
+        zsize=findheader(hd,'Dim_3','integer');
+        if isempty(zsize)
+            zsize=1;
+        end
+
+        datatype=findheader(hd,'DataType','string');
+           switch datatype
+                   case 'UnsignedByte',
+                           datatype='uint8';
+                           nbytes=1;
+                   case 'UnsignedShort',
+                           datatype='uint16';
+                           nbytes=2;
+                   case {'UnsignedInteger','UnsignedLong'}
+                           datatype='uint32';
+                           nbytes=4;
+                   case {'Float','FloatValue','FLOATVALUE','Real'}
+                           datatype='float32';
+                           nbytes=4;
+                   case 'DoubleValue'
+                           datatype='float64';
+                           nbytes=8;
+    %etcetera
+           end
+
+           if isempty(who('rows'))
+                   rows=1:xsize;
+           end
+           if isempty(who('columns'))
+                   columns=1:ysize;
+        end
+        if isempty(who('layers'))
+                   layers=1:zsize;
+        end
+
+        if zsize==1
+           fseek(fid,nbytes*(rows(1)-1+(columns(1)-1)*xsize),0);
+               im=fread(fid,[length(rows),length(columns)],sprintf('%d*%s',length(rows),datatype),nbytes*(xsize-length(rows)));
+        else
+            j=1;
+            for i=layers
+                fseek(fid,1024+nbytes*(rows(1)-1+(columns(1)-1)*xsize+xsize*ysize*(i-1)),-1);
+                im(:,:,j)=fread(fid,[length(rows),length(columns)],sprintf('%d*%s',length(rows),datatype),nbytes*(xsize-length(rows)));        
+                j=j+1;
+            end
+        end    
+        st=fclose(fid);
+
+    end
+end % if usepm