From: Loriane Weber Date: Wed, 20 Apr 2016 08:51:59 +0000 (+0200) Subject: useful functions for simulations, created by ESRF people mainly (free to use) X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=commitdiff_plain;h=99b9890c11d31f9ae22dd481873e01c708175073 useful functions for simulations, created by ESRF people mainly (free to use) --- diff --git a/utilities_ESRF/cut.m b/utilities_ESRF/cut.m new file mode 100644 index 0000000..007e8e8 --- /dev/null +++ b/utilities_ESRF/cut.m @@ -0,0 +1,60 @@ +## 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 +## +## 2012-10-14 P. Cloetens +## * 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 diff --git a/utilities_ESRF/edfread.m b/utilities_ESRF/edfread.m new file mode 100644 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 +## +## 2007-03-03 P. Cloetens +## * 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 diff --git a/utilities_ESRF/edfwrite.m b/utilities_ESRF/edfwrite.m new file mode 100644 index 0000000..d69c77d --- /dev/null +++ b/utilities_ESRF/edfwrite.m @@ -0,0 +1,50 @@ +## 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 diff --git a/utilities_ESRF/findheader.m b/utilities_ESRF/findheader.m new file mode 100644 index 0000000..a24aca1 --- /dev/null +++ b/utilities_ESRF/findheader.m @@ -0,0 +1,127 @@ +## 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 +## . + +## 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 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 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 diff --git a/utilities_ESRF/headerstring.m b/utilities_ESRF/headerstring.m new file mode 100644 index 0000000..a7d7083 --- /dev/null +++ b/utilities_ESRF/headerstring.m @@ -0,0 +1,29 @@ +## 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 diff --git a/utilities_ESRF/impad.m b/utilities_ESRF/impad.m new file mode 100644 index 0000000..d8151f9 --- /dev/null +++ b/utilities_ESRF/impad.m @@ -0,0 +1,142 @@ +## 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 . + +## -*- 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 +## 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 diff --git a/utilities_ESRF/isoctave.m b/utilities_ESRF/isoctave.m new file mode 100644 index 0000000..e2110f8 --- /dev/null +++ b/utilities_ESRF/isoctave.m @@ -0,0 +1,15 @@ +% function [ ret ] = isoctave () +% equal to 1 when executed from octave, to 0 when executed from non-octave + +% Author: P. Cloetens +% +% 2006-06-04 P. Cloetens +% Initial revision + +function [ ret ] = isoctave () + if exist('OCTAVE_VERSION','builtin') + ret = 1; + else + ret = 0; + end +end diff --git a/utilities_ESRF/readedfheader.m b/utilities_ESRF/readedfheader.m new file mode 100644 index 0000000..132f428 --- /dev/null +++ b/utilities_ESRF/readedfheader.m @@ -0,0 +1,11 @@ +% 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 diff --git a/utilities_ESRF/writeheader.m b/utilities_ESRF/writeheader.m new file mode 100644 index 0000000..7a4ebcf --- /dev/null +++ b/utilities_ESRF/writeheader.m @@ -0,0 +1,44 @@ +## 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 + + +