--- /dev/null
+## Copyright (C) 2012 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
+
+## cut
+## function imcut = cut(im, nn, mn[, cnn, cmn])
+## cuts a sub-matrix out of a larger matrix
+## cuts in the center of the original matrix, except if new center is specified
+## NO CHECKING of validity indices sub-matrix!
+##
+## arguments:
+## argument 1: original matrix im
+## argument 2: number of rows in result
+## argument 3: number of columns in result
+## argument 4: center row around which to cut ( default: center )
+## argument 5: center column around which to cut ( default: center )
+##
+## examples:
+## im_roi = cut(im, 1024, 1024) -> cut center 1024x1024 pixels
+## im_roi = cut(im, 1024, 1024, 600.5, 700.5) -> cut 1024x1024 pixels around pixels (600-601, 700-701)
+
+## Author: P. Cloetens <cloetens@esrf.eu>
+##
+## 2012-10-14 P. Cloetens <cloetens@esrf.eu>
+## * Initial revision
+
+function imcut = cut(im, nn, mn, cnn, cmn)
+ if !exist('cnn', 'var')
+ cnn = [];
+ endif
+ if !exist('cmn', 'var')
+ cmn = [];
+ endif
+
+ [n,m] = size(im);
+ if isempty(cnn)
+ cnn = (n+1)/2;
+ endif
+ if isempty(cmn)
+ cmn = (m+1)/2;
+ endif
+
+ rb = round(0.5+cnn-nn/2);
+ re = nn+rb-1;
+ cb = round(0.5+cmn-mn/2);
+ ce = mn+cb-1;
+ imcut = im(rb:re, cb:ce);
+endfunction
--- /dev/null
+## 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
--- /dev/null
+## function count=edfwrite(filename,matrix,datatype,head)
+## writes an image in esrf data format
+##
+## datatype can be
+## uint8 = 8 bits = UnsignedByte
+## uint16 = 16 bits = UnsignedShort
+## uint32 = 32 bits = UnsignedLong = UnsignedInteger
+## float32 = = Float = Real
+## float64 = = DoubleValue = Double
+## writes bigendian files
+##
+## head (optional) is a structure whose fields are written to the file
+##
+## 30.09.2009 (HSu): Modified to write header structures
+function count=edfwrite(filename,matrix,datatype,varargin)
+
+ switch nargin
+ case 3
+ head=[];
+ case 4
+ head=varargin{1};
+ end
+
+ switch datatype
+ case 'uint8',
+ esrfdatatype='UnsignedByte';
+ nbytes=1;
+ case 'uint16',
+ esrfdatatype='UnsignedShort';
+ nbytes=2;
+ case 'uint32',
+ esrfdatatype='UnsignedLong';
+ nbytes=4;
+ case 'float32',
+ esrfdatatype='Float';
+ nbytes=4;
+ case 'float64',
+ esrfdatatype='DoubleValue';
+ nbytes=8;
+ end
+
+ if isempty(head) || isstruct(head)
+ head=writeheader(size(matrix),esrfdatatype,nbytes,1,head);
+ end
+
+ fid=fopen(filename,'w','b');
+ fwrite(fid,head);
+ fwrite(fid,matrix,datatype);
+ fclose(fid);
+endfunction
--- /dev/null
+## Copyright (C) 2010 Peter 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 Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## findheader
+## function [res, pos] = findheader(hd, what, kind)
+## read a given info from header of edf-files (.edf), info-file (.info) or parameter file (.par)
+##
+## arguments
+## argument 1 : header or info string
+## argument 2 : parameter to be read (can be 'date')
+## argument 3 : kind|type of parameter
+## possible values are: 'motor', 'counter', 'integer', 'float', 'list' or 'string'
+
+## Author: Peter Cloetens cloetens@esrf.eu
+## Created: 2010-12-10
+## * main change from historical version: avoid any assumption on offset, change help
+
+function [res, pos] = findheader(hd, what, kind)
+
+ if nargin < 3
+ disp('3 input arguments required:')
+ help findheader
+ return
+ end
+
+ if strcmp(what,'date')
+ pos = findstr(hd,what);
+ if !isempty(pos)
+ pos_equal = findstr(hd(pos:end),'=');
+ pos_equal = pos(1)+pos_equal(1)-1;
+ pos_semicolon = findstr(hd(pos_equal:end),';');
+ pos_semicolon = pos_equal+pos_semicolon(1)-1;
+ res = hd(pos_equal+2:pos_semicolon-2);
+ else % what is not in header
+ res = [];
+ end
+ elseif strcmp(kind,'motor')
+ pos = findstr(hd,'motor_mne');
+ if !isempty(pos)
+ pos_equal = findstr(hd(pos:end),'=');
+ pos_equal = pos(1)+pos_equal(1)-1;
+ pos_semicolon = findstr(hd(pos_equal:end),';');
+ pos_semicolon = pos_equal+pos_semicolon(1)-1;
+ pos = pos_equal+1;
+ posend = pos_semicolon-1;
+ motnum = 1;
+ res = []; % in case motor does not exist
+ while pos<posend; % looking number of motor
+ motor = sscanf(hd(pos:posend),'%s',1);
+ if strcmp(motor, what) % motor was found -> taking corresponding position
+ pos = findstr(hd,'motor_pos');
+ pos_equal = findstr(hd(pos:end),'=');
+ pos_equal = pos(1)+pos_equal(1)-1;
+ res = sscanf(hd(pos_equal+1:end),'%g',motnum);
+ res = res(motnum);
+ break
+ endif
+ motnum++;
+ pos += findstr(hd(pos:posend), motor)(1) + length(motor);
+ endwhile
+ else
+ res = [];
+ endif
+ elseif strcmp(kind,'counter')
+ pos=findstr(hd,'counter_mne');
+ if not(isempty(pos))
+ pos_equal = findstr(hd(pos:end),'=');
+ pos_equal = pos(1)+pos_equal(1)-1;
+ pos_semicolon = findstr(hd(pos_equal:end),';');
+ pos_semicolon = pos_equal+pos_semicolon(1)-1;
+ pos = pos_equal+1;
+ posend = pos_semicolon-1;
+ motnum = 1;
+ res = []; % in case counter does not exist
+ while pos<posend; % looking number of motor
+ motor = sscanf(hd(pos:posend),'%s',1);
+ if strcmp(motor,what) % motor was found -> taking corresponding position
+ pos = findstr(hd,'counter_pos');
+ pos_equal = findstr(hd(pos:end),'=');
+ pos_equal = pos(1)+pos_equal(1)-1;
+ res = sscanf(hd(pos_equal+1:end),'%g',motnum);
+ res = res(motnum);
+ break
+ endif
+ motnum++;
+ pos += findstr(hd(pos:posend), motor)(1) + length(motor);
+ endwhile
+ else
+ res = [];
+ endif
+ else # kind is not 'motor' 'counter' 'date'
+ pos=findstr(hd,what);
+ if not(isempty(pos))
+ pos = pos(1);
+ pos_equal = findstr(hd(pos:end),'=');
+ pos_equal = pos+pos_equal(1);
+ switch kind
+ case 'integer',
+ res = sscanf(hd(pos_equal:end),'%i');
+ case 'float',
+ res = sscanf(hd(pos_equal:end),'%f');
+ case 'string',
+ res = sscanf(hd(pos_equal:end),'%s',1);
+ case 'list',
+ pos_end = findstr(hd(pos_equal:end),'}');
+ res = hd(pos_equal:pos_equal+pos_end-1);
+ otherwise,
+ res = [];
+ endswitch
+ else
+ res = [];
+ endif
+ endif
+endfunction
--- /dev/null
+## function str=headerstring(what,value,kind)
+## used by writeheader
+## returns a string for incorporation in an edf-header
+## what is the information item (e.g. 'Dim_1')
+## value of the item
+## kind is string, char, integer, logical, float or double
+## origin: peter
+##
+## 30.09.2009 (HSu): Updated to work with double and char data types
+## 2011-08-29 PC
+## * add logical data type
+
+function str=headerstring(what,value,kind)
+
+ str=[what ' = '];
+
+ switch kind
+ case {'string', 'char'},
+ str=[str sprintf('%s ',value')];
+ case {'integer', 'logical'},
+ str=[str sprintf('%i ',value')];
+ case {'float', 'double'}
+ str=[str sprintf('%g ',value')];
+ otherwise
+ printf("Unknown data type in header %s\n",kind);
+ end
+
+ str = [str sprintf(';\n')];
+end
--- /dev/null
+## Copyright (C) 2000 Teemu Ikonen
+##
+## 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, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} impad(@var{A}, @var{xpad}, @var{ypad}, [@var{padding}, [@var{const}]])
+## Pad (augment) a matrix for application of image processing algorithms.
+##
+## Pads the input image @var{A} with @var{xpad}(1) elements from left,
+## @var{xpad}(2), elements from right, @var{ypad}(1) elements from above
+## and @var{ypad}(2) elements from below.
+## Values of padding elements are determined from the optional arguments
+## @var{padding} and @var{const}. @var{padding} is one of
+##
+## @table @samp
+## @item "zeros"
+## pad with zeros (default)
+##
+## @item "ones"
+## pad with ones
+##
+## @item "constant"
+## pad with a value obtained from the optional fifth argument const
+##
+## @item "symmetric"
+## pad with values obtained from @var{A} so that the padded image mirrors
+## @var{A} starting from edges of @var{A}
+##
+## @item "reflect"
+## same as symmetric, but the edge rows and columns are not used in the padding
+##
+## @item "replicate"
+## pad with values obtained from A so that the padded image
+## repeates itself in two dimensions
+##
+## @end table
+## @end deftypefn
+
+## Author: Teemu Ikonen <tpikonen@pcu.helsinki.fi>
+## Created: 5.5.2000
+## Keywords: padding image processing
+
+## A nice test matrix for padding:
+## A = 10*[1:5]' * ones(1,5) + ones(5,1)*[1:5]
+
+function retim = impad(im, xpad, ypad, padding = "zeros", const = 1)
+ ## Input checking
+ if (!ismatrix(im) || ndims(im) < 2 || ndims(im) > 3)
+ error("impad: first input argument must be an image");
+ endif
+ if (isscalar(xpad) && xpad >= 0)
+ xpad(2) = xpad;
+ elseif (!isreal(xpad) || numel(xpad) != 2 || any(xpad < 0))
+ error("impad: xpad must be a positive scalar or 2-vector");
+ endif
+ if (isscalar(ypad) && ypad >= 0)
+ ypad(2) = ypad;
+ elseif (!isreal(ypad) || numel(ypad) != 2 || any(ypad < 0))
+ error("impad: ypad must be a positive scalar or 2-vector");
+ endif
+ if (!isscalar(const))
+ error("impad: fifth input argument must be a scalar");
+ endif
+
+ origx = size(im,2);
+ origy = size(im,1);
+ retx = origx + xpad(1) + xpad(2);
+ rety = origy + ypad(1) + ypad(2);
+ cl = class(im);
+
+ for iter = size(im,3):-1:1
+ A = im(:,:,iter);
+ switch (lower(padding))
+ case "zeros"
+ retval = zeros(rety, retx, cl);
+ retval(ypad(1)+1 : ypad(1)+origy, xpad(1)+1 : xpad(1)+origx) = A;
+ case "ones"
+ retval = ones(rety, retx, cl);
+ retval(ypad(1)+1 : ypad(1)+origy, xpad(1)+1 : xpad(1)+origx) = A;
+ case "constant"
+ retval = const.*ones(rety, retx, cl);
+ retval(ypad(1)+1 : ypad(1)+origy, xpad(1)+1 : xpad(1)+origx) = A;
+ case "replicate"
+ y1 = origy-ypad(1)+1;
+ x1 = origx-xpad(1)+1;
+ if (y1 < 1 || x1 < 1 || ypad(2) > origy || xpad(2) > origx)
+ error("impad: too large padding for this padding type");
+ else
+ yrange1 = y1:origy;
+ yrange2 = 1:ypad(2);
+ xrange1 = x1:origx;
+ xrange2 = 1:xpad(2);
+ retval = [ A(yrange1, xrange1), A(yrange1, :), A(yrange1, xrange2);
+ A(:, xrange1), A, A(:, xrange2);
+ A(yrange2, xrange1), A(yrange2, :), A(yrange2, xrange2) ];
+ endif
+ case "symmetric"
+ y2 = origy-ypad(2)+1;
+ x2 = origx-xpad(2)+1;
+ if (ypad(1) > origy || xpad(1) > origx || y2 < 1 || x2 < 1)
+ error("impad: too large padding for this padding type");
+ else
+ yrange1 = 1 : ypad(1);
+ yrange2 = y2 : origy;
+ xrange1 = 1 : xpad(1);
+ xrange2 = x2 : origx;
+ retval = [ fliplr(flipud(A(yrange1, xrange1))), flipud(A(yrange1, :)), fliplr(flipud(A(yrange1, xrange2)));
+ fliplr(A(:, xrange1)), A, fliplr(A(:, xrange2));
+ fliplr(flipud(A(yrange2, xrange1))), flipud(A(yrange2, :)), fliplr(flipud(A(yrange2, xrange2))) ];
+ endif
+ case "reflect"
+ y2 = origy-ypad(2);
+ x2 = origx-xpad(2);
+ if (ypad(1)+1 > origy || xpad(1)+1 > origx || y2 < 1 || x2 < 1)
+ error("impad: too large padding for this padding type");
+ else
+ yrange1 = 2 : ypad(1)+1;
+ yrange2 = y2 : origy-1;
+ xrange1 = 2 : xpad(1)+1;
+ xrange2 = x2 : origx-1;
+ retval = [ fliplr(flipud(A(yrange1, xrange1))), flipud(A(yrange1, :)), fliplr(flipud(A(yrange1, xrange2)));
+ fliplr(A(:, xrange1)), A, fliplr(A(:, xrange2));
+ fliplr(flipud(A(yrange2, xrange1))), flipud(A(yrange2, :)), fliplr(flipud(A(yrange2, xrange2))) ];
+ endif
+ otherwise
+ error("impad: unknown padding type");
+ endswitch
+
+ retim(:,:,iter) = retval;
+ endfor
+endfunction
--- /dev/null
+% function [ ret ] = isoctave ()
+% equal to 1 when executed from octave, to 0 when executed from non-octave
+
+% Author: P. Cloetens <cloetens@esrf.fr>
+%
+% 2006-06-04 P. Cloetens <cloetens@esrf.fr>
+% Initial revision
+
+function [ ret ] = isoctave ()
+ if exist('OCTAVE_VERSION','builtin')
+ ret = 1;
+ else
+ ret = 0;
+ end
+end
--- /dev/null
+% function hd=readedfheader(fid)
+
+function hd=readedfheader(fid)
+
+headerlength = 512;
+closing = '}';
+hd = fscanf(fid,'%c',headerlength);
+
+while not(strcmp(hd(end-1),closing))
+ hd = [hd fscanf(fid,'%c',headerlength)];
+end
--- /dev/null
+## function header=writeheader(sizes,datatype,nbytes,nimage)
+## used by edfwrite
+## sizes is size(image to be written), cf. Dim_1 and Dim_2
+## datatype is string, cf. DataType
+## nbytes is number of bytes used for each written element
+## nimage is probably number of the image, used for writing multiple images in a single file?
+## origin: peter
+## 2011-05-18 PC
+## * correct error in number of blanks to add to get multiple of 512 bytes header
+## * avoid sprintf, use double quotes
+
+function header = writeheader(sizes,datatype,nbytes,nimage,headstruct)
+
+ header = sprintf("%s\n","{");
+ header = [header headerstring("HeaderID","EH:000001:000000:000000","string")];
+ header = [header headerstring("Image",nimage,"integer")];
+ header = [header headerstring("ByteOrder","HighByteFirst","string")];
+ header = [header headerstring("DataType",datatype,"string")];
+ header = [header headerstring("Dim_1",sizes(1),"integer")];
+ header = [header headerstring("Dim_2",sizes(2),"integer")];
+ header = [header headerstring("Size",nbytes*sizes(1)*sizes(2),"integer")];
+ header = [header headerstring("Date",date,"string")];
+
+ if ((nargin == 5) && (isstruct(headstruct) == 1))
+ ## If the structure exists, write the fields to the file
+ if (!isempty(headstruct))
+ names = fieldnames(headstruct);
+ for ii=1:length(names),
+ header = [header sprintf("%s",headerstring(names{ii},getfield(headstruct,names{ii}),...
+ class(getfield(headstruct,names{ii}))))];
+ endfor
+ endif
+ endif
+
+ ## Fill to next multiple of 512
+ lnspace = mod(length(header)+2, 512);
+ if lnspace > 0
+ lnspace = 512-lnspace;
+ endif
+ header = [header blanks(lnspace) "}\n"];
+endfunction
+
+
+
git://git.creatis.insa-lyon.fr / CreaPhase.git / commitdiff