text:
function varargout = bbGt( action, varargin )
% Bounding box (bb) annotations struct, evaluation and sampling routines.
%
% bbGt gives access to two types of routines:
% (1) Data structure for storing bb image annotations.
% (2) Routines for evaluating the Pascal criteria for object detection.
%
% The bb annotation stores bb for objects of interest with additional
% information per object, such as occlusion information. The underlying
% data structure is simply a Matlab stuct array, one struct per object.
% This annotation format is an alternative to the annotation format used
% for the PASCAL object challenges (in addition routines for loading PASCAL
% format data are provided, see bbLoad()).
%
% Each object struct has the following fields:
%  lbl  - a string label describing object type (eg: 'pedestrian')
%  bb   - [l t w h]: bb indicating predicted object extent
%  occ  - 0/1 value indicating if bb is occluded
%  bbv  - [l t w h]: bb indicating visible region (may be [0 0 0 0])
%  ign  - 0/1 value indicating bb was marked as ignore
%  ang  - [0-360] orientation of bb in degrees
%
% Note: although orientation (angle) is stored for each bb, for now it is
% not being used during evaluation or sampling.
%
% bbGt contains a number of utility functions, accessed using:
%  outputs = bbGt( 'action', inputs );
% The list of functions and help for each is given below. Also, help on
% individual subfunctions can be accessed by: "help bbGt>action".
%
%%% (1) Data structure for storing bb image annotations.
% Create annotation of n empty objects.
%   objs = bbGt( 'create', [n] );
% Save bb annotation to text file.
%   objs = bbGt( 'bbSave', objs, fName )
% Load bb annotation from text file and filter.
%   [objs,bbs] = bbGt( 'bbLoad', fName, [pLoad] )
% Get object property 'name' (in a standard array).
%   vals = bbGt( 'get', objs, name )
% Set object property 'name' (with a standard array).
%   objs = bbGt( 'set', objs, name, vals )
% Draw an ellipse for each labeled object.
%   hs = draw( objs, pDraw )
%
%%% (2) Routines for evaluating the Pascal criteria for object detection.
% Get all corresponding files in given directories.
%   [fs,fs0] = bbGt('getFiles', dirs, [f0], [f1] )
% Copy corresponding files into given directories.
%   fs = bbGt( 'copyFiles', fs, dirs )
% Load all ground truth and detection bbs in given directories.
%   [gt0,dt0] = bbGt( 'loadAll', gtDir, [dtDir], [pLoad] )
% Evaluates detections against ground truth data.
%   [gt,dt] = bbGt( 'evalRes', gt0, dt0, [thr], [mul] )
% Display evaluation results for given image.
%   [hs,hImg] = bbGt( 'showRes' I, gt, dt, varargin )
% Compute ROC or PR based on outputs of evalRes on multiple images.
%   [xs,ys,ref] = bbGt( 'compRoc', gt, dt, roc, ref )
% Extract true or false positives or negatives for visualization.
%   [Is,scores,imgIds] = bbGt( 'cropRes', gt, dt, imFs, varargin )
% Computes (modified) overlap area between pairs of bbs.
%   oa = bbGt( 'compOas', dt, gt, [ig] )
% Optimized version of compOas for a single pair of bbs.
%   oa = bbGt( 'compOa', dt, gt, ig )
%
% USAGE
%  varargout = bbGt( action, varargin );
%
% INPUTS
%  action     - string specifying action
%  varargin   - depends on action, see above
%
% OUTPUTS
%  varargout  - depends on action, see above
%
% EXAMPLE
%
% See also bbApply, bbLabeler, bbGt>create, bbGt>bbSave, bbGt>bbLoad,
% bbGt>get, bbGt>set, bbGt>draw, bbGt>getFiles, bbGt>copyFiles,
% bbGt>loadAll, bbGt>evalRes, bbGt>showRes,  bbGt>compRoc, bbGt>cropRes,
% bbGt>compOas, bbGt>compOa
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.26
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

%#ok<*DEFNU>
varargout = cell(1,max(1,nargout));
[varargout{:}] = feval(action,varargin{:});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function objs = create( n )
% Create annotation of n empty objects.
%
% USAGE
%  objs = bbGt( 'create', [n] )
%
% INPUTS
%  n      - [1] number of objects to create
%
% OUTPUTS
%  objs   - annotation of n 'empty' objects
%
% EXAMPLE
%  objs = bbGt('create')
%
% See also bbGt
o=struct('lbl','','bb',[0 0 0 0],'occ',0,'bbv',[0 0 0 0],'ign',0,'ang',0);
if(nargin<1 || n==1), objs=o; return; end; objs=o(ones(n,1));
end

function objs = bbSave( objs, fName )
% Save bb annotation to text file.
%
% USAGE
%  objs = bbGt( 'bbSave', objs, fName )
%
% INPUTS
%  objs   - objects to save
%  fName  - name of text file
%
% OUTPUTS
%  objs   - objects to save
%
% EXAMPLE
%
% See also bbGt, bbGt>bbLoad
vers=3; fid=fopen(fName,'w'); assert(fid>0);
fprintf(fid,'%% bbGt version=%i\n',vers);
objs=set(objs,'bb',round(get(objs,'bb')));
objs=set(objs,'bbv',round(get(objs,'bbv')));
objs=set(objs,'ang',round(get(objs,'ang')));
for i=1:length(objs)
  o=objs(i); bb=o.bb; bbv=o.bbv;
  fprintf(fid,['%s' repmat(' %i',1,11) '\n'],o.lbl,...
    bb,o.occ,bbv,o.ign,o.ang);
end
fclose(fid);
end

function [objs,bbs] = bbLoad( fName, varargin )
% Load bb annotation from text file and filter.
%
% FORMAT: Specify 'format' to indicate the format of the ground truth.
% format=0 is the default format (created by bbSave/bbLabeler). format=1 is
% the PASCAL VOC format. Loading ground truth in this format requires
% 'VOCcode/' to be in directory path. It's part of VOCdevkit available from
% the PASCAL VOC: http://pascallin.ecs.soton.ac.uk/challenges/VOC/. Objects
% labeled as either 'truncated' or 'occluded' using the PASCAL definitions
% have the 'occ' flag set to true. Objects labeled as 'difficult' have the
% 'ign' flag set to true. 'class' is used for 'lbl'. format=2 is the
% ImageNet detection format and requires the ImageNet Dev Kit.
%
% FILTERING: After loading, the objects can be filtered. First, only
% objects with lbl in lbls or ilbls or returned. For each object, obj.ign
% is set to 1 if it was already at 1, if its label was in ilbls, or if any
% object property is outside of the specified range. The ignore flag is
% used during training and testing so that objects with certain properties
% (such as very small or heavily occluded objects) are excluded. The range
% for each property is a two element vector, [0 inf] by default; a property
% value v is inside the range if v>=rng(1) && v<=rng(2). Tested properties
% include height (h), width (w), area (a), aspect ratio (ar), orientation
% (o), extent x-coordinate (x), extent y-coordinate (y), and fraction
% visible (v). The last property is computed as the visible object area
% divided by the total area, except if o.occ==0, in which case v=1, or
% all(o.bbv==o.bb), which indicates the object may be barely visible, in
% which case v=0 (note that v~=1 in this case).
%
% RETURN: In addition to outputting the objs, bbLoad() can return the
% corresponding bounding boxes (bbs) in an [nx5] array where each row is of
% the form [x y w h ignore], [x y w h] is the bb and ignore=obj.ign. For
% oriented bbs, the extent of the bb is returned, where the extent is the
% smallest axis aligned bb containing the oriented bb. If the oriented bb
% was labeled as a rectangle as opposed to an ellipse, the tightest bb will
% usually increase slightly in size due to the corners of the rectangle
% sticking out beyond the ellipse bounds. The 'ellipse' flag controls how
% an oriented bb is converted to a regular bb. Specifically, set ellipse=1
% if an ellipse tightly delineates the object and 0 if a rectangle does.
% Finally, if 'squarify' is not empty the (non-ignore) bbs are converted to
% a fixed aspect ratio using bbs=bbApply('squarify',bbs,squarify{:}).
%
% USAGE
%  [objs,bbs] = bbGt( 'bbLoad', fName, [pLoad] )
%
% INPUTS
%  fName    - name of text file
%  pLoad    - parameters (struct or name/value pairs)
%   .format   - [0] gt format 0:default, 1:PASCAL, 2:ImageNet
%   .ellipse  - [1] controls how oriented bb is converted to regular bb
%   .squarify - [] controls optional reshaping of bbs to fixed aspect ratio
%   .lbls     - [] return objs with these labels (or [] to return all)
%   .ilbls    - [] return objs with these labels but set to ignore
%   .hRng     - [] range of acceptable obj heights
%   .wRng     - [] range of acceptable obj widths
%   .aRng     - [] range of acceptable obj areas
%   .arRng    - [] range of acceptable obj aspect ratios
%   .oRng     - [] range of acceptable obj orientations (angles)
%   .xRng     - [] range of x coordinates of bb extent
%   .yRng     - [] range of y coordinates of bb extent
%   .vRng     - [] range of acceptable obj occlusion levels
%
% OUTPUTS
%  objs     - loaded objects
%  bbs      - [nx5] array containg ground truth bbs [x y w h ignore]
%
% EXAMPLE
%
% See also bbGt, bbGt>bbSave

% get parameters
df={'format',0,'ellipse',1,'squarify',[],'lbls',[],'ilbls',[],'hRng',[],...
  'wRng',[],'aRng',[],'arRng',[],'oRng',[],'xRng',[],'yRng',[],'vRng',[]};
[format,ellipse,sqr,lbls,ilbls,hRng,wRng,aRng,arRng,oRng,xRng,yRng,vRng]...
  = getPrmDflt(varargin,df,1);

% load objs
if( format==0 )
  % load objs stored in default format
  fId=fopen(fName);
  if(fId==-1), error(['unable to open file: ' fName]); end; v=0;
  try v=textscan(fId,'%% bbGt version=%d'); v=v{1}; catch, end %#ok<CTCH>
  if(isempty(v)), v=0; end
  % read in annotation (m is number of fields for given version v)
  if(all(v~=[0 1 2 3])), error('Unknown version %i.',v); end
  frmt='%s %d %d %d %d %d %d %d %d %d %d %d';
  ms=[10 10 11 12]; m=ms(v+1); frmt=frmt(1:2+(m-1)*3);
  in=textscan(fId,frmt); for i=2:m, in{i}=double(in{i}); end; fclose(fId);
  % create objs struct from read in fields
  n=length(in{1}); objs=create(n);
  for i=1:n, objs(i).lbl=in{1}{i}; objs(i).occ=in{6}(i); end
  bb=[in{2} in{3} in{4} in{5}]; bbv=[in{7} in{8} in{9} in{10}];
  for i=1:n, objs(i).bb=bb(i,:); objs(i).bbv=bbv(i,:); end
  if(m>=11), for i=1:n, objs(i).ign=in{11}(i); end; end
  if(m>=12), for i=1:n, objs(i).ang=in{12}(i); end; end
elseif( format==1 )
  % load objs stored in PASCAL VOC format
  if(exist('PASreadrecord.m','file')~=2)
    error('bbLoad() requires the PASCAL VOC code.'); end
  os=PASreadrecord(fName); os=os.objects;
  n=length(os); objs=create(n);
  if(~isfield(os,'occluded')), for i=1:n, os(i).occluded=0; end; end
  for i=1:n
    bb=os(i).bbox; bb(3)=bb(3)-bb(1); bb(4)=bb(4)-bb(2); objs(i).bb=bb;
    objs(i).lbl=os(i).class; objs(i).ign=os(i).difficult;
    objs(i).occ=os(i).occluded || os(i).truncated;
    if(objs(i).occ), objs(i).bbv=bb; end
  end
elseif( format==2 )
  if(exist('VOCreadxml.m','file')~=2)
    error('bbLoad() requires the ImageNet dev code.'); end
  os=VOCreadxml(fName); os=os.annotation;
  if(isfield(os,'object')), os=os.object; else os=[]; end
  n=length(os); objs=create(n);
  for i=1:n
    bb=os(i).bndbox; bb=str2double({bb.xmin bb.ymin bb.xmax bb.ymax});
    bb(3)=bb(3)-bb(1); bb(4)=bb(4)-bb(2); objs(i).bb=bb;
    objs(i).lbl=os(i).name;
  end
else error('bbLoad() unknown format: %i',format);
end

% only keep objects whose lbl is in lbls or ilbls
if(~isempty(lbls) || ~isempty(ilbls)), K=true(n,1);
  for i=1:n, K(i)=any(strcmp(objs(i).lbl,[lbls ilbls])); end
  objs=objs(K); n=length(objs);
end

% filter objs (set ignore flags)
for i=1:n, objs(i).ang=mod(objs(i).ang,360); end
if(~isempty(ilbls)), for i=1:n, v=objs(i).lbl;
    objs(i).ign = objs(i).ign || any(strcmp(v,ilbls)); end; end
if(~isempty(xRng)),  for i=1:n, v=objs(i).bb(1);
    objs(i).ign = objs(i).ign || v<xRng(1) || v>xRng(2); end; end
if(~isempty(xRng)),  for i=1:n, v=objs(i).bb(1)+objs(i).bb(3);
    objs(i).ign = objs(i).ign || v<xRng(1) || v>xRng(2); end; end
if(~isempty(yRng)),  for i=1:n, v=objs(i).bb(2);
    objs(i).ign = objs(i).ign || v<yRng(1) || v>yRng(2); end; end
if(~isempty(yRng)),  for i=1:n, v=objs(i).bb(2)+objs(i).bb(4);
    objs(i).ign = objs(i).ign || v<yRng(1) || v>yRng(2); end; end
if(~isempty(wRng)),  for i=1:n, v=objs(i).bb(3);
    objs(i).ign = objs(i).ign || v<wRng(1) || v>wRng(2); end; end
if(~isempty(hRng)),  for i=1:n, v=objs(i).bb(4);
    objs(i).ign = objs(i).ign || v<hRng(1) || v>hRng(2); end; end
if(~isempty(oRng)),  for i=1:n, v=objs(i).ang; if(v>180), v=v-360; end
    objs(i).ign = objs(i).ign || v<oRng(1) || v>oRng(2); end; end
if(~isempty(aRng)),  for i=1:n, v=objs(i).bb(3)*objs(i).bb(4);
    objs(i).ign = objs(i).ign || v<aRng(1) || v>aRng(2); end; end
if(~isempty(arRng)), for i=1:n, v=objs(i).bb(3)/objs(i).bb(4);
    objs(i).ign = objs(i).ign || v<arRng(1) || v>arRng(2); end; end
if(~isempty(vRng)),  for i=1:n, o=objs(i); bb=o.bb; bbv=o.bbv; %#ok<ALIGN>
    if(~o.occ || all(bbv==0)), v=1; elseif(all(bbv==bb)), v=0; else
      v=(bbv(3)*bbv(4))/(bb(3)*bb(4)); end
    objs(i).ign = objs(i).ign || v<vRng(1) || v>vRng(2); end
end

% finally get extent of each bounding box (not trivial if ang~=0)
if(nargout<=1), return; end; if(n==0), bbs=zeros(0,5); return; end
bbs=double([reshape([objs.bb],4,[]); [objs.ign]]'); ign=bbs(:,5)==1;
for i=1:n, bbs(i,1:4)=bbExtent(bbs(i,1:4),objs(i).ang,ellipse); end
if(~isempty(sqr)), bbs(~ign,:)=bbApply('squarify',bbs(~ign,:),sqr{:}); end

  function bb = bbExtent( bb, ang, ellipse )
    % get bb that fully contains given oriented bb
    if(~ang), return; end
    if( ellipse ) % get bb that encompases ellipse (tighter)
      x=bbApply('getCenter',bb); a=bb(4)/2; b=bb(3)/2; ang=ang-90;
      rx=(a*cosd(ang))^2+(b*sind(ang))^2; rx=abs(rx/sqrt(rx));
      ry=(a*sind(ang))^2+(b*cosd(ang))^2; ry=abs(ry/sqrt(ry));
      bb=[x(1)-rx x(2)-ry 2*rx 2*ry];
    else % get bb that encompases rectangle (looser)
      c=cosd(ang); s=sind(ang); R=[c -s; s c]; rs=bb(3:4)/2;
      x0=-rs(1); x1=rs(1); y0=-rs(2); y1=rs(2); pc=bb(1:2)+rs;
      p=[x0 y0; x1 y0; x1 y1; x0 y1]*R'+pc(ones(4,1),:);
      x0=min(p(:,1)); x1=max(p(:,1)); y0=min(p(:,2)); y1=max(p(:,2));
      bb=[x0 y0 x1-x0 y1-y0];
    end
  end
end

function vals = get( objs, name )
% Get object property 'name' (in a standard array).
%
% USAGE
%  vals = bbGt( 'get', objs, name )
%
% INPUTS
%  objs   - [nx1] struct array of objects
%  name   - property name ('lbl','bb','occ',etc.)
%
% OUTPUTS
%  vals   - [nxk] array of n values (k=1 or 4)
%
% EXAMPLE
%
% See also bbGt, bbGt>set
nObj=length(objs); if(nObj==0), vals=[]; return; end
switch name
  case 'lbl', vals={objs.lbl}';
  case 'bb',  vals=reshape([objs.bb]',4,[])';
  case 'occ', vals=[objs.occ]';
  case 'bbv', vals=reshape([objs.bbv]',4,[])';
  case 'ign', vals=[objs.ign]';
  case 'ang', vals=[objs.ang]';
  otherwise, error('unkown type %s',name);
end
end

function objs = set( objs, name, vals )
% Set object property 'name' (with a standard array).
%
% USAGE
%  objs = bbGt( 'set', objs, name, vals )
%
% INPUTS
%  objs   - [nx1] struct array of objects
%  name   - property name ('lbl','bb','occ',etc.)
%  vals   - [nxk] array of n values (k=1 or 4)
%
% OUTPUTS
%  objs   - [nx1] struct array of updated objects
%
% EXAMPLE
%
% See also bbGt, bbGt>get
nObj=length(objs);
switch name
  case 'lbl', for i=1:nObj, objs(i).lbl=vals{i}; end
  case 'bb',  for i=1:nObj, objs(i).bb=vals(i,:); end
  case 'occ', for i=1:nObj, objs(i).occ=vals(i); end
  case 'bbv', for i=1:nObj, objs(i).bbv=vals(i,:); end
  case 'ign', for i=1:nObj, objs(i).ign=vals(i); end
  case 'ang', for i=1:nObj, objs(i).ang=vals(i); end
  otherwise, error('unkown type %s',name);
end
end

function hs = draw( objs, varargin )
% Draw an ellipse for each labeled object.
%
% USAGE
%  hs = bbGt( 'draw', objs, pDraw )
%
% INPUTS
%  objs       - [nx1] struct array of objects
%  pDraw      - parameters (struct or name/value pairs)
%   .col        - ['g'] color or [nx1] array of colors
%   .lw         - [2] line width
%   .ls         - ['-'] line style
%
% OUTPUTS
%  hs     - [nx1] handles to drawn graphic objects
%
% EXAMPLE
%
% See also bbGt
dfs={'col',[],'lw',2,'ls','-'};
[col,lw,ls]=getPrmDflt(varargin,dfs,1);
n=length(objs); hold on; hs=zeros(n,4);
if(isempty(col)), if(n==1), col='g'; else col=hsv(n); end; end
tProp={'FontSize',10,'color','w','FontWeight','bold',...
  'VerticalAlignment','bottom'};
for i=1:n
  bb=objs(i).bb; ci=col(i,:);
  hs(i,1)=text(bb(1),bb(2),objs(i).lbl,tProp{:});
  x=bbApply('getCenter',bb); r=bb(3:4)/2; a=objs(i).ang/180*pi-pi/2;
  [hs(i,2),hs(i,3),hs(i,4)]=plotEllipse(x(2),x(1),r(2),r(1),a,ci,[],lw,ls);
end; hold off;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fs,fs0] = getFiles( dirs, f0, f1 )
% Get all corresponding files in given directories.
%
% The first dir in 'dirs' serves as the baseline dir. getFiles() returns
% all files in the baseline dir and all corresponding files in the
% remaining dirs to the files in the baseline dir, in the same order. Two
% files are in correspondence if they have the same base name (regardless
% of extension). For example, given a file named "name.jpg", a
% corresponding file may be named "name.txt" or "name.jpg.txt". Every file
% in the baseline dir must have a matching file in the remaining dirs.
%
% USAGE
%  [fs,fs0] = bbGt('getFiles', dirs, [f0], [f1] )
%
% INPUTS
%   dirs      - {1xm} list of m directories
%   f0        - [1] index of first file in baseline dir to use
%   f1        - [inf] index of last file in baseline dir to use
%
% OUTPUTS
%   fs        - {mxn} list of full file names in each dir
%   fs0       - {1xn} list of file names without path or extensions
%
% EXAMPLE
%
% See also bbGt

if(nargin<2 || isempty(f0)), f0=1; end
if(nargin<3 || isempty(f1)), f1=inf; end
m=length(dirs); assert(m>0); sep=filesep;

for d=1:m, dir1=dirs{d}; dir1(dir1=='\')=sep; dir1(dir1=='/')=sep;
  if(dir1(end)==sep), dir1(end)=[]; end; dirs{d}=dir1; end

[fs0,fs1] = getFiles0(dirs{1},f0,f1,sep);
n1=length(fs0); fs=cell(m,n1); fs(1,:)=fs1;
for d=2:m, fs(d,:)=getFiles1(dirs{d},fs0,sep); end

  function [fs0,fs1] = getFiles0( dir1, f0, f1, sep )
    % get fs1 in dir1 (and fs0 without path or extension)
    fs1=dir([dir1 sep '*']); fs1={fs1.name}; fs1=fs1(3:end);
    fs1=fs1(f0:min(f1,end)); fs0=fs1; n=length(fs0);
    if(n==0), error('No files found in baseline dir %s.',dir1); end
    for i=1:n, fs1{i}=[dir1 sep fs0{i}]; end
    n=length(fs0); for i=1:n, f=fs0{i};
      f(find(f=='.',1,'first'):end)=[]; fs0{i}=f; end
  end

  function fs1 = getFiles1( dir1, fs0, sep )
    % get fs1 in dir1 corresponding to fs0
    n=length(fs0); fs1=cell(1,n); i2=0; i1=0;
    fs2=dir(dir1); fs2={fs2.name}; n2=length(fs2);
    eMsg='''%s'' has no corresponding file in %s.';
    for i0=1:n, r=length(fs0{i0}); match=0;
      while(i2<n2), i2=i2+1; if(strcmpi(fs0{i0},fs2{i2}(1:min(end,r))))
          i1=i1+1; fs1{i1}=fs2{i2}; match=1; break; end; end
      if(~match), error(eMsg,fs0{i0},dir1); end
    end
    for i1=1:n, fs1{i1}=[dir1 sep fs1{i1}]; end
  end
end

function fs = copyFiles( fs, dirs )
% Copy corresponding files into given directories.
%
% Useful for splitting data into training, validation and testing sets.
% See also bbGt>getFiles for obtaining a set of corresponding files.
%
% USAGE
%  fs = bbGt( 'copyFiles', fs, dirs )
%
% INPUTS
%   fs        - {mxn} list of full file names in each dir
%   dirs      - {1xm} list of m target directories
%
% OUTPUTS
%   fs        - {mxn} list of full file names of copied files
%
% EXAMPLE
%
% See also bbGt, bbGt>getFiles
[m,n]=size(fs); assert(numel(dirs)==m); if(n==0), return; end
for d=1:m
  if(~exist(dirs{d},'dir')), mkdir(dirs{d}); end
  for i=1:n, f=fs{d,i}; j=[0 find(f=='/' | f=='\')]; j=j(end);
    fs{d,i}=[dirs{d} '/' f(j+1:end)]; copyfile(f,fs{d,i}); end
end
end

function [gt0,dt0] = loadAll( gtDir, dtDir, pLoad )
% Load all ground truth and detection bbs in given directories.
%
% Loads each ground truth (gt) annotation in gtDir and the corresponding
% detection (dt) in dtDir. gt and dt files must correspond according to
% getFiles(). Alternatively, dtDir may be a filename of a single text file
% that contains the detection results across all images.
%
% Each dt should be a text file where each row contains 5 numbers
% representing a bb (left/top/width/height/score). If dtDir is a text file,
% it should contain the detection results across the full set of images. In
% this case each row in the text file should have an extra leading column
% specifying the image id: (imgId/left/top/width/height/score).
%
% The output of this function can be used in bbGt>evalRes().
%
% USAGE
%  [gt0,dt0] = bbGt( 'loadAll', gtDir, [dtDir], [pLoad] )
%
% INPUTS
%  gtDir      - location of ground truth
%  dtDir      - [] optional location of detections
%  pLoad      - {} params for bbGt>bbLoad() (determine format/filtering)
%
% OUTPUTS
%  gt0        - {1xn} loaded ground truth bbs (each is a mx5 array of bbs)
%  dt0        - {1xn} loaded detections (each is a mx5 array of bbs)
%
% EXAMPLE
%
% See also bbGt, bbGt>getFiles, bbGt>evalRes

% get list of files
if(nargin<2), dtDir=[]; end
if(nargin<3), pLoad={}; end
if(isempty(dtDir)), fs=getFiles({gtDir}); gtFs=fs(1,:); else
  dtFile=length(dtDir)>4 && strcmp(dtDir(end-3:end),'.txt');
  if(dtFile), dirs={gtDir}; else dirs={gtDir,dtDir}; end
  fs=getFiles(dirs); gtFs=fs(1,:);
  if(dtFile), dtFs=dtDir; else dtFs=fs(2,:); end
end

% load ground truth
persistent keyPrv gtPrv; key={gtDir,pLoad}; n=length(gtFs);
if(isequal(key,keyPrv)), gt0=gtPrv; else gt0=cell(1,n);
  for i=1:n, [~,gt0{i}]=bbLoad(gtFs{i},pLoad); end
  gtPrv=gt0; keyPrv=key;
end

% load detections
if(isempty(dtDir) || nargout<=1), dt0=cell(0); return; end
if(iscell(dtFs)), dt0=cell(1,n);
  for i=1:n, dt1=load(dtFs{i},'-ascii');
    if(numel(dt1)==0), dt1=zeros(0,5); end; dt0{i}=dt1(:,1:5); end
else
  dt1=load(dtFs,'-ascii'); if(numel(dt1)==0), dt1=zeros(0,6); end
  ids=dt1(:,1); assert(max(ids)<=n);
  dt0=cell(1,n); for i=1:n, dt0{i}=dt1(ids==i,2:6); end
end

end

function [gt,dt] = evalRes( gt0, dt0, thr, mul )
% Evaluates detections against ground truth data.
%
% Uses modified Pascal criteria that allows for "ignore" regions. The
% Pascal criteria states that a ground truth bounding box (gtBb) and a
% detected bounding box (dtBb) match if their overlap area (oa):
%  oa(gtBb,dtBb) = area(intersect(gtBb,dtBb)) / area(union(gtBb,dtBb))
% is over a sufficient threshold (typically .5). In the modified criteria,
% the dtBb can match any subregion of a gtBb set to "ignore". Choosing
% gtBb' in gtBb that most closely matches dtBb can be done by using
% gtBb'=intersect(dtBb,gtBb). Computing oa(gtBb',dtBb) is equivalent to
%  oa'(gtBb,dtBb) = area(intersect(gtBb,dtBb)) / area(dtBb)
% For gtBb set to ignore the above formula for oa is used.
%
% Highest scoring detections are matched first. Matches to standard,
% (non-ignore) gtBb are preferred. Each dtBb and gtBb may be matched at
% most once, except for ignore-gtBb which can be matched multiple times.
% Unmatched dtBb are false-positives, unmatched gtBb are false-negatives.
% Each match between a dtBb and gtBb is a true-positive, except matches
% between dtBb and ignore-gtBb which do not affect the evaluation criteria.
%
% In addition to taking gt/dt results on a single image, evalRes() can take
% cell arrays of gt/dt bbs, in which case evaluation proceeds on each
% element. Use bbGt>loadAll() to load gt/dt for multiple images.
%
% Each gt/dt output row has a flag match that is either -1/0/1:
%  for gt: -1=ignore,  0=fn [unmatched],  1=tp [matched]
%  for dt: -1=ignore,  0=fp [unmatched],  1=tp [matched]
%
% USAGE
%  [gt, dt] = bbGt( 'evalRes', gt0, dt0, [thr], [mul] )
%
% INPUTS
%  gt0  - [mx5] ground truth array with rows [x y w h ignore]
%  dt0  - [nx5] detection results array with rows [x y w h score]
%  thr  - [.5] the threshold on oa for comparing two bbs
%  mul  - [0] if true allow multiple matches to each gt
%
% OUTPUTS
%  gt   - [mx5] ground truth results [x y w h match]
%  dt   - [nx6] detection results [x y w h score match]
%
% EXAMPLE
%
% See also bbGt, bbGt>compOas, bbGt>loadAll

% get parameters
if(nargin<3 || isempty(thr)), thr=.5; end
if(nargin<4 || isempty(mul)), mul=0; end

% if gt0 and dt0 are cell arrays run on each element in turn
if( iscell(gt0) && iscell(dt0) ), n=length(gt0);
  assert(length(dt0)==n); gt=cell(1,n); dt=gt;
  for i=1:n, [gt{i},dt{i}] = evalRes(gt0{i},dt0{i},thr,mul); end; return;
end

% check inputs
if(isempty(gt0)), gt0=zeros(0,5); end
if(isempty(dt0)), dt0=zeros(0,5); end
assert( size(dt0,2)==5 ); nd=size(dt0,1);
assert( size(gt0,2)==5 ); ng=size(gt0,1);

% sort dt highest score first, sort gt ignore last
[~,ord]=sort(dt0(:,5),'descend'); dt0=dt0(ord,:);
[~,ord]=sort(gt0(:,5),'ascend'); gt0=gt0(ord,:);
gt=gt0; gt(:,5)=-gt(:,5); dt=dt0; dt=[dt zeros(nd,1)];

% Attempt to match each (sorted) dt to each (sorted) gt
oa = compOas( dt(:,1:4), gt(:,1:4), gt(:,5)==-1 );
for d=1:nd
  bstOa=thr; bstg=0; bstm=0; % info about best match so far
  for g=1:ng
    % if this gt already matched, continue to next gt
    m=gt(g,5); if( m==1 && ~mul ), continue; end
    % if dt already matched, and on ignore gt, nothing more to do
    if( bstm~=0 && m==-1 ), break; end
    % compute overlap area, continue to next gt unless better match made
    if(oa(d,g)<bstOa), continue; end
    % match successful and best so far, store appropriately
    bstOa=oa(d,g); bstg=g; if(m==0), bstm=1; else bstm=-1; end
  end; g=bstg; m=bstm;
  % store type of match for both dt and gt
  if(m==-1), dt(d,6)=m; elseif(m==1), gt(g,5)=m; dt(d,6)=m; end
end

end

function [hs,hImg] = showRes( I, gt, dt, varargin )
% Display evaluation results for given image.
%
% USAGE
%  [hs,hImg] = bbGt( 'showRes', I, gt, dt, varargin )
%
% INPUTS
%  I          - image to display, image filename, or []
%  gt         - first output of evalRes()
%  dt         - second output of evalRes()
%  varargin   - additional parameters (struct or name/value pairs)
%   .evShow     - [1] if true show results of evaluation
%   .gtShow     - [1] if true show ground truth
%   .dtShow     - [1] if true show detections
%   .cols       - ['krg'] colors for ignore/mistake/correct
%   .gtLs       - ['-'] line style for gt bbs
%   .dtLs       - ['--'] line style for dt bbs
%   .lw         - [3] line width
%
% OUTPUTS
%  hs         - handles to bbs and text labels
%  hImg       - handle for image graphics object
%
% EXAMPLE
%
% See also bbGt, bbGt>evalRes
dfs={'evShow',1,'gtShow',1,'dtShow',1,'cols','krg',...
  'gtLs','-','dtLs','--','lw',3};
[evShow,gtShow,dtShow,cols,gtLs,dtLs,lw]=getPrmDflt(varargin,dfs,1);
% optionally display image
if(ischar(I)), I=imread(I); end
if(~isempty(I)), hImg=im(I,[],0); title(''); end
% display bbs with or w/o color coding based on output of evalRes
hold on; hs=cell(1,1000); k=0;
if( evShow )
  if(gtShow), for i=1:size(gt,1), k=k+1;
      hs{k}=bbApply('draw',gt(i,1:4),cols(gt(i,5)+2),lw,gtLs); end; end
  if(dtShow), for i=1:size(dt,1), k=k+1;
      hs{k}=bbApply('draw',dt(i,1:5),cols(dt(i,6)+2),lw,dtLs); end; end
else
  if(gtShow), k=k+1; hs{k}=bbApply('draw',gt(:,1:4),cols(3),lw,gtLs); end
  if(dtShow), k=k+1; hs{k}=bbApply('draw',dt(:,1:5),cols(3),lw,dtLs); end
end
hs=[hs{:}]; hold off;
end

function [xs,ys,score,ref] = compRoc( gt, dt, roc, ref )
% Compute ROC or PR based on outputs of evalRes on multiple images.
%
% ROC="Receiver operating characteristic"; PR="Precision Recall"
% Also computes result at reference points (ref):
%  which for ROC curves is the *detection* rate at reference *FPPI*
%  which for PR curves is the *precision* at reference *recall*
% Note, FPPI="false positive per image"
%
% USAGE
%  [xs,ys,score,ref] = bbGt( 'compRoc', gt, dt, roc, ref )
%
% INPUTS
%  gt         - {1xn} first output of evalRes() for each image
%  dt         - {1xn} second output of evalRes() for each image
%  roc        - [1] if 1 compue ROC else compute PR
%  ref        - [] reference points for ROC or PR curve
%
% OUTPUTS
%  xs         - x coords for curve: ROC->FPPI; PR->recall
%  ys         - y coords for curve: ROC->TP; PR->precision
%  score      - detection scores corresponding to each (x,y)
%  ref        - recall or precision at each reference point
%
% EXAMPLE
%
% See also bbGt, bbGt>evalRes

% get additional parameters
if(nargin<3 || isempty(roc)), roc=1; end
if(nargin<4 || isempty(ref)), ref=[]; end
% convert to single matrix, discard ignore bbs
nImg=length(gt); assert(length(dt)==nImg);
gt=cat(1,gt{:}); gt=gt(gt(:,5)~=-1,:);
dt=cat(1,dt{:}); dt=dt(dt(:,6)~=-1,:);
% compute results
if(size(dt,1)==0), xs=0; ys=0; score=0; ref=ref*0; return; end
m=length(ref); np=size(gt,1); score=dt(:,5); tp=dt(:,6);
[score,order]=sort(score,'descend'); tp=tp(order);
fp=double(tp~=1); fp=cumsum(fp); tp=cumsum(tp);
if( roc )
  xs=fp/nImg; ys=tp/np; xs1=[-inf; xs]; ys1=[0; ys];
  for i=1:m, j=find(xs1<=ref(i)); ref(i)=ys1(j(end)); end
else
  xs=tp/np; ys=tp./(fp+tp); xs1=[xs; inf]; ys1=[ys; 0];
  for i=1:m, j=find(xs1>=ref(i)); ref(i)=ys1(j(1)); end
end
end

function [Is,scores,imgIds] = cropRes( gt, dt, imFs, varargin )
% Extract true or false positives or negatives for visualization.
%
% USAGE
%  [Is,scores,imgIds] = bbGt( 'cropRes', gt, dt, imFs, varargin )
%
% INPUTS
%  gt         - {1xN} first output of evalRes() for each image
%  dt         - {1xN} second output of evalRes() for each image
%  imFs       - {1xN} name of each image
%  varargin   - additional parameters (struct or name/value pairs)
%   .dims       - ['REQ'] target dimensions for extracted windows
%   .pad        - [0] padding amount for cropping
%   .type       - ['fp'] one of: 'fp', 'fn', 'tp', 'dt'
%   .n          - [100] max number of windows to extract
%   .show       - [1] figure for displaying results (or 0)
%   .fStr       - ['%0.1f'] label{i}=num2str(score(i),fStr)
%   .embed      - [0] if true embed dt/gt bbs into cropped windows
%
% OUTPUTS
%  Is         - [dimsxn] extracted image windows
%  scores     - [1xn] detection score for each bb unless 'fn'
%  imgIds     - [1xn] image id for each cropped window
%
% EXAMPLE
%
% See also bbGt, bbGt>evalRes
dfs={'dims','REQ','pad',0,'type','fp','n',100,...
  'show',1,'fStr','%0.1f','embed',0};
[dims,pad,type,n,show,fStr,embed]=getPrmDflt(varargin,dfs,1);
N=length(imFs); assert(length(gt)==N && length(dt)==N);
% crop patches either in gt or dt according to type
switch type
  case 'fn', bbs=gt; keep=@(bbs) bbs(:,5)==0;
  case 'fp', bbs=dt; keep=@(bbs) bbs(:,6)==0;
  case 'tp', bbs=dt; keep=@(bbs) bbs(:,6)==1;
  case 'dt', bbs=dt; keep=@(bbs) bbs(:,6)>=0;
  otherwise, error('unknown type: %s',type);
end
% create ids that will map each bb to correct name
ms=zeros(1,N); for i=1:N, ms(i)=size(bbs{i},1); end; cms=[0 cumsum(ms)];
ids=zeros(1,sum(ms)); for i=1:N, ids(cms(i)+1:cms(i+1))=i; end
% flatten bbs and keep relevent subset
bbs=cat(1,bbs{:}); K=keep(bbs); bbs=bbs(K,:); ids=ids(K); n=min(n,sum(K));
% reorder bbs appropriately
if(~strcmp(type,'fn')), [~,ord]=sort(bbs(:,5),'descend'); else
  if(size(bbs,1)<n), ord=randperm(size(bbs,1)); else ord=1:n; end; end
bbs=bbs(ord(1:n),:); ids=ids(ord(1:n));
% extract patches from each image
if(n==0), Is=[]; scores=[]; imgIds=[]; return; end;
Is=cell(1,n); scores=zeros(1,n); imgIds=zeros(1,n);
if(any(pad>0)), dims1=dims.*(1+pad); rs=dims1./dims; dims=dims1; end
if(any(pad>0)), bbs=bbApply('resize',bbs,rs(1),rs(2)); end
for i=1:N
  locs=find(ids==i); if(isempty(locs)), continue; end; I=imread(imFs{i});
  if( embed )
    if(any(strcmp(type,{'fp','dt'}))), bbs1=gt{i};
    else bbs1=dt{i}(:,[1:4 6]); end
    I=bbApply('embed',I,bbs1(bbs1(:,5)==0,1:4),'col',[255 0 0]);
    I=bbApply('embed',I,bbs1(bbs1(:,5)==1,1:4),'col',[0 255 0]);
  end
  Is1=bbApply('crop',I,bbs(locs,1:4),'replicate',dims);
  for j=1:length(locs), Is{locs(j)}=Is1{j}; end;
  scores(locs)=bbs(locs,5); imgIds(locs)=i;
end; Is=cell2array(Is);
% optionally display
if(~show), return; end; figure(show); pMnt={'hasChn',size(Is1{1},3)>1};
if(isempty(fStr)), montage2(Is,pMnt); title(type); return; end
ls=cell(1,n); for i=1:n, ls{i}=int2str2(imgIds(i)); end
if(~strcmp(type,'fn'))
  for i=1:n, ls{i}=[ls{i} '/' num2str(scores(i),fStr)]; end; end
montage2(Is,[pMnt 'labels' {ls}]); title(type);
end

function oa = compOas( dt, gt, ig )
% Computes (modified) overlap area between pairs of bbs.
%
% Uses modified Pascal criteria with "ignore" regions. The overlap area
% (oa) of a ground truth (gt) and detected (dt) bb is defined as:
%  oa(gt,dt) = area(intersect(dt,dt)) / area(union(gt,dt))
% In the modified criteria, a gt bb may be marked as "ignore", in which
% case the dt bb can can match any subregion of the gt bb. Choosing gt' in
% gt that most closely matches dt can be done using gt'=intersect(dt,gt).
% Computing oa(gt',dt) is equivalent to:
%  oa'(gt,dt) = area(intersect(gt,dt)) / area(dt)
%
% USAGE
%  oa = bbGt( 'compOas', dt, gt, [ig] )
%
% INPUTS
%  dt       - [mx4] detected bbs
%  gt       - [nx4] gt bbs
%  ig       - [nx1] 0/1 ignore flags (0 by default)
%
% OUTPUTS
%  oas      - [m x n] overlap area between each gt and each dt bb
%
% EXAMPLE
%  dt=[0 0 10 10]; gt=[0 0 20 20];
%  oa0 = bbGt('compOas',dt,gt,0)
%  oa1 = bbGt('compOas',dt,gt,1)
%
% See also bbGt, bbGt>evalRes
m=size(dt,1); n=size(gt,1); oa=zeros(m,n);
if(nargin<3), ig=zeros(n,1); end
de=dt(:,[1 2])+dt(:,[3 4]); da=dt(:,3).*dt(:,4);
ge=gt(:,[1 2])+gt(:,[3 4]); ga=gt(:,3).*gt(:,4);
for i=1:m
  for j=1:n
    w=min(de(i,1),ge(j,1))-max(dt(i,1),gt(j,1)); if(w<=0), continue; end
    h=min(de(i,2),ge(j,2))-max(dt(i,2),gt(j,2)); if(h<=0), continue; end
    t=w*h; if(ig(j)), u=da(i); else u=da(i)+ga(j)-t; end; oa(i,j)=t/u;
  end
end
end

function oa = compOa( dt, gt, ig )
% Optimized version of compOas for a single pair of bbs.
%
% USAGE
%  oa = bbGt( 'compOa', dt, gt, ig )
%
% INPUTS
%  dt       - [1x4] detected bb
%  gt       - [1x4] gt bb
%  ig       - 0/1 ignore flag
%
% OUTPUTS
%  oa       - overlap area between gt and dt bb
%
% EXAMPLE
%  dt=[0 0 10 10]; gt=[0 0 20 20];
%  oa0 = bbGt('compOa',dt,gt,0)
%  oa1 = bbGt('compOa',dt,gt,1)
%
% See also bbGt, bbGt>compOas
w=min(dt(3)+dt(1),gt(3)+gt(1))-max(dt(1),gt(1)); if(w<=0),oa=0; return; end
h=min(dt(4)+dt(2),gt(4)+gt(2))-max(dt(2),gt(2)); if(h<=0),oa=0; return; end
i=w*h; if(ig),u=dt(3)*dt(4); else u=dt(3)*dt(4)+gt(3)*gt(4)-i; end; oa=i/u;
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function fem2d_pack_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests the BASIS_11_**_TEST routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  BASIS_11_T3_TEST tests BASIS_11_T3.\n' );
  fprintf ( 1, '  BASIS_11_T4_TEST tests BASIS_11_T4.\n' );
  fprintf ( 1, '  BASIS_11_T6_TEST tests BASIS_11_T6.\n' );

  basis_11_t3_test ( );

  basis_11_t4_test ( );

  basis_11_t6_test ( );

  return
end

author: johannesgerer

repo: jburkardt-m


function result = c_strcmp(str1,str2)
%C_STRCMP Compare strings like C function "strcmp"
% Function C_STRCMP compares two strings, and returns
% a -1 if str1 < str2, a 0 if str1 == str2, and a
% +1 if str1 > str2.

% Define variables:
%   diff      -- Logical array of string differences
%   msg       -- Error message
%   result    -- Result of function
%   str1      -- First string to compare
%   str2      -- Second string to compare
%   strings   -- Padded array of strings
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    02/25/07    S. J. Chapman        Original code

% Check for a legal number of input arguments.
msg = nargchk(2,2,nargin);
error(msg);

% Check to see if the arguments are strings
if ~(isstr(str1) & isstr(str2))
   error('Both str1 and str2 must both be strings!')
else

   % Pad strings
   strings = strvcat(str1,str2);

   % Compare strings
   diff = strings(1,:) ~= strings(2,:);
   if sum(diff) == 0

      % Strings match, so return a zero!
      result = 0;
   else
      % Find first difference between strings
      ival = find(diff);
      if strings(1,ival(1)) > strings(2,ival(1))
         result = 1;
      else
         result = -1;
      end
   end
end


text:
function outputImage = imfill(this, locations, nD)
% Fills image slice-wise; mimicks imfill in matlab functionality
%
%   Y = MrImage()
%   filledY = Y.imfill(locations, nD)
%
% This is a method of class MrImage.
%
% IN
%   locations
%           array of 2D coordinates or string 'holes' to fill all holes
%   nD      dimensionality to perform operation
%           '2d' = slicewise application, separate 2d images
%           '3d' = as volume
%
% OUT
%   outputImage
%           MrImage where data matrix is inflated
%
% EXAMPLE
%   Y = MrImage();
%   filledY = Y.imfill()
%   filledY = Y.imfill('holes')
%
%
%   See also MrImage imfill MrImage.imerode perform_unary_operation

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-08-04
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% update geometry-header with flipping left-right
if nargin < 2
    locations = 'holes';
end

if nargin < 3
    nD = '2d';
end

if isreal(this)
    outputImage = this.perform_unary_operation(...
        @(x) imfill(x, locations), nD);
else
    outputImage = this.abs.perform_unary_operation(...
        @(x) imfill(x, locations), nD);
end

author: translationalneuromodeling

repo: tapas


text:
function [sys,x0,str,ts] = scrqp(t,x,u,flag,T0,nr,nq,np)
% FBFW feedforward feedback controller
% Output of the controller is calculated follows:
%
%            R(z^-1)                Q(z^-1)
% U(z^-1) = ---------- * W(z^-1) - --------- * Y(z^-1)
%            P(z^-1)                P(z^-1)
%
% Orders of polynomials R, Q and P are determined by parameters nr, nq and np.
%
% [sys,x0,str,ts] = scrqp(t,x,u,flag,T0,nr,nq,np)
% States: x(1:nr) ... w(k-1), w(k-2), ...
%         x(nr+1:nr+nq) ... y(k-1), y(k-2), ...
%         x(nr+nq+1:nr+nq+np-1) ... u(k-2), u(k-3), ... (u(k-1) is input - see below)
% Inputs: u(1) ... w(k)
%         u(2) ... y(k)
%         u(2+1:2+nr+1) ... r0, r1, ....
%         u(2+nr+1+1:2+nr+1+nq+1) ... q0, q1, ...
%         u(2+nr+1+nq+1+1:2+nr+1+nq+1+np+1) ... p0, p1, ...
%         u(2+nr+1+nq+1+np+1+1) ... u(k-1) real input to the controlled system, can differ to u(k) computed
%                             in previous loop (case of saturation, ...)
%         T0 ... sample time
%         nr ... odrer of polynomial R (R has nr+1 coefficients)
%         nq ... odrer of polynomial R (R has nq+1 coefficients)
%         np ... odrer of polynomial R (R has np+1 coefficients)

if flag == 0    %initialization
   %zero initial conditions
   x0 = zeros(nr+nq+np-1,1);
   %information about this function
   sys(1) = 0;  %Number of continuous states.
   sys(2) = nr+nq+np-1;  %Number of discrete states.
   sys(3) = 1;  %Number of outputs.
   sys(4) = 2 + nr+1 + nq+1 + np+1 + 1;  %Number of inputs.
   sys(5) = 0;  %Reserved for root finding. Must be zero.
   sys(6) = 1;  %Direct feedthrough flag (1=yes, 0=no)
   sys(7) = 1;  %Number of sample times. This is the number of rows in TS.
   ts = [T0 0];
   str = [];

elseif (flag==2 | flag==3)
    W = [u(1); x(1:nr)];             %[w(k), w(k-1), w(k-2), ...]
    Y = [u(2); x(nr+1:nr+nq)];       %[y(k), y(k-1), y(k-2), ...]
    U1 = [u(2+nr+1+nq+1+np+1+1); x(nr+nq+1:nr+nq+np-1)];      %[u(k-1), u(k-2), u(k-3),...]
    R = u(2+1:2+nr+1);                      % [r0, r1, ...]
    Q = u(2+nr+1+1:2+nr+1+nq+1);            % [q0, q1, ...]
    p0 = u(2+nr+1+nq+1+1);
    P1 = u(2+nr+1+nq+1+2:2+nr+1+nq+1+np+1);  % [p0, p1, ...]

    if (flag==2)   %update states
        sys = [W(1:nr); Y(1:nq); U1(1:np-1)];  %[y(k-1) y(k-2) u(k-2)]

    else     %flag == 3   compute output
        u = 1/p0*( R'*W - Q'*Y - P1'*U1);
        sys = u;
    end

else

    sys=[];
end

author: Sable

repo: mcbench-benchmarks


text:
% test202 Read data from cube.db.

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-25 Created (Johan Meijdam, VORtech)

function test202()
   OtStart();
   htable = OtTableNew();
   OtTableSetName(htable, [testdir() 'cube.DB']);
   if OtTableOpen(htable) ~= 1
      error('OtTableOpen failed')
   end
   OtTableLast(htable);

   disp('------');
   disp('FIELDS');
   disp('------');
   arr = OtTableFields(htable);
   disp(arr);
   disp('------------');
   disp(['VALUES (' num2str(OtTableRecordCount(htable)) ')']);
   disp('------------');
   while ~OtTableBof(htable)
      arr = OtTableGetAll(htable);
      disp(arr);
      OtTablePrior(htable);
   end
   OtTableClose(htable);
   OtTableFree(htable);
   OtStop();
end


text:
function timewindow = get_timewindow(utdnum_stop, numMins, utdnum_start)
% GET_TIMEWINDOW
%
% Usage:  timewindow = get_timewindow([utdnum, [numMins, [utdnum_start]]])
%
% get_timewindow returns the Matlab datenumbers corresponding to the start and end of a time window
% This time window is aligned such that if the numMins of data requested is 10, the start and end times
% are aligned at 0, 10, 20, 30, 40 or 50 minutes past the hour, in such a way that the end time is the
% most recent value that does not exceed dnum.
%
% If numMins is not given, it is set to 10.
%
% Example 1:
% 	[timewindow.start,timewindow.stop]=gettartAndEndTimes(datenum(2007,04,23,17,13,13), 15)
%
%	In this case timewindow.start would be datenum(2007,04,23,16,45,00)
%	and timewindow.stop would be datenum(2007,04,23,17,00,00)
%
% Example 2:
% 	[timewindow.start,timewindow.stop]=getStartAndEndTimes(datenum(2007,04,23,17,13,13), 10)
%
%	In this case timewindow.start would be datenum(2007,04,23,17,00,00)
%	and timewindow.stop would be datenum(2007,04,23,17,10,00)
%
% Glenn Thompson, 2007

global PARAMS

if nargin == 0
	utdnum_stop = utnow();
end

if nargin < 2
	numMins = 10;
end

if nargin > 3
	disp('get_timewindow([utdnum_stop [, numMins, [utdnum_start]]]')
	return;
end

timewindow.stop  = boundaryBeforeDnum(utdnum_stop, numMins);
timewindow.start = timewindow.stop - numMins/1440;

if (nargin == 3)
	snum = boundaryBeforeDnum(utdnum_start, numMins);
	snum_array = snum: numMins/1440: timewindow.start;
	enum_array = snum_array + numMins/1440;
	timewindow.start = snum_array;
	timewindow.stop = enum_array;
end

%%%%%%%%%%%%%%%%%%%
function dnum1=boundaryBeforeDnum(dnum, numMins)
dayFraction		=	rem(dnum,1);
minutesIntoThisDay	=	dayFraction*1440;
dnum1 			= 	dnum - rem(minutesIntoThisDay, numMins)/1440;

function dnum1=boundaryAfterDnum(dnum, numMins)
dnum1 = boundaryBeforeDnum(dnum, numMins) + numMins/1440;

author: geoscience-community-codes

repo: GISMO


text:
function v = calcAxisVolume(odf,axis,radius,varargin)
% amount of orientations with a specific misorientation axis
%
% Syntax
%   vol = calcAxisVolume(odf,axis,radius)
%
% Input
%  odf  - @SO3Fun
%  axis - @vector3d / @Miller
%  radius - double
%
% Output
%  vol - volumeportion of all axes within the specified radius around axis
%
% See also
% plotAxisDistribution

% get resolution for quadrature
res = get_option(varargin,'resolution',min(radius/5,2.5*degree));

% find fundamental region
sym = properGroup(disjoint(odf.CS,odf.SS));
if odf.antipodal || check_option(varargin,'antipodal')
  sym = sym.Laue;
end

% define a grid for quadrature
h = equispacedS2Grid(sym.fundamentalSector,'resolution',res,varargin{:});

% find those within the ball
ind = angle(Miller(axis,sym),Miller(h,sym)) < radius+1e-5;
h = h(ind);

% remember volume of the ball
vol = nnz(ind)/numel(ind);

% compute axis distrubtion
w = calcAxisDistribution(odf,h,varargin{:});

v = min(mean(w) * vol,1);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function y=fcn_isConnected()
%% About:
% Method to check the connection with the robot.

%% Return value:
% y: true, if a connection between the robot and the PC is active.
%     false, if no connection is present.

    global isConnected;
    y=isConnected;
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function testing(video)

disp(video);

previousMethod = 'MSCNN'; % BasicCNN or MSCNN

load(fullfile([previousMethod 'net/'], video, 'net-epoch-20'));
net.layers{end} = struct('name', 'data_hat_sigmoid', ...
    'type', 'sigmoid');
net = vl_simplenn_move(net,'gpu');

load meanPixel;
meanPixel(:,:,4) =0;

half_size = 15;

imgDir = ['../SBMIDataset/' video '/input'];
resDir = fullfile([previousMethod '-result'], video);

grayDir = ['../' previousMethod '/result/', video];

mkdir(resDir);
images = [dir([imgDir '/*.jpg']); dir([imgDir '/*.png'])];
grayimages = dir(fullfile(grayDir, '*.png'));

for kk = 1 : numel(images)
    %fprintf('%d\n', kk);
    imagename = images(kk).name;
    grayname = grayimages(kk).name;
    im = single(imread(fullfile(imgDir, imagename)));
    if size(im,1) > 400 || size(im,2) > 400
        im = imresize(im, 0.5, 'nearest');
    end

    im(:,:,4) = single(imread(fullfile(grayDir, grayname)));

    im_large = padarray(im, [half_size, half_size], 'symmetric');
    im_large = bsxfun(@minus, im_large, meanPixel);
    im_large = gpuArray(im_large);

    %     tic;
    A = vl_simplenn(net, im_large);
    B = gather(A(end).x);

    %     toc;

    map_im = uint8(B * 255);
    imwrite(map_im, fullfile(resDir, imagename));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [f,fp] = stft( x, sz, hp, pd, w)

if isreal( x)

	% Defaults
	if nargin < 5
		w = 1;
	end
	if nargin < 4
		pd = 0;
	end
	if nargin < 3
		hp = sz/2;
	end

        extra = (length(x)-sz)/hp;
        padding = ceil(extra)*hp + sz - length(x);
	x = [x zeros( 1, padding)];
%	x = [zeros( 1, sz+pd) x zeros( 1, sz+pd)];

	% Pack frames into matrix
	s = zeros( sz, (length(x)-sz)/hp);
	j = 1;
	for i = sz:hp:length( x)
		s(:,j) = w .* x((i-sz+1):i).';
		j = j + 1;
	end

	% FFT it
	f = fft( s, sz+pd);

	% Chop redundant part
	f = f(1:end/2+1,:);

	% Return phase component if asked to
	if nargout == 2
		fp = angle( f);
		fp = cos( fp) + sqrt(-1)*sin( fp);
	end

% Inverse transform
else

	% Defaults
	if nargin < 5
		w = 1;
	end
	if nargin < 4
		pd = 0;
	end
	if nargin < 3
		hp = sz/2;
	end

	% Ignore padded part
	if length( w) == sz
		w = [w; zeros( pd, 1)];
	end

	% Overlap add/window/replace conjugate part
	f = zeros( 1, (size(x,2)-1)*hp+sz+pd);
	v = 1:sz+pd;
	for i = 1:size( x,2)
		f((i-1)*hp+v) = f((i-1)*hp+v) + ...
			(w .* real( ifft( [x(:,i); conj( x(end-1:-1:2,i))])))';
	end

	% Norm for overlap
	f = f / (sz/hp);
	f = f(sz+pd+1:end-sz-2*pd);
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function e = trimSegment(s,imSize,mrg)

% TRIMSEGMENT  Trim segment at image borders
%   TRIMSEGMENT(S,IMSIZE) trims the segment S at the image borders
%   specified by IMSIZE. S is a 4-vector containing the two segment's
%   endpoints. IMSIZE is a 2-vector with the image dimensions in pixels,
%   IMSIZE = [HSIZE,VSIZE]. The output segment has the same orientation as
%   the input one.
%
%   TRIMSEGMENT(...,MRG) restricts the image size to be smaller in MRG
%   pixels at its four borders.
%
%   See also PINHOLESEGMENT.

%   (c) 2008-2009 Joan Sola @ LAAS-CNRS

% input segment's endpoints
a = s(1:2);
b = s(3:4);

% image witdh and height
[w,h] = split(imSize);

if nargin < 3
    insq = inSquare([a b],[0 w 0 h]);
else
    insq = inSquare([a b],[0 w 0 h],mrg);
end

if all(insq) % both endpoints are in the image

    e = s; % return the segment unchanged

else % at least one endpoint is out of the image

    H = pp2hmgLin(a,b); % homogeneous line
    L = [1; 0; 0];  % left image border
    R = [1; 0;-w];  % right
    T = [0; 1; 0];  % top
    B = [0; 1;-h];  % bottom

    % intersections of infinite line with infinite borders
    HL = intersectHmgLin(H,L,1);
    HR = intersectHmgLin(H,R,1);
    HT = intersectHmgLin(H,T,1);
    HB = intersectHmgLin(H,B,1);

    % bring to image borders
    i = 1;
    if inInterval(HL(2),[0,h])
        e(i:i+1,1) = HL;
        i = 3;
    end
    if inInterval(HR(2),[0,h])
        e(i:i+1,1) = HR;
        i = 3;
    end
    if inInterval(HT(1),[0,w])
        e(i:i+1,1) = HT;
        i = 3;
    end
    if inInterval(HB(1),[0,w])
        e(i:i+1,1) = HB;
    end

    if insq(1) % endpoint a is in the image

        p = e(1:2);
        q = e(3:4);

        u = b - a;
        v = p - a;
        if any(u./v > 0)
            e(1:2) = a;
            e(3:4) = p;
        else
            e(1:2) = a;
            e(3:4) = q;
        end

    elseif insq(2) % endpoint b is in the image

        p = e(1:2);
        q = e(3:4);

        u = a - b;
        v = p - b;
        if any(u./v > 0)
            e(1:2) = p;
            e(3:4) = b;
        else
            e(1:2) = q;
            e(3:4) = b;
        end

    else % no endpoint is inside the image

        if i == 1 % no intersection with image borders
            % Segment is not visible
            e = [];

        else
            p = e(1:2);
            q = e(3:4);

            if i==3 && any((p-a)./(b-p) > 0)  % segment is visible
                % check orientations

                u = b - a;
                v = q - p;
                if any(u./v < 0)
                    e = e([3 4 1 2]); % match orientations
                end
            else % segment is not visible
                e = [];
            end
        end

    end
end

return

%% test
imsize = [10 10];
s{1}  = [1 2 3 4]';
s{2}  = [-1 2 3 4]';
s{3}  = [1 -2 3 4]';
s{4}  = [1 2 -3 4]';
s{5}  = [1 2 3 -4]';
s{6}  = [1 2 3 11]';
s{7}  = [1 2 11 4]';
s{8}  = [1 11 3 4]';
s{9}  = [11 2 3 4]';
s{10} = [-1 -2 13 14]';
s{11} = [13 14 -1 -2]';

for i=1:numel(s)
    s{i}'
    (trimSegment(s{i},imsize))'
end

%% test
lmin = 10;
mrg = 0;
imSize = [640;480];
cla
lh = line('color','c','linewidth',3);
yh = line('color','r','linestyle','--')

%%
s = randn(4,1)*600

set(lh,'xdata',s([1,3]),'ydata',s([2,4]))

sv = trimSegment(s,imSize)

if isempty(sv)
    set(yh,'xdata',[],'ydata',[])
else
    set(yh,'xdata',sv([1,3]),'ydata',sv([2,4]))
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
close all;
clear all;
clc;
rng('default');
png_export = true;
pdf_export = false;
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

mf = spx.graphics.Figures();

% Signal space
N = 256;
% Number of measurements
M = 64;
% Sparsity level
K = 10;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K);
% Generate bi-uniform signals
x = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
% Measurement vectors
y = Phi.apply(x);
% GOMP solver instance
solver = spx.pursuit.single.GOMP(Phi, K);
%solver.Verbose = true;
% Solve the sparse recovery problem
result = solver.solve(y);
% Solution vector
z = result.z;

stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
spx.commons.sparse.print_recovery_performance(stats);

mf.new_figure('GOMP solution');
subplot(411);
stem(x, '.');
title('Sparse vector');
subplot(412);
stem(z, '.');
title('Recovered sparse vector');
subplot(413);
stem(abs(x - z), '.');
title('Recovery error');
subplot(414);
stem(y, '.');
title('Measurement vector');

author: indigits

repo: sparse-plex


text:
function title = p04_title ( title )

%*****************************************************************************80
%
%% P04_TITLE returns the title of problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'F(X) = EXP ( X ) - 1 / ( 10 * X )^2';

  return
end


text:
function [ res ] = run_KCF_RANSAC( seq, res_path, bSaveImage )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
project_dir = fileparts(fileparts(fileparts(mfilename('fullpath'))));
addpath(fullfile(project_dir, 'trackers', 'KCF'));
addpath(fullfile(project_dir, 'trackers', 'KCF', 'mex'));

kernel_type = 'gaussian';
kernel.type = kernel_type;

features.gray = false;
features.hog = false;

padding = 1.5;  %extra area surrounding the target
lambda = 1e-4;  %regularization
output_sigma_factor = 0.1;  %spatial bandwidth (proportional to target)

interp_factor = 0.02;

kernel.sigma = 0.5;

kernel.poly_a = 1;
kernel.poly_b = 9;

features.hog = true;
features.hog_orientations = 9;
cell_size = 4;

target_sz = seq.init_rect(1,[4,3]);
pos = seq.init_rect(1,[2,1]) + floor(target_sz/2);
img_files = seq.s_frames;
video_path = [];
visualization = 0;

%call tracker function with all the relevant parameters
positions = tracker_ransac(video_path, img_files, pos, target_sz, ...
    padding, kernel, lambda, output_sigma_factor, interp_factor, ...
    cell_size, features, visualization);

%return results to benchmark, in a workspace variable
rects = [positions(:,2) - target_sz(2)/2, positions(:,1) - target_sz(1)/2];
rects(:,3) = target_sz(2);
rects(:,4) = target_sz(1);
res.type = 'rect';
res.res = rects;

end

author: flyers

repo: drone-tracking


text:
function [ zfilt ] = gaussfilt( t,z,sigma )
%Apply a Gaussian filter to a time series
%   Inputs: t = independent variable, z = data at points t, and
%       sigma = standard deviation of Gaussian filter to be applied.
%   Outputs: zfilt = filtered data.
%
%   written by James Conder. Aug 22, 2013

n = length(z);  % number of data
zfilt = 0*z;    % initalize output vector

%%% get distances between points for proper weighting
w = 0*t;
w(2:end-1) = 0.5*(t(3:end)-t(1:end-2));
w(1) = t(2)-t(1);
w(end) = t(end)-t(end-1);

%%% check if sigma smaller than data spacing
iw = find(w > 2*sigma, 1);
if ~isempty(iw)
    disp('WARNING: sigma smaller than half node spacing')
    disp('May lead to unstable result')
    iw = w > 2.5*sigma;
    w(iw) = 2.5*sigma;
    % this correction leaves some residual for spacing between 2-3sigma.
    % otherwise ok.
    % In general, using a Gaussian filter with sigma less than spacing is
    % a bad idea anyway...
end

%%% loop over points
a = 1/(sqrt(2*pi)*sigma);
sigma2 = sigma*sigma;
for i = 1:n
    filter = a*exp(-0.5*((t - t(i)).^2)/(sigma2));
    zfilt(i) = sum(w.*z.*filter);
end

%%% clean-up edges - mirror data for correction
ss = 2.4*sigma;   % distance from edge that needs correcting

% left edge
tedge = min(t);
iedge = find(t < tedge + ss);
nedge = length(iedge);
for i = 1:nedge;
    dist = t(iedge(i)) - tedge;
    include = find( t > t(iedge(i)) + dist);
    filter = a*exp(-0.5*((t(include) - t(iedge(i))).^2)/(sigma2));
    zfilt(iedge(i)) = zfilt(iedge(i)) + sum(w(include).*filter.*z(include));
end

% right edge
tedge = max(t);
iedge = find(t > tedge - ss);
nedge = length(iedge);
for i = 1:nedge;
    dist = tedge - t(iedge(i));
    include = find( t < t(iedge(i)) - dist);
    filter = a*exp(-0.5*((t(include) - t(iedge(i))).^2)/(sigma2));
    zfilt(iedge(i)) = zfilt(iedge(i)) + sum(w(include).*filter.*z(include));
end

end


text:
function [total_time, ts, X] = TrajectoryPlanning(path, decomp, time_allocation)
    speed   = time_allocation.avg_speed;
    acc     = time_allocation.acc;

    disp(['The planned speed is : ', num2str(speed)]);
    tic
    if strcmp(time_allocation.type, 'averageSpeed')
        [ts, total_time] = averageSpeed_ta(path0, speed);
    elseif strcmp(time_allocation.type, 'trapzoidSpeed')
        [ts, total_time] = trapezoidalSpeed_ta(path, speed, acc);
    end
    disp('TimeAllocation time is :');
    toc
    disp(['time management: total_time is ', num2str(total_time), 'seconds']);
    disp(['Split time is : ', num2str(ts)]);

    tic
    %  use Ax=b get Trajectory planning. ===========================
    %  X = Ax_equal_b(path, total_time, ts);

    % use 'Quadratic Programming' and SFC(which make by Ax < b) get Trajectory planning. ========================
    % use SFC make Inequality constraints
    X = QPbyUseSFC(path, ts, decomp);

    disp('generator trajectory time is :');
    toc
end

author: LenaShengzhen

repo: AerialRobotics


text:
function path = grFindGeodesicPath(nodes, edges, ind0, ind1, edgeWeights)
%GRFINDGEODESICPATH Find a geodesic path between two nodes in the graph.
%
%   PATH = grFindGeodesicPath(NODES, EDGES, NODE1, NODE2, WEIGHTS)
%   NODES and EDGES defines the graph, NODE1 and NODE2 are indices of the
%   node extremities, and WEIGHTS is the set of weights associated to each
%   edge.
%   The function returns a set of edge indices.
%
%
%   See also
%   grFindMaximalLengthPath

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-05-22, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% ensure weights are defined
if ~exist('edgeWeights', 'var')
    edgeWeights = ones(size(edges, 1), 1);
end

% check indices limits
nNodes = size(nodes, 1);
if max(ind0) > nNodes
    error('Start index exceed number of nodes in the graph');
end
if max(ind1) > nNodes
    error('End index exceed number of nodes in the graph');
end

% find a vertex opposite to the first extremity
dists = grPropagateDistance(nodes, edges, ind0, edgeWeights);

% iterate on neighbors of current node: choose next neighbor with smallest
% cumulated weight, until we are back on source node
path = [];
while true
    % find neighbor with lowest cumulated distance
    neighs = grAdjacentNodes(edges, ind1);
    neighDists = dists(neighs);
    indN = find(neighDists == min(neighDists), 1);
    ind2 = neighs(indN);

    if isempty(ind2)
        warning('graphs:grFindGeodesicPath', ...
            'No neighbor node found for node %d, graph may be not connected', ind1);
        break;
    end

    % add edge index to the path
    indE = find(sum(ismember(edges, [ind1 ind2]), 2) == 2, 1);
    path = [path indE]; %#ok<AGROW>

    % test if path is finished or not
    if ind2 == ind0
        break;
    end
    ind1 = ind2;
end

% reverse path direction
path = path(end:-1:1);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
clear; close all; clc;

% Datasets = {'TUM', 'KITTI', 'Tanks_and_Temples', 'CPC'};
Datasets = {'TUM'};

Methods = {'SIFT-RT-RANSAC'};

Errors = cell(length(Methods),length(Datasets));
Inlier_rates = cell(length(Methods),length(Datasets));
Numbers = cell(length(Methods),length(Datasets));
for d = 1 : length(Datasets)
    dataset = Datasets{d};
    for m = 1 : length(Methods)
        method = Methods{m};

        results_dir = ['../Results/' dataset '/'];
        filename = [results_dir method '.mat'];
        Results = importdata(filename);

        Error = -ones(length(Results), 1);
        Inlier_rate = -ones(length(Results), 2);
        Number =  zeros(length(Results), 2);

        for idx = 1 : length(Results)

            if Results{idx}.status ~=0
                Results{idx}.sgd_error = -1;
                Results{idx}.inlier_rate = [0,0];
                continue;
            end

            F1 = Results{idx}.F_gt;
            F2 = Results{idx}.F_hat;
            size1 = Results{idx}.size_l;
            size2 = Results{idx}.size_r;
            X1 = Results{idx}.X_l';
            X2 = Results{idx}.X_r';
            inliers = Results{idx}.inliers;

            if isfield(Results{idx}, 'sgd_error') ~= 1 || Results{idx}.sgd_error < 0
                Results{idx}.sgd_error = ComputeNormlizedSGD(F1, F2, size1, size2);
            end
            Error(idx) = Results{idx}.sgd_error;

            if isfield(Results{idx}, 'inlier_rate') ~= 1 || isempty(Results{idx}.inlier_rate) == 1
                Results{idx}.inlier_rate = ComputeInlierRate(F1, X1, X2, inliers, size1, size2, 0.003);
            end
            Inlier_rate(idx,:) = Results{idx}.inlier_rate;
            Number(idx,:) = [length(Results{idx}.inliers), sum(Results{idx}.inliers)];
        end

        save(filename, 'Results');

        mask = Error < 0;
        Error(mask) = [];
        Errors{m, d} = Error;
        Inlier_rate(mask,:) = [];
        Inlier_rates{m, d} = Inlier_rate;
        Number(mask, :) = [];
        Numbers{m, d} = Number;
    end
end

% Recall---(Error)
num_pairs = 1000;
X = linspace(0,0.2,20);
for d = 1 : length(Datasets)
    dataset = Datasets{d};
    Y = zeros(length(Methods), length(X));
    for m = 1 : length(Methods)
       method = Methods{m};
       for t = 1 : length(X)
           Y(m, t) = sum(Errors{m,d} < X(t)) / num_pairs;
       end
    end
    figure;
    h = plot(X,Y,'linewidth',3);
    ylim([0 1]);
    legend(h, Methods, 'Location', 'SouthEast');
    title(dataset);
    xlabel('NSGD Threshold');
    ylabel('Recall');
end

threshold = 0.05;
for d = 1 : length(Datasets)
    dataset = Datasets{d};
    disp(['Dataset : ' dataset]);
    disp('method recall inlier_rate_before inlier_rate_after');

    for m = 1 : length(Methods)
       method = Methods{m};

       recall = sum(Errors{m,d} < threshold) / num_pairs;

       [meanInlierRate] = mean(Inlier_rates{m,d});
       before_rate = meanInlierRate(1);
       after_rate = meanInlierRate(2);

       fprintf(sprintf('%s %f %f %f\n', method, recall, before_rate, after_rate));
    end
end

author: JiawangBian

repo: FM-Bench


text:
clear
patchSize = 8;

I = double(rgb2gray(imread('160068.jpg')))/255;

% load blurring kernel (you can download the kernels from Dilip Krishnan's
% website
% load kernels.mat
% K = kernel1;
K = fspecial('motion',10,45);
% K = fspecial('gaussian',[5 5],1);
noiseSD = 0.01;
patchSize = 8;

% convolve with kernel and add noise
ks = floor((size(K, 1) - 1)/2);
yorig = I;
y = conv2(yorig, K, 'valid');
y = y + noiseSD*randn(size(y));
y = double(uint8(y .* 255))./255;

% code excerpt taken from Krishnan et al.

% edgetaper to better handle circular boundary conditions
y = padarray(y, [1 1]*ks, 'replicate', 'both');
for a=1:4
  y = edgetaper(y, K);
end

noiseI = y;

% load GMM model
load GSModel_8x8_200_2M_noDC_zeromean.mat

% uncomment this line if you want the total cost calculated
% LogLFunc = @(Z) GMMLogL(Z,GS);

% initialize prior function handle
excludeList = [];
prior = @(Z,patchSize,noiseSD,imsize) aprxMAPGMM(Z,patchSize,noiseSD,imsize,GS,excludeList);

% comment this line if you want the total cost calculated
LogLFunc = [];

% deblur
tic
[cleanI,psnr,~] = EPLLhalfQuadraticSplitDeblur(noiseI,64/noiseSD^2,K,patchSize,50*[1 2 4 8 16 32 64],1,prior,I,LogLFunc);
toc

% output result
figure(1);
imshow(I); title('Original');
figure(2);
imshow(noiseI); title('Corrupted Image');
figure(3);
imshow(cleanI); title('Restored Image');

author: lbasek

repo: image-denoising-benchmark


text:
% Assume a target positioned at x = 1, travelling with speed v = 0.1
state = [1;0.1;2;0;3;0.2;4;0.3];

% Create an instance of a 2D Ornstein-Uhlenbeck model
ou = OrnsteinUhlenbeckModelX('NumDims',4,...
                             'VelocityErrVariance',0.1,...
                             'DampingCoefficient',0.1,...
                             'TimestepDuration',duration(1,0,1));

% View the transition matrix and process covariance matrices
F = ou.feval();
Q = ou.covar();

% Predict the target's position and velocity after the interval has passed
newState  = ou.feval(state);

% Do the same as above, but this time add process noise to the prediction
newState2 = ou.feval(state,true);

% Generate 50 random noise samples from the dynamic model
noise = ou.random(50);

% Check how likely the predictions we made are
lik = ou.pdf(newState,state);
lik2 = ou.pdf(newState2,state); % HINT: newState2 should be less likely

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [sfmL,sfm,names] = spm2delta(SPM,varargin)
% [d_hires,d_atTR,names] = spm2delta(SPM,varargin)
% or
% [onsets,delta_atTR,names] = spm2delta(SPM,varargin)
%
% Example: for SPM99
% [dL,d] = spm2delta([],Sess);

if length(varargin) > 0,
    Sess = varargin{1};,
    SPM.Sess = Sess;
else
    % SPM2
    Sess = SPM.Sess;
end

if iscell(Sess)

    ons = {};

    for i = 1:length(Sess)

        [sfc,sfm{i}] = downsample_delta(Sess{i}.sf,16);

        for j = 1:size(sfm{i},2)                % save onsets
            ons{end+1} = find(sfm{i}(:,j)) - 1; % in TRs, starting with 0
        end

        %sfmL{i} = cat(2,Sess{i}.sf{:});        % delta hires
    end

    sfm = cat(1,sfm{:});
    sfmL = ons;         % save onsets
    %sfmL = cat(1,sfmL{:}); % save delta_hires
    names = SPM.Sess{1}.name;
else

    % new spm2 way?
    ons = [];
    wh = [];, for i = 1:length(Sess), if isempty(Sess(i).U), wh(i) = 1;,end,end
        Sess(find(wh)) = [];

    delta = cell(1,length(Sess(1).U));

    for i = 1:length(Sess)
        sessons = [];
        for j = 1:length(Sess(i).U)
            ons{end+1} = Sess(i).U(j).ons/SPM.xY.RT;      % all onsets
            sessons{end+1} = Sess(i).U(j).ons/SPM.xY.RT;  % onsets this session
        end

        % add delta functions for this session to overall
        sessdelt = onsets2delta(sessons,size(SPM.xX.X,1) ./ length(Sess));
        for k = 1:size(sessdelt,2)

            % just in case some sessions do not have some
                    % regressors
            if k > length(delta)
            	delta{k} = zeros(size(delta{1}));
            end

            try
                delta{k} = cat(1,[delta{k}; sessdelt(:,k)]);
            catch
                 delta{k} = cat(1,[delta{k}; sessdelt(:,k)']);
            end
        end
    end
    sfmL = ons;
    sfm = delta;

    names = cat(2,SPM.Sess(1).U(:).name);
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function F = uminus(F)
%UMINUS   Unary minus for CHEBFUN3 objects.
%   -F negates the CHEBFUN3 object F.
%
%   G = UMINUS(F) is called for the syntax '-F'.
%
% See also CHEBFUN3/UPLUS and CHEBFUN3/MINUS.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

F.core = -F.core;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function A = sim_genTVMVARStructure(Aproto,Nl,ndisc,p)
% generate time-varying MVAR structure from a prototype structure
%
% Inputs:
%
%       Aproto:   Prototype structure returned by sim_genVARModelFromEq().
%       Nl:       Number of sample points
%       ndisc:    Number of 'startup' points we will discard when generating data
%       p:        MVAR model order
%
% Outputs:
%
%       A:        Cell array of dimension Nl+ndisc+p containing VAR
%                 coefficient matrices for each time point
%
% See Also: sim_genVARModelFromEq()
%
% References:
%
%
% Author: Tim Mullen, May 2011, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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 3 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

    A = repmat({Aproto},1,Nl+ndisc+p);

    if iscell(Aproto)
        % find all the function handles
        numidx = find(cellfun(@(x)isnumeric(x),Aproto));
        funidx = setdiff_bc(1:numel(Aproto),numidx);

        % funidx = find(cellfun(@(x)strcmpi(class(x),'function_handle'),Aproto));

        % construct inline objects
        for fun = 1:length(funidx)
            fx{fun} = inline(Aproto{funidx(fun)});
        end

        % evaluate each function handle
        for t=1:Nl+ndisc+p;

            if ~mod(t,1000)
                fprintf('%d/%d - ',t,Nl+ndisc+p);
            end

            for fun = 1:length(funidx)
                A{t}{funidx(fun)} = double(fx{fun}(t));
            end

            A{t} = cell2mat(A{t});
        end
    end
end

repo: REST


text:
function res = is_symbol(x)
% check whether an expression is a symbol
res = isequal(exp_head(x),@Symbol);


text:
function value = p02_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P02_F evaluates the integrand for problem 02.
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integrand:
%
%    ( sum ( 2 * x(1:dim_num) - 1 ) )**4
%
%  Exact Integral:
%
%    DIM_NUM * (5*DIM_NUM-2) / 15
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  value(1:point_num) = 0.0;

  for point = 1 : point_num
    value(point) = ( sum ( 2.0 * x(1:dim_num,point) - 1.0 ) )^4;
  end

  p02_i4 ( 'I', '#', point_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = prtRvUtilMvtLogPdf(x,mu,Sigma,dof)
% y = prtRvUtilMvtLogPdf(x,mu,Sigma,dof)
% xxx Need Help xxx

d = size(x,2);

[R,err] = cholcov(Sigma,0);
if err ~= 0
    error('mvtLogPdf:BadCovariance', ...
        'SIGMA must be symmetric and positive definite.');
end

% Create array of standardized data
xRinv = bsxfun(@minus,x,mu(:)') / R;

c = gammaln((d+dof)*0.5) - (d*0.5)*log(dof*pi) - gammaln(dof*0.5) - 0.5*log(det(Sigma));

y = c - (dof+d)*0.5*log(1 +1/dof*sum(xRinv.^2, 2));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function S = struct_mfile_reference(mfile,struct_name,strs)
% S = struct_mfile_reference(mfile,struct_name,strs)

% AF 11/1/01

pattern = [struct_name '.'];
if (~isempty(mfile))
   strs = grep(mfile,pattern);
end
S = [];

for i = 1:length(strs)
   inds = strfind(strs{i},pattern);
   for j = inds(:)'
      start_field = j+length(pattern);
      k = start_field+1;
      while ((k <= length(strs{i})) & (isvarname(strs{i}(start_field:k))))
         k = k+1;
      end
      field_name = strs{i}(start_field:k-1);
      if (~isempty(field_name))
         eval(['S.' field_name ' = struct_mfile_reference(''' mfile ''',''' struct_name '.' field_name ''',strs);']);
      end
   end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   example 4 DOF planar robot.
%
%   Author: Arturo Gil. Universidad Miguel Hern�ndez de Elche.
%   email: arturo.gil@umh.es date:   05/03/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function robot = parameters()

robot.name='Example 4DOF planar arm';

%kinematic data DH parameters
robot.DH.theta='[q(1) q(2) q(3) q(4)]';
robot.DH.d='[0  0  0  0]';
robot.DH.a='[0.4  0.3  0.2  0.1]';
robot.DH.alpha='[0  0  0  0]';

%number of degrees of freedom
robot.DOF = 4;

%Jacobian matrix
% if not defined, compute it with manipulator_jacobian
robot.J=[];
%the proper rows of J to compute manipulability
robot.selJ =[1,2,6];
robot.kind=['R' 'R' 'R' 'R'];

%Function name to compute inverse kinematic
robot.inversekinematic_fn = 'inverse_kinematics_4dofplanar(robot, T, q)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-180) deg2rad(180); %Axis 1, minimum, maximum
                deg2rad(-180) deg2rad(180);
                deg2rad(-180) deg2rad(180)]; %Axis 2, minimum, maximum

            %minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-90) deg2rad(90); %Axis 1, minimum, maximum
                deg2rad(-90) deg2rad(90);
                deg2rad(-180) deg2rad(180)]; %Axis 2, minimum, maximum

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = []; %empty, not available

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time
% end effectors maximum velocity
robot.linear_velmax = 0; %m/s, example, not available

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

%this is an ad hoc transformation between the piece and the robot's end
%effector
% robot.Tcoupling=[-1 0 0  0;
%      0 1 0  0;
%      0 0 -1 0;
%      0 0 0 1];
% Tc1 = [cos(-pi/2) 0 sin(-pi/2) 0;
%        0 1 0 0;
%        -sin(-pi/2) 0 cos(-pi/2) 0;
%        0  0  0  1       ];
% Tc2 = [1 0 0 0;
%        0 cos(pi) -sin(pi) 0;
%        0  sin(pi) cos(pi) 0;
%        0 0 0 1];
% robot.Tcoupling =Tc1*Tc2;
robot.Tcoupling = eye(4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GRAPHICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read graphics files
robot.graphical.has_graphics=1;
robot.graphical.color = [255 20 40]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-3.5 3.5 -3.5 3.5 0 1]
robot = read_graphics(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DYNAMIC PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[1 1 1] ;

%COM of each link with respect to own reference system
robot.dynamics.r_com=[-0.5      0         0; %(rx, ry, rz) link 1
                      -0.5      0         0;
                      -0.5      0         0];%(rx, ry, rz) link 2

%link masses
m1 = robot.dynamics.masses(1);
m2 = robot.dynamics.masses(2);
m3 = robot.dynamics.masses(3);

a=eval(robot.DH.a);
L1 = a(1);
L2 = a(2);
L3 = a(3);

%Momentos de inercia de cada eslabon.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, por cada fila
robot.dynamics.Inertia=[0   m1*L1^2/3   m1*L1^2/3    0	0	0;
                        0   m2*L2^2/3   m2*L2^2/3    0	0	0;
                        0   m3*L3^2/3   m3*L3^2/3    0	0	0];

%Actuator rotor inertia
robot.motors.Inertia=[0 0 0];
%Reduction ratio motor/joint speed
robot.motors.G=[1 1 1];
%consider friction
robot.friction=0;
%Viscous friction factor, motor referred
%robot.B = [1e-3  1e-3 1e-3];
robot.motors.Viscous = [10  10 10];
%Coulomb friction, motor referred
robot.motors.Coulomb = [0 0;
            0 0;
            0 0];

%SPECIAL PARAMETERS TO SOLVE THE INVERSE KINEMATICS
robot.parameters.step_time=0.1;
%Error in XYZ to stop inverse kinematics
robot.parameters.epsilonXYZ=0.01;
%Error in Quaternion to stop inverse kinematics.
robot.parameters.epsilonQ=0.01;
robot.parameters.stop_iterations=500;

repo: ARTE


text:
function Obj = SOFAconvertTHK2SOFA(miroObj)
%SOFAconvertTHK2SOFA - converts from miroObj to SOFA format
%   Usage: OBJ=SOFAconvertTHK2SOFA(miroObj)
%
%   SOFAconvertTHK2SOFA(miroObj) converts the HRIRs, BRIRs, and DRIRs (VariSphear array measurements) described in miroObj to SOFA. miroObj is the miro object saved at the Technische Hochschule Koeln, provided by Benjamin Bernschuetz.
%
%   Input parameters:
%     miroObj : HRTF data in miro format
%
%   Output parameters:
%     Obj : New SOFA object (SOFA format)
%
%   Reference to the source format: http://www.audiogroup.web.th-koeln.de/FILES/miro_documentation.pdf
%   Reference to the source coordinate system: http://www.audiogroup.web.th-koeln.de/SOFiA_wiki/COORDINATES.html

% #Author: Tim Lübeck, TH Köln (2018)
% #Author: Piotr Majdak
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
%
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences;
% Licensed under the EUPL, Version 1.2 or – as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%%
%if isoctave,
if exist('OCTAVE_VERSION','builtin') ~= 0
    error(['Octave is not able to convert THK to SOFA, use Matlab instead.']);
end

%% Check if miro.m class file is available, if not download file from server
% miro class might also be included in other toolboxes, eg. AKtools
if exist('miro','class') ~= 8
    % download miro.m
    disp('Downloading miro.m from TH Köln server...')
    url = 'http://audiogroup.web.th-koeln.de/FILES/miro.m';
    basepath=which('SOFAstart');
    basepath=basepath(1:end-12); % Kill the function name from the path.
    target=[basepath filesep 'helpers' filesep 'miro.m'];
    websave (target,url);
end

%% Get an empty conventions structure depending on miro format
if ( strcmp (miroObj.type,'HRIR') )
    Obj = SOFAgetConventions('SimpleFreeFieldHRIR');
    Obj.GLOBAL_Title = 'HRIR';
    Obj.GLOBAL_RoomType = 'free field';
    typeFlag = 0;
elseif ( strcmp (miroObj.type,'BRIR') )
    Obj = SOFAgetConventions('MultiSpeakerBRIR');
    Obj.GLOBAL_Title = 'BRIR';
    Obj.GLOBAL_RoomType = 'reverberant';
    typeFlag = 1;
else %( strcmp (miroObj.type,'MICARRAY') )
    Obj = SOFAgetConventions('SingleRoomDRIR');
    Obj.GLOBAL_Title = 'DRIR';
    Obj.GLOBAL_RoomType = 'reverberant';
    typeFlag = 2;
end

%% Fill global attributes
Obj.GLOBAL_ListenerShortName = miroObj.name;
Obj.GLOBAL_AuthorContact = miroObj.contact;
Obj.GLOBAL_Comment = [miroObj.context,' / Sampling Grid: ',miroObj.quadGrid];
Obj.GLOBAL_History = SOFAappendText(Obj,'GLOBAL_History','Converted from the miro file format');
Obj.GLOBAL_License = 'CC 3.0 BY-SA';
Obj.GLOBAL_Organization = 'Technische Hochschule Koeln, Germany';
Obj.GLOBAL_Author = miroObj.engineer;
Obj.GLOBAL_Origin = 'http://audiogroup.web.th-koeln.de';
Obj.GLOBAL_DateCreated = datestr(datenum(miroObj.date),'yyyy-mm-dd HH:MM:SS');
Obj.GLOBAL_DatabaseName='THK';
Obj.GLOBAL_ListenerDescription = miroObj.microphone;
Obj.GLOBAL_ReceiverDescription = [miroObj.microphone '; ' miroObj.micPreamp];
Obj.GLOBAL_SourceDescription = miroObj.source;
Obj.GLOBAL_EmitterDescription = miroObj.source;
Obj.GLOBAL_RoomDescription = [miroObj.location,' / avgAirTemp: ',num2str(miroObj.avgAirTemp),' / avgRelHumidity: ',num2str(miroObj.avgRelHumidity)];

%% Set miroObj to degree mode
miroObj.shutUp = 1;
miroObj = setDEG(miroObj);

%% Get miroObject data and convert to SOFA structure
%M: number of measurements;
%R: number of receivers;
%N: number of data samples describing one measurement. Data is a function of N;
%E: number of emitters;
%C: coordinate dimension, always three with the meaning
%FIR : [M R N];

if (typeFlag == 0 || typeFlag == 1)  %BRIR or HRIR
    irChOne = zeros(miroObj.returnTaps,miroObj.nIr);
    irChTwo = zeros(miroObj.returnTaps,miroObj.nIr);
    for channel = 1 : miroObj.nIr
        IR = getIR( miroObj, channel );
        irChOne(:, channel) = IR(:,1);  %[N M]
        irChTwo(:, channel) = IR(:,2);  %[N M]
    end
    if (typeFlag == 0) %HRIR
        Obj.Data.IR = irChOne; % irChOne is [N M]
        Obj.Data.IR(:,:,2) = irChTwo;
        Obj.Data.IR = shiftdim(Obj.Data.IR,1); % convert from [N M R] to [M R N]
    else % BRIR
        Obj.Data.IR = zeros(size(irChOne,2), 2, 1, size( irChTwo,1)); %[N 2 1 M]
        Obj.Data.IR(:,1,:) = shiftdim(shiftdim(irChOne,-2),3);
        Obj.Data.IR(:,2,:) = shiftdim(shiftdim(irChTwo,-2),3);
    end
else % DRIR
    irData = zeros(miroObj.returnTaps,miroObj.nIr);
    for channel = 1 : miroObj.nIr
        irData(:,channel) = miroObj.getIR(channel);
    end
    Obj.Data.IR             = zeros(1, miroObj.nIr ,miroObj.returnTaps); %[M R N]
    for R = 1 : miroObj.nIr
        Obj.Data.IR(1,R,:) = irData(:,R);
    end
    Obj.Data.Delay = zeros(1,miroObj.nIr);
end
Obj.Data.SamplingRate = miroObj.fs;

%% Fill the mandatory source emitter variables
Obj.ListenerPosition = [0 0 0]; % for BRIR and HRIR listener in center
Obj.ReceiverPosition = [0 +miroObj.radius 0; 0 -miroObj.radius 0];  % for HRIR and BRIR ears as receiver

if (typeFlag == 0)%HRIR
    Obj.ListenerView = [1 0 0];
    Obj.ListenerUp = [0 0 1];
    Obj.SourcePosition = [...
        miroObj.azimuth' ... % azimuth angle in a range of (0-360°(. Whereas AZ=0° is defined to be the front direction and AZ=180° to be the rear direction.
        90-miroObj.elevation' ... % elevation angle in range of (0-180°). EL=0 points upwards, EL=90° points to the horizontal plane and EL=180° points downwards.
        miroObj.sourceDistance*ones(size(miroObj.azimuth'))]; % radius in meters
elseif (typeFlag == 1) %BRIR
    Obj.SourcePosition = [0 0 0]; % default edit manually!
    Obj.EmitterPosition = [miroObj.sourceDistance 0 0];   % default position is center, otherwise define manually
    Obj.EmitterPosition_Type  = 'cartesian';
    Obj.EmitterUp = [0 0 1];
    Obj.EmitterView = [-1 0 0];
    Obj.ListenerView = [miroObj.azimuth',  ...          % see HRIR definitions
                        90-miroObj.elevation',  ...
                        zeros(size(miroObj.azimuth'))]; %miroObj.sourceDistance*ones(size(miroObj.azimuth'))];
    Obj.ListenerView_Type = 'spherical';
    Obj.ListenerView_Units = 'degree, degree, metre';
    Obj.ListenerUp = [0 0 1];
else %DRIR
    Obj.SourcePosition  = [1,0,0]; % default edit manually!
    Obj.EmitterPosition = [0,0,0];
    Obj.ListenerPosition    = [0 0 0];
    Obj.ListenerView        = [1 0 0];
    Obj.ListenerUp          = [0 0 1];
    Obj.ListenerView_Type = 'spherical';
    Obj.ListenerView_Units = 'degree, degree, metre';
    Obj.ReceiverPosition = [...
                            miroObj.azimuth' ...
                            90-miroObj.elevation' ...
                            miroObj.radius*ones(size(miroObj.azimuth'))];
    Obj.ReceiverPosition_Type = 'spherical';
    Obj.ReceiverPosition_Units = 'degree, degree, metre';
end

%% Update dimensions
Obj=SOFAupdateDimensions(Obj);

author: sofacoustics

repo: SOFAtoolbox


text:
function [varargout]=curvePathOrderFix(V)

% function [Vf,indFix]=curvePathOrderFix(V)
% ------------------------------------------------------------------------
% This function unscrambles a scrambles curve using distances. The function
% assumes that the proper point order can be resolved by choosing the next
% nearest point.
%
% Change log:
% 2022/03/21 Updated description and documentation
% ------------------------------------------------------------------------

%%

%Deriving distance matrix
try
    D=dist(V,V');
catch
    D=distND(V,V);
end
D(eye(size(D))==1)=nan;

%Start loop to connect all points based on nearest distances
indFix=ones(1,size(V,1));
numPoints=size(V,1);
for qIter=2:1:numPoints
    [~,indMin]=gnanmin(D,[],2);
    indFix(qIter)=indMin(indFix(qIter-1));
    D(:,indFix)=NaN;
end
Vf=V(indFix,:);

% Collect output
varargout{1}=Vf;
varargout{2}=indFix;

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON