text:
function bool = inside_contour(pos, contour)

npos = size(pos,1);
ncnt = size(contour,1);
x = pos(:,1);
y = pos(:,2);

minx = min(x);
miny = min(y);
maxx = max(x);
maxy = max(y);

bool = true(npos,1);
bool(x<minx) = false;
bool(y<miny) = false;
bool(x>maxx) = false;
bool(y>maxy) = false;

% the summed angle over the contour is zero if the point is outside, and 2*pi if the point is inside the contour
% leave some room for inaccurate f
critval = 0.1;

% the remaining points have to be investigated with more attention
sel = find(bool);
for i=1:length(sel)
  contourx = contour(:,1) - pos(sel(i),1);
  contoury = contour(:,2) - pos(sel(i),2);
  angle = atan2(contoury, contourx);
  % angle = unwrap(angle);
  angle = my_unwrap(angle);
  total = sum(diff(angle));
  bool(sel(i)) = (abs(total)>critval);
end

function x = my_unwrap(x)
% this is a faster implementation of the MATLAB unwrap function
% with hopefully the same functionality
d    = diff(x);
indx = find(abs(d)>pi);
for i=indx(:)'
  if d(i)>0
    x((i+1):end) = x((i+1):end) - 2*pi;
  else
    x((i+1):end) = x((i+1):end) + 2*pi;
  end
end


text:
%--- help for ts/expanding ---
%
%  Applies a function to an expanding window of the time series
%
%  ::
%
%     h = expanding(db,func);
%     h = expanding(db,func,varargin);
%
%  Args:
%
%     db (ts object): time series object to get data
%
%     func: function that will apply to the expanding window of the time
%       series. The function is recursively applied to the data 1:t for
%       t=1,2,3,...,T, where T is the number of observations
%
%     varargin: additional arguments to func
%
%

author: jmaih


text:
 function X = dtft(x, omega, varargin)
%function X = dtft(x, omega [,options])
%|
%| Compute d-dimensional DTFT of signal x at frequency locations omega
%|
%| in
%|	x	[(Nd) L]	signal values
%|	omega	[M dd]		frequency locations (radians), dd = numel(Nd)
%|
%| option
%|	'n_shift' [dd 1]	use [0:N-1]-n_shift (default [0 .. 0])
%|	'how'			'outer' (default) big outer product
%|				'loop' reduce memory use (slower)
%|				'arrayfun' uses arrayfun()
%|
%| out
%|	X	[M L]		DTFT values
%|
%| Requires enough memory to store M * prod(Nd) size matrices (for testing)
%|
%| Copyright 2001-9-17, Jeff Fessler, University of Michigan
%| 2013-03-22, Daniel Weller added arrayfun version and other improvements
%| 2013-03-27, JF converted to vararg

if nargin == 1 && streq(x, 'test'), dtft_test(0), return, end
if nargin == 1 && streq(x, 'time'), dtft_test(1), return, end
if nargin < 2, ir_usage(), end

arg.n_shift = 0;
arg.how = 'outer';
arg = vararg_pair(arg, varargin);

dd = size(omega, 2);
Nd = size(x);

if numel(arg.n_shift) == 1
	arg.n_shift = repmat(arg.n_shift, dd);
end
n_shift = arg.n_shift;
if numel(n_shift) ~= dd
	fail 'n_shift size bad')
end

if dd == 1 && numel(Nd) == 2 && Nd(2) == 1 % 1D
	Nd = Nd(1);
end

if length(Nd) == dd % just one image
	x = x(:);
elseif length(Nd) == dd+1 % multiple images
	Nd = Nd(1:(end-1));
	x = reshapee(x, prod(Nd), []); % [*Nd L]
else
	error 'bad input signal size'
end

% dsw alternative to the loop:
% nn = arrayfun(@(nd,nshift) (0:(nd-1))-nshift, Nd, n_shift, 'UniformOutput', false);
for id=1:dd
	nn{id} = [0:(Nd(id)-1)] - n_shift(id);
end

% nn = ndgrid_jf('cell', nn);
if dd > 1
	[nn{:}] = ndgrid(nn{:});
end

switch arg.how
case 'outer'
	X = dtft_outer(x, omega, nn);
case 'loop'
	X = dtft_loop(x, omega, Nd, nn);
case 'arrayfun'
	X = dtft_arrayfun(x, omega, nn);
otherwise
	fail('unknown how "%s"', arg.how)
end

% dtft_outer()
function X = dtft_outer(x, omega, nn);
X = 0;
dd = ncol(omega);
for id=1:dd % add up phases
	X = X + omega(:,id) * col(nn{id})'; % [M *Nd]
end
X = exp(-1i*X) * x;

% dtft_loop()
% loop way: slower but less memory
function X = dtft_loop(x, omega, Nd, nn);
M = nrow(omega);
X = zeros(numel(x)/prod(Nd),M); % [L M]
%t1 = col(nn{1})';
if ncol(omega) > 3
	fail 'only up to 3d done'
end
if ncol(omega) < 3
	for dd = (ncol(omega)+1) : 3
		nn{dd} = 0; % dummy 0's
	end
	omega(1,3) = 0; % trick: make '3d'
end
t1 = nn{1}(:)';
t2 = col(nn{2})';
t3 = col(nn{3})';
for mm=1:M
	tmp = omega(mm,1)*t1 + omega(mm,2)*t2 + omega(mm,3)*t3;
	X(:,mm) = exp(-1i * tmp) * x;
end
X = X.'; % [M L]

% dtft_arrayfun()
% by Dan Weller, 2013-03-27
function X = dtft_arrayfun(x, omega, nn);
nn = cellfun(@(x) col(x).',nn,'UniformOutput',false); % make row vectors
nn = cat(1,nn{:}); % [dd *Nd]

M = nrow(omega);
X = arrayfun(@(m) exp((-1i*omega(m,:)) * nn) * x, 1:M, 'UniformOutput', false); % each cell is [1 L]
X = cat(1,X{:}); % [M L]

% X = exp(-1i*(omega * nn)) * x; % [M L]

% dtft_test
% simple test
function dtft_test(do_time)
Nd = [4 6 5] * 2^(1+do_time);
n_shift = [1 3 2];
rng(0), x = randn(Nd); % test signal
o1 = 2*pi*[0:(Nd(1)-1)]'/Nd(1); % test with uniform frequency locations
o2 = 2*pi*[0:(Nd(2)-1)]'/Nd(2);
o3 = 2*pi*[0:(Nd(3)-1)]'/Nd(3);
[o1 o2 o3] = ndgrid(o1, o2, o3);
om = [o1(:) o2(:) o3(:)];
cpu etic
Xd = dtft(x, om, 'n_shift', n_shift);
cpu etoc outer
cpu etic
Xl = dtft(x, om, 'n_shift', n_shift, 'how', 'loop');
cpu etoc loop
cpu etic
Xa = dtft(x, om, 'n_shift', n_shift, 'how', 'arrayfun');
cpu etoc arrayfun
printm('loop max %% difference = %g', max_percent_diff(Xl,Xd))
printm('afun max %% difference = %g', max_percent_diff(Xa,Xd))
Xf = fftn(x);
Xf = Xf(:) .* exp(1i * (om * n_shift(:))); % phase shift
printm('fftn max %% difference = %g', max_percent_diff(Xf,Xd))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function test_optimization_test03 ( )

%*****************************************************************************80
%
%% TEST_OPTIMIZATION_TEST03 tries Compass Search on each problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2012
%
%  Author:
%
%    John Burkardt
%
  m = 2;
  n = 1000;
  delta_tol = 0.000001;
  k_max = 20000;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_OPTIMIZATION_TEST03\n' );
  fprintf ( 1, '  For each problem, using dimension M = 2\n' );
  fprintf ( 1, '  try compass search.\n' );
%
%  Get the number of problems.
%
  problem_num = p00_problem_num ( );

  for problem = 1 : problem_num

    seed = 123456789;

    [ a, b ] = p00_ab ( problem, m );
    [ x0, seed ] = r8col_uniform ( m, 1, a, b, seed );
    fx = p00_f ( problem, m, 1, x0 );
    delta_init = 0.3 * norm ( x0 ) / m;
    delta_init = max ( delta_init, 1000.0 * delta_tol );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem %2d  DELTA_INIT = %14g\n', problem, delta_init );
    fprintf ( 1, '  Initial:  %14g  %14g  %14g\n', x0, fx );
    [ x, fx, k ] = p00_compass_search ( problem, m, x0, delta_tol, delta_init, ...
      k_max );
    fprintf ( 1, '  Final:    %14g  %14g  %14g  Steps = %d\n', x, fx, k );

    know = 0;
    while ( 1 )
      [ know, x ] = p00_sol ( problem, m, know );
      if ( know == 0 )
        break
      end
      fx = p00_f ( problem, m, 1, x );
      fprintf ( 1, '  Exact:    %14g  %14g  %14g\n', x, fx );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function D = addsmallobjectlabel(D, height, width, method)
%
% D = addsmallobjectlabel(D, height, width)
%
% Add the 'smallobject' label to any object that both dimensions are
% smaller than 32 pixels:
%
% D = addsmallobjectlabel(D, 32, 32)
%
% You can specify two decisions:
% D = addsmallobjectlabel(D, 32, 32, method)
%
% method = [and] | or
%
%
% Then we can use the LMquery to remove all the small objects.

if nargin<4
    method = 'and';
end

for i = 1:length(D)
    if isfield(D(i).annotation, 'object')
        Nobjects = length(D(i).annotation.object);
        bb = LMobjectboundingbox(D(i).annotation, 1:Nobjects)';

        for j = 1:Nobjects
            H = bb(4,:)-bb(2,:); % height of each annotated object
            W = bb(3,:)-bb(1,:); % width of each annotated object

            switch method
                case 'and'
                    small = find((H<height)&(W<width));
                case 'or'
                    small = find((H<height)|(W<width));
                otherwise
                    error('Invalid method')
            end

            for n = small
                D(i).annotation.object(n).name = [strtrim(D(i).annotation.object(n).name) ' smallobject'];
            end
        end
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
function [xm xcov] = UT(Xi, W, noiseCov)
%
%
[n, kmax] = size(Xi);

xm = 0;
for k=1:kmax
  xm = xm + W(k)*Xi(:, k);
end

xcov = zeros(n, n);
for k=1:kmax
  xcov = xcov + W(k)*(Xi(:, k) - xm)*(Xi(:, k) - xm)';
end
xcov = xcov + noiseCov;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function H=plot(x,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% KDE plotting function
%
%  plot(kde)            -- plot a KDE with various features.
%  plot(kde,style)      --  style is of the form
%  plot(kde,dim,style)  --    [STR1,STR2,...] where
%
%  style is of the form [STR1 STR2 STR3...] where
%
%  STR1 is the style to plot the line (1D) or kernel locations (2+D).  Note
%     that options must be specified in lowercase, e.g. 'ro' for red circles.
%     Default style is '-b'
%  STRN is a style for various optional plot features:
%    'W' : show kernel weights (by color: black = low, color in STR1 = high)
%    'S' : show relative kernel sizes, as circles around each center
%    'Bs': show KD-tree structure / bounding boxes, 's' is e.g. '-b' for
%             bounding boxes of solid blue lines
%    'Nd':show d levels of the bounding boxes (d an integer) (default all)
%  The default style is '.b'
%    Example styles:   '*kS-b'  -- black centers, blue variance circles
%                      '.gWS-kB-rN3' -- green dots, colored by weights,
%                          with black variances and 3 levels of bounding
%                          boxes in red.
%                      'S-k' -- just plot variances (no centers), in black
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

dim = []; args = [];
if (nargin == 1)
  dim = 1:getDim(x);
  if (getDim(x) == 1) args = 'b-'; else args = 'b.'; end;
elseif (nargin == 2)
  dim = [1:getDim(x)];
  args = varargin{1};
else
  dim = varargin{1};
  args = varargin{2};
end;
%
% Separate the arguments into:  [plot info][TAG taginfo][TAG taginfo] ...
%   where TAG is one of 'W','B','S'.
%
F = find(args ~= lower(args)); Fmin = min([F,length(args)+1]);
argsPlot = args(1:Fmin-1);
argsKDE = args(Fmin:length(args));

if (getDim(x) == 1)
  pts = getPoints(x);
  N = 200;   range = [min(pts),max(pts)];
  range(1) = range(1) - .05*(range(2)-range(1));
  range(2) = range(2) + .05*(range(2)-range(1));
  H=draw1D(x,linspace(range(1),range(2),N),argsPlot,argsKDE);
else
  H=drawAllPairs(x,dim,argsPlot,argsKDE);
end;

%
% Internal functions
%
function e=draw1D(x, bins, style, myStyle)
  y = evaluate(x,bins);
  e=plot(bins,y,style);

  mx = max(y);
  wts = getWeights(x);

  if(strfind(myStyle,'W'))
    holdf = ishold; hold on;
    subStyle = extract(myStyle,'W');
    subStyle=subStyle(2:end); if (length(subStyle)==0) subStyle = '^'; end;
    etmp = stem(getPoints(x), wts, subStyle);
    e = [e;etmp];
    if (~holdf) hold off; end;
  end

  if(strfind(myStyle,'S'))
    holdf = ishold; hold on;
    subStyle = extract(myStyle,'S');
    subStyle=subStyle(2:end); if (length(subStyle)==0) subStyle = '--b'; end;

    bw = getBW(x); pts = getPoints(x); type = getType(x);
    for i=1:length(pts)
      xtmp = kde(pts(i), bw(i), 1, type);
      etmp = plot(bins, evaluate(xtmp, bins) * wts(i), subStyle);
      e = [e;etmp];
    end

    % for plotting only gaussian kernels faster:
%    gr = repmat(bins, length(wts), 1);
%    wts = repmat(wts', 1, size(gr, 2));
%    pts = repmat(getPoints(x)', 1, size(gr,2));
%    bw = repmat(getBW(x)', 1, size(gr, 2));
%    etmp = plot( bins, (wts.*(1./(sqrt(2*pi)*bw)).*exp(-(gr-pts).^2./(2*bw.^2)))', subStyle);

    e = [e;etmp];
    if (~holdf) hold off; end;
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function e=drawAllPairs(x,dims,style,myStyle)
  pts = getPoints(x);
  holdf = ishold;
  e = [];
  Nout = length(dims);
  PlotI2 = triu(repmat(dims ,[Nout,1]), 1);
  PlotI1 = triu(repmat((dims)',[1,Nout]), 1);
  PlotI1 = PlotI1(find(PlotI1)); PlotI2 = PlotI2(find(PlotI2));
  Ncol = fix(sqrt(length(PlotI2))); Nrow = ceil(length(PlotI2)/Ncol);
  for iT=1:length(PlotI2)                   %  output all dimension pairs:
    subplot(Nrow,Ncol,iT);
    holdfs = ishold;
    if (holdf) hold on; end;
    etmp = drawPair(x,[PlotI1(iT),PlotI2(iT)],style,myStyle);
    if (~holdfs) hold off; end;
    e = [e;etmp];
  end;
  drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function etmp = drawPair(x,dims,style,myStyle)
  pts = getPoints(x);
  etmp = [];

  if (strfind(myStyle,'W'))  %plot weight info with color
    subStyle = extract(myStyle,'W');
    %wts = getWeights(x);  wts = wts - min(wts);
    %if (max(wts)==0) wts = wts+1; else wts = .75*wts/max(wts) + .25; end;
    wts = getWeights(x);
    wts = .75*wts/max(wts) + .25;
    holdf = ishold;
    for i=1:size(pts,2)
      etmp2 = plot(pts(dims(1),i),pts(dims(2),i),style);   % plot each location and
      hold on;
      etmp = [etmp;etmp2]; ctmp = get(etmp2,'Color');      %  its color in prop.
      set(etmp2,'Color',[1 1 1] - wts(i)*([1 1 1] - ctmp));%  to its weight
    end;
    if (~holdf) hold off; end;
  else                                                    % otherwise, just plot
    if (length(style)~=0)
      etmp = plot(pts(dims(1),:),pts(dims(2),:),style);   %   all locations
  end; end;

  if (strfind(myStyle,'S'))  %plot BW info with circles
    subStyle = extract(myStyle,'S');
    subStyle=subStyle(2:end); if (length(subStyle)==0) subStyle = '-b'; end;
    pts = getPoints(x); sig=getBW(x);
	meanX = pts(dims(1),:);
	meanY = pts(dims(2),:);
	sigX = sig(dims(1),:);
	sigY = sig(dims(2),:);

	holdf = ishold;	hold on;
	theta = linspace(0,2*pi,100);
	for (i = 1:length(meanX))
      etmp2 = plot(meanX(i)+sigX(i)*cos(theta), meanY(i)+sigY(i)*sin(theta),subStyle);
      etmp = [etmp;etmp2];
	end
	if (~holdf) hold off; end;
  end;

  if (strfind(myStyle,'B'))  %plot bounding box info with rectangles
    levels = [];
    if (isempty(strfind(myStyle,'N'))) levels = ceil(log2(getNpts(x)));
    else
      subStyle = extract(myStyle,'N');              % get # of balls if spec'd
      levels = sscanf(subStyle(2:end),'%d');
      levels = min(ceil(log2(getNpts(x))), levels);
    end;
    subStyle = extract(myStyle,'B');              % get plot style if spec'd
    subStyle=subStyle(2:end);
    if(isempty(subStyle))     subStyle = '-b';    end

    N = getNpts(x);
    indices = []; tmp = [1];
    for i=1:levels
      indices = [indices tmp];
      % get rid of NO_CHILD right children
      rt = double(x.rightch(tmp)) + 1;
      rt = rt .* (rt < N+1) + 1 * (rt > N);
      tmp = [double(x.leftch(tmp))+1   rt];
    end

    rX = x.ranges(dims(1),indices);
    rY = x.ranges(dims(2),indices);
    nodesX = x.centers(dims(1),indices);
    nodesY = x.centers(dims(2),indices);
    leavesX = x.centers(dims(1),N+1:end);
    leavesY = x.centers(dims(2),N+1:end);

    squaresX(1,:) = nodesX + rX;
    squaresX(2,:) = nodesX + rX;
    squaresX(3,:) = nodesX - rX;
    squaresX(4,:) = nodesX - rX;
    squaresX(5,:) = nodesX + rX;

    squaresY(1,:) = nodesY + rY;
    squaresY(2,:) = nodesY - rY;
    squaresY(3,:) = nodesY - rY;
    squaresY(4,:) = nodesY + rY;
    squaresY(5,:) = nodesY + rY;

    holdf = ishold;	hold on;
    for i = 1:length(indices)
      etmp2 = plot(squaresX, squaresY, subStyle);
      etmp = [etmp;etmp2];
    end
%   axis square
    if (~holdf) hold off; end;

  end;

  titlestr = ['Dim ',int2str(dims(1)),' v. ',int2str(dims(2))];
  title(titlestr);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function substr=extract(str,tag)
  substr = [];
  loc = strfind(str,tag);
  if (loc)
    substr = str(loc(1):length(str));
    loc = find(substr(2:end) ~= lower(substr(2:end)));
    if (loc) substr = substr(1:loc(1)); end;
  end;

author: ShapeNet

repo: RenderForCNN


text:
function [SDR,ISR,SIR,SAR,perm]=bss_eval_images(ie,i)

% BSS_EVAL_IMAGES Ordering and measurement of the separation quality for
% estimated source spatial image signals in terms of true source, spatial
% (or filtering) distortion, interference and artifacts.
%
% [SDR,ISR,SIR,SAR,perm]=bss_eval_images(ie,i)
%
% Inputs:
% ie: nsrc x nsampl x nchan matrix containing estimated source images
% i: nsrc x nsampl x nchan matrix containing true source images
%
% Outputs:
% SDR: nsrc x 1 vector of Signal to Distortion Ratios
% ISR: nsrc x 1 vector of source Image to Spatial distortion Ratios
% SIR: nsrc x 1 vector of Source to Interference Ratios
% SAR: nsrc x 1 vector of Sources to Artifacts Ratios
% perm: nsrc x 1 vector containing the best ordering of estimated source
% images in the mean SIR sense (estimated source image number perm(j)
% corresponds to true source image number j)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2007-2008 Emmanuel Vincent
% This software is distributed under the terms of the GNU Public License
% version 3 (http://www.gnu.org/licenses/gpl.txt)
% If you find it useful, please cite the following reference:
% Emmanuel Vincent, Hiroshi Sawada, Pau Bofill, Shoji Makino and Justinian
% P. Rosca, "First stereo audio source separation evaluation campaign:
% data, algorithms and results," In Proc. Int. Conf. on Independent
% Component Analysis and Blind Source Separation (ICA), pp. 552-559, 2007.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Errors %%%
if nargin<2, error('Not enough input arguments.'); end
[nsrc,nsampl,nchan]=size(ie);
[nsrc2,nsampl2,nchan2]=size(i);
if nsrc2~=nsrc, error('The number of estimated source images and reference source images must be equal.'); end
if nsampl2~=nsampl, error('The estimated source images and reference source images must have the same duration.'); end
if nchan2~=nchan, error('The estimated source images and reference source images must have the same number of channels.'); end

%%% Performance criteria %%%
% Computation of the criteria for all possible pair matches
SDR=zeros(nsrc,nsrc);
ISR=zeros(nsrc,nsrc);
SIR=zeros(nsrc,nsrc);
SAR=zeros(nsrc,nsrc);
for jest=1:nsrc,
    for jtrue=1:nsrc,
        [s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(reshape(ie(jest,:,:),nsampl,nchan).',i,jtrue,512);
        [SDR(jest,jtrue),ISR(jest,jtrue),SIR(jest,jtrue),SAR(jest,jtrue)]=bss_image_crit(s_true,e_spat,e_interf,e_artif);
    end
end
% Selection of the best ordering
perm=perms(1:nsrc);
nperm=size(perm,1);
meanSIR=zeros(nperm,1);
for p=1:nperm,
    meanSIR(p)=mean(SIR((0:nsrc-1)*nsrc+perm(p,:)));
end
[meanSIR,popt]=max(meanSIR);
perm=perm(popt,:).';
SDR=SDR((0:nsrc-1).'*nsrc+perm);
ISR=ISR((0:nsrc-1).'*nsrc+perm);
SIR=SIR((0:nsrc-1).'*nsrc+perm);
SAR=SAR((0:nsrc-1).'*nsrc+perm);

return;

function [s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se,s,j,flen)

% BSS_DECOMP_MTIFILT Decomposition of an estimated source image into four
% components representing respectively the true source image, spatial (or
% filtering) distortion, interference and artifacts, derived from the true
% source images using multichannel time-invariant filters.
%
% [s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se,s,j,flen)
%
% Inputs:
% se: nchan x nsampl matrix containing the estimated source image (one row per channel)
% s: nsrc x nsampl x nchan matrix containing the true source images
% j: source index corresponding to the estimated source image in s
% flen: length of the multichannel time-invariant filters in samples
%
% Outputs:
% s_true: nchan x nsampl matrix containing the true source image (one row per channel)
% e_spat: nchan x nsampl matrix containing the spatial (or filtering) distortion component
% e_interf: nchan x nsampl matrix containing the interference component
% e_artif: nchan x nsampl matrix containing the artifacts component

%%% Errors %%%
if nargin<4, error('Not enough input arguments.'); end
[nchan2,nsampl2]=size(se);
[nsrc,nsampl,nchan]=size(s);
if nchan2~=nchan, error('The number of channels of the true source images and the estimated source image must be equal.'); end
if nsampl2~=nsampl, error('The duration of the true source images and the estimated source image must be equal.'); end

%%% Decomposition %%%
% True source image
s_true=[reshape(s(j,:,:),nsampl,nchan).',zeros(nchan,flen-1)];
% Spatial (or filtering) distortion
e_spat=project(se,s(j,:,:),flen)-s_true;
% Interference
e_interf=project(se,s,flen)-s_true-e_spat;
% Artifacts
e_artif=[se,zeros(nchan,flen-1)]-s_true-e_spat-e_interf;

return;

function sproj=project(se,s,flen)

% SPROJ Least-squares projection of each channel of se on the subspace
% spanned by delayed versions of the channels of s, with delays between 0
% and flen-1

[nsrc,nsampl,nchan]=size(s);
s=reshape(permute(s,[3 1 2]),nchan*nsrc,nsampl);

%%% Computing coefficients of least squares problem via FFT %%%
% Zero padding and FFT of input data
s=[s,zeros(nchan*nsrc,flen-1)];
se=[se,zeros(nchan,flen-1)];
fftlen=2^nextpow2(nsampl+flen-1);
sf=fft(s,fftlen,2);
sef=fft(se,fftlen,2);
% Inner products between delayed versions of s
G=zeros(nchan*nsrc*flen);
for k1=0:nchan*nsrc-1,
    for k2=0:k1,
        ssf=sf(k1+1,:).*conj(sf(k2+1,:));
        ssf=real(ifft(ssf));
        ss=toeplitz(ssf([1 fftlen:-1:fftlen-flen+2]),ssf(1:flen));
        G(k1*flen+1:k1*flen+flen,k2*flen+1:k2*flen+flen)=ss;
        G(k2*flen+1:k2*flen+flen,k1*flen+1:k1*flen+flen)=ss.';
    end
end
% Inner products between se and delayed versions of s
D=zeros(nchan*nsrc*flen,nchan);
for k=0:nchan*nsrc-1,
    for i=1:nchan,
        ssef=sf(k+1,:).*conj(sef(i,:));
        ssef=real(ifft(ssef,[],2));
        D(k*flen+1:k*flen+flen,i)=ssef(:,[1 fftlen:-1:fftlen-flen+2]).';
    end
end

%%% Computing projection %%%
% Distortion filters
C=G\D;
C=reshape(C,flen,nchan*nsrc,nchan);
% Filtering
sproj=zeros(nchan,nsampl+flen-1);
for k=1:nchan*nsrc,
    for i=1:nchan,
        sproj(i,:)=sproj(i,:)+fftfilt(C(:,k,i).',s(k,:));
    end
end

return;

function [SDR,ISR,SIR,SAR]=bss_image_crit(s_true,e_spat,e_interf,e_artif)

% BSS_IMAGE_CRIT Measurement of the separation quality for a given source
% image in terms of true source, spatial (or filtering) distortion,
% interference and artifacts.
%
% [SDR,ISR,SIR,SAR]=bss_image_crit(s_true,e_spat,e_interf,e_artif)
%
% Inputs:
% s_true: nchan x nsampl matrix containing the true source image (one row per channel)
% e_spat: nchan x nsampl matrix containing the spatial (or filtering) distortion component
% e_interf: nchan x nsampl matrix containing the interference component
% e_artif: nchan x nsampl matrix containing the artifacts component
%
% Outputs:
% SDR: Signal to Distortion Ratio
% ISR: source Image to Spatial distortion Ratio
% SIR: Source to Interference Ratio
% SAR: Sources to Artifacts Ratio

%%% Errors %%%
if nargin<4, error('Not enough input arguments.'); end
[nchant,nsamplt]=size(s_true);
[nchans,nsampls]=size(e_spat);
[nchani,nsampli]=size(e_interf);
[nchana,nsampla]=size(e_artif);
if ~((nchant==nchans)&&(nchant==nchani)&&(nchant==nchana)), error('All the components must have the same number of channels.'); end
if ~((nsamplt==nsampls)&&(nsamplt==nsampli)&&(nsamplt==nsampla)), error('All the components must have the same duration.'); end

%%% Energy ratios %%%
% SDR
SDR=10*log10(sum(sum(s_true.^2))/sum(sum((e_spat+e_interf+e_artif).^2)));
% ISR
ISR=10*log10(sum(sum(s_true.^2))/sum(sum(e_spat.^2)));
% SIR
SIR=10*log10(sum(sum((s_true+e_spat).^2))/sum(sum(e_interf.^2)));
% SAR
SAR=10*log10(sum(sum((s_true+e_spat+e_interf).^2))/sum(sum(e_artif.^2)));

return;

author: posenhuang

repo: deeplearningsourceseparation


text:
function borderCoordsRoi=dtiFindBorderBetweenRois(roi1_img, roi2_img, dt6File, minDist, borderRoiName)
%Find a shared border separating two ROIs
%
% borderCoords=dtiFindBorderBetweenRois(roi1_img, roi2_img, dt6File, [minDist=1.87])
%
% Given two nifti images with ROIs, find a set of points where the two
% ROIs are touching within the distance set by minDist parameter.
%
% Example:
%        dt6File=fullfile(pathtoChild, 'dt6.mat');
%        roi1_img=fullfile(pathtoChildInFreesurfer, 'aparc_aseg1015.nii');
%        roi2_img=fullfile(pathtoChildInFreesurfer, 'aparc_aseg1030.nii');
%        borderCoordsRoi=dtiFindBorderBetweenRois(roi1_img, roi2_img, dt6File, 5, borderRoiName);
%        borderCoordsRoiFile=fullfile(pathtoChild, 'ROIs', 'border1015to1030');
%        borderRoiName='borderRoi1toRoi2';
%        dtiWriteRoi(borderCoordsRoi, borderCoordsRoiFile);
%
% (c) Vistalab

% HISTORY:
% 03/2010 ER wrote it

if ~exist('minDist', 'var') || isempty(minDist)
minDist=1.87; %mm
end

roi1=dtiImportRoiFromNifti(roi1_img, dt6File);
roi2=dtiImportRoiFromNifti(roi2_img, dt6File);

[indices, bestSq]=nearpoints(roi1.coords', roi2.coords');
borderCoords=(roi1.coords(bestSq<(minDist^2), :)+roi2.coords(indices(bestSq<(minDist^2)), :))./2;
borderCoordsRoi=dtiNewRoi(borderRoiName, [], borderCoords);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ok = all_set_item(item)

% function ok = all_set_item(item)
% Perform within-item all_set check. For branches, this is always true.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: all_set_item.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok
ok = true;

author: spm


text:
function p02_demo ( iteration_max, h )

%*****************************************************************************80
%
%% P02_DEMO runs the 2D demo problem #2, with mesh size H.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    06 February 2006
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer ITERATION_MAX, the maximum number of iterations that DISTMESH
%    should take.  (The program might take fewer iterations if it detects convergence.)
%
%    Input, real H, the mesh spacing parameter.
%
  if ( nargin < 1 )
    iteration_max = 200;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P02_DEMO - Note:\n' );
    fprintf ( 1, '  No value of ITERATION_MAX was supplied.\n' );
    fprintf ( 1, '  The default value ITERATION_MAX = %d will be used.\n', ...
      iteration_max );
  end

  if ( nargin < 2 )
    h = 0.10;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P02_DEMO - Note:\n' );
    fprintf ( 1, '  No value of H was supplied.\n' );
    fprintf ( 1, '  The default value H = %f will be used.\n', h );
  end
%
%  Put the random number generator into a fixed initial state.
%
  rand ( 'state', 111 );
%
%  Set the rendering method for the current figure to Z-buffering.
%
  set ( gcf, 'rend', 'z' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 2:\n' );
  fprintf ( 1, '  Unit circle with a hole, h = %f\n', h )

  fd = @p02_fd;
  fh = @p02_fh;
  box = [-1.0, -1.0; 1.0, 1.0 ];
  fixed = [];

  [ p, t ] = distmesh_2d ( fd, fh, h, box, iteration_max, fixed );

  post_2d ( p, t, fh )
%
%  Write a PostScript image of the triangulation.
%
  [ node_num, junk ] = size ( p );
  [ tri_num, junk ] = size ( t );
  p = p';
  t = t';
  node_show = 0;
  triangle_show = 1;

  triangulation_order3_plot ( 'p02_mesh.eps', node_num, p, tri_num, ...
    t, node_show, triangle_show );
%
%  Write a text file containing the nodes.
%
  r8mat_write ( 'p02_nodes.txt', 2, node_num, p );
%
%  Write a text file containing the triangles.
%
  i4mat_write ( 'p02_elements.txt', 3, tri_num, t );

  return
end


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters all_theta
%   p = PREDICT(all_theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(all_theta'*x) >= 0.5, predict 1)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

% train each row for all classes and keep the label with the highest probability
htheta = X * all_theta';
[temp, p] = max(htheta, [], 2);

% =========================================================================

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
function OUTPUT = swap_spaces(STRING, varargin)
%replaces spaces with underscores and vice versa
%SWAP_SPACES replaces spaces with underscores or underscores with spaces
% for desciptions and file names. If there are underscores AND spaces
% present in STRING, either give an option or only underscores will be replaced.
%
%Syntax
%    swap_spaces(textstring [, option])
%
%    Valid options are:
%     * '_' replaces all underscores with spaces
%     * ' ' replaces all spaces with underscores
%
%Example
%    f = swap_spaces('text_text')
%    f =
%        'text text'
%
%    g = swap_spaces('text text')
%    g =
%        'text_text'
%
%    h = swap_spaces('text text_text', '_')
%    h =
%        'text text text'
%See also
%    regexprep, strrep

p	=	inputParser;
validOptions = {'_', ' '};
p.addRequired('STRING',			@(x) all(ischar(x), ~isempty(x)) );
p.addOptional('OPT',	'',		@(x) any(strcmp(x, validOptions)) );
p.parse(STRING, varargin{:});

if isempty(p.Results.OPT)
	if strfind(STRING, '_')
		OUTPUT	=	strrep(STRING, '_', ' ');
	else
		OUTPUT	=	strrep(STRING, ' ', '_');
	end
elseif strcmp(p.Results.OPT, '_')
	OUTPUT	=	strrep(STRING, '_', ' ');
else
	OUTPUT	=	strrep(STRING, ' ', '_');
end

% Copyright 2009-2011 Alexandra Heidsieck <aheidsieck@tum.de>,
%                     IMETUM, Technische Universitaet Muenchen

author: Sable

repo: mcbench-benchmarks


text:
function h=plotEllipse(ra,rb,ang,x0,y0,C,Nb)
% Ellipse adds ellipses to the current plot
%
% ELLIPSE(ra,rb,ang,x0,y0) adds an ellipse with semimajor axis of ra,
% a semimajor axis of radius rb, a semimajor axis of ang, centered at
% the point x0,y0.
%
% The length of ra, rb, and ang should be the same.
% If ra is a vector of length L and x0,y0 scalars, L ellipses
% are added at point x0,y0.
% If ra is a scalar and x0,y0 vectors of length M, M ellipse are with the same
% radii are added at the points x0,y0.
% If ra, x0, y0 are vectors of the same length L=M, M ellipses are added.
% If ra is a vector of length L and x0, y0 are  vectors of length
% M~=L, L*M ellipses are added, at each point x0,y0, L ellipses of radius ra.
%
% ELLIPSE(ra,rb,ang,x0,y0,C)
% adds ellipses of color C. C may be a string ('r','b',...) or the RGB value.
% If no color is specified, it makes automatic use of the colors specified by
% the axes ColorOrder property. For several circles C may be a vector.
%
% ELLIPSE(ra,rb,ang,x0,y0,C,Nb), Nb specifies the number of points
% used to draw the ellipse. The default value is 300. Nb may be used
% for each ellipse individually.
%
% h=ELLIPSE(...) returns the handles to the ellipses.
%
% as a sample of how ellipse works, the following produces a red ellipse
% tipped up at a 45 deg axis from the x axis
% ellipse(1,2,pi/8,1,1,'r')
%
% note that if ra=rb, ELLIPSE plots a circle
%

% written by D.G. Long, Brigham Young University, based on the
% CIRCLES.m original
% written by Peter Blattner, Institute of Microtechnology, University of
% Neuchatel, Switzerland, blattner@imt.unine.ch

% Check the number of input arguments

if nargin<1,
  ra=[];
end;
if nargin<2,
  rb=[];
end;
if nargin<3,
  ang=[];
end;

%if nargin==1,
%  error('Not enough arguments');
%end;

if nargin<5,
  x0=[];
  y0=[];
end;

if nargin<6,
  C=[];
end

if nargin<7,
  Nb=[];
end

% set up the default values

if isempty(ra),ra=1;end;
if isempty(rb),rb=1;end;
if isempty(ang),ang=0;end;
if isempty(x0),x0=0;end;
if isempty(y0),y0=0;end;
if isempty(Nb),Nb=300;end;
if isempty(C),C=get(gca,'colororder');end;

% work on the variable sizes

x0=x0(:);
y0=y0(:);
ra=ra(:);
rb=rb(:);
ang=ang(:);
Nb=Nb(:);

if isstr(C),C=C(:);end;

if length(ra)~=length(rb),
  error('length(ra)~=length(rb)');
end;
if length(x0)~=length(y0),
  error('length(x0)~=length(y0)');
end;

% how many inscribed elllipses are plotted

if length(ra)~=length(x0)
  maxk=length(ra)*length(x0);
else
  maxk=length(ra);
end;

% drawing loop

for k=1:maxk

  if length(x0)==1
    xpos=x0;
    ypos=y0;
    radm=ra(k);
    radn=rb(k);
    if length(ang)==1
      an=ang;
    else
      an=ang(k);
    end;
  elseif length(ra)==1
    xpos=x0(k);
    ypos=y0(k);
    radm=ra;
    radn=rb;
    an=ang;
  elseif length(x0)==length(ra)
    xpos=x0(k);
    ypos=y0(k);
    radm=ra(k);
    radn=rb(k);
    an=ang(k)
  else
    rada=ra(fix((k-1)/size(x0,1))+1);
    radb=rb(fix((k-1)/size(x0,1))+1);
    an=ang(fix((k-1)/size(x0,1))+1);
    xpos=x0(rem(k-1,size(x0,1))+1);
    ypos=y0(rem(k-1,size(y0,1))+1);
  end;

  co=cos(an);
  si=sin(an);
  the=linspace(0,2*pi,Nb(rem(k-1,size(Nb,1))+1,:)+1);
%  x=radm*cos(the)*co-si*radn*sin(the)+xpos;
%  y=radm*cos(the)*si+co*radn*sin(the)+ypos;
  h(k)=line(radm*cos(the)*co-si*radn*sin(the)+xpos,radm*cos(the)*si+co*radn*sin(the)+ypos);
  set(h(k),'color',C(rem(k-1,size(C,1))+1,:));

end;

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function sigmergt
%SIGMERGT Unit test for the function SIGMERGE.

%	O. Lemoine - March 1996.

N=256;

s1=noisecg(N);
s2=amgauss(N);
R=10;
sig=sigmerge(s1,s2,R);
E1=norm(s1)^2;
E2=norm(s2)^2;
h=sqrt(E1/(E2*10^(R/10)));
err=sig-s1-h*s2;
if any(abs(err)>sqrt(eps)),
 error('sigmerge test 1 failed');
end

s1=noisecu(N);
s2=fmlin(N);
R=-2.4;
sig=sigmerge(s1,s2,R);
E1=norm(s1)^2;
E2=norm(s2)^2;
h=sqrt(E1/(E2*10^(R/10)));
err=sig-s1-h*s2;
if any(abs(err)>sqrt(eps)),
 error('sigmerge test 2 failed');
end

N=245;

s1=noisecg(N);
s2=amgauss(N);
R=10;
sig=sigmerge(s1,s2,R);
E1=norm(s1)^2;
E2=norm(s2)^2;
h=sqrt(E1/(E2*10^(R/10)));
err=sig-s1-h*s2;
if any(abs(err)>sqrt(eps)),
 error('sigmerge test 3 failed');
end

s1=noisecu(N);
s2=fmlin(N);
R=-2.4;
sig=sigmerge(s1,s2,R);
E1=norm(s1)^2;
E2=norm(s2)^2;
h=sqrt(E1/(E2*10^(R/10)));
err=sig-s1-h*s2;
if any(abs(err)>sqrt(eps)),
 error('sigmerge test 4 failed');
end

author: HeLiangHIT


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   [JX, Jq, Jphi]=compute_jacobians_5R(robot, q)
%
%   Compute Jacobians for the for the 5R planar parallel robot
%
%   JX is the 4x2 Jacobian with respect to the cartesian coordinates (x, y)
%   Jq is the 4x2 Jacobian with respect to the active joints (q1, q2)
%   Jphi is the 4x2 Jacobian with respect to the passive joints (phi1, phi2)
%
%   Author: Arturo Gil Aparicio arturo.gil@umh.es
%   Date: 13/09/2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 [JX, Jq, Jphi]=compute_jacobians_5R(robot, q)

% The geometric parameters of this robot can be retrieved from
% the two DOF planar arms that form it
a1=eval(robot.robot1.DH.a);
a2=eval(robot.robot2.DH.a);

%Link lengths
l1=abs(a1(1));
l2=abs(a1(2));
l3=abs(a2(1));
l4=abs(a2(2));
%distance along X that separate both arms
L=robot.L;

%retrieve the active and passive joints from the vector q
q1=q(1);
phi1=q(2);
q2=q(3);
phi2=q(4);

%The loop closure equations for this robot can be written as:
%
% \begin{bmatrix}
% \Gamma_1\\
% \Gamma_2\\
% \Gamma_3\\
% \Gamma_4
% \end{bmatrix}=
% \begin{bmatrix}
% x-l_1\cos(q_1)-l_2\cos(q_1+\varphi_1)\\
% y-l_1\sin(q_1)-l_2\sin(q_1+\varphi_1)\\
% x-l_4\cos(q_2)-l_3\cos(q_2+\varphi_2)-L\\
% y-l_4\sin(q_2)-l_3\sin(q_2+\varphi_2)
% \end{bmatrix}

%In this matrix, JX11 is the partial derivative of Gamma1 with respect to x
%                JX12 is the partial derivative of Gamma1 with respect to y
%                JX21 is the partial derivative of Gamma2 with respect to x
%...
JX=[1   0;
    0   1;
    1   0;
    0   1];

%In this matrix, Jq11 is the partial der. of Gamma1 with respect to q1
%                Jq12 is the partial der. of Gamma1 with respect to q2
%                Jq21 is the partial der. of Gamma2 with respect to q1
%etc...
Jq = [l1*sin(q1)+l2*sin(q1+phi1)              0;
      -l1*cos(q1)-l2*cos(q1+phi1)             0;
      0                           l4*sin(q2)+l3*sin(q2+phi2);
      0                           -l4*cos(q2)-l3*cos(q2+phi2)];

%In this matrix, Jphi11 is the partial der. of Gamma1 with respect to phi1
%                Jphi12 is the partial der. of Gamma1 with respect to phi2
%                Jphi21 is the partial der. of Gamma2 with respect to phi1
%etc...
Jphi = [l2*sin(q1+phi1)              0;
       -l2*cos(q1+phi1)             0;
      0                           l3*sin(q2+phi2);
      0                           -l3*cos(q2+phi2)];

repo: ARTE


text:
function PopDec = EvolALG(Problem,PCheby,Dec,model,IFEs)
% Solution update in ParEGO, where a solution with the best expected
% improvement is re-evaluated

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    Off   = [OperatorGA(Problem,Dec(TournamentSelection(2,size(Dec,1),PCheby),:));OperatorGA(Problem,Dec,{0,0,1,20})];
    N     = size(Off,1);
    EI    = zeros(N,1);
    Gbest = min(PCheby);
    E0    = inf;
    while IFEs > 0
        drawnow('limitrate');
        for i = 1 : N
            [y,~,mse] = predictor(Off(i,:),model);
            s         = sqrt(mse);
            EI(i)     = -(Gbest-y)*normcdf((Gbest-y)/s)-s*normpdf((Gbest-y)/s);
        end
        [~,index] = sort(EI);
        if EI(index(1)) < E0
            Best = Off(index(1),:);
            E0   = EI(index(1));
        end
        Parent = Off(index(1:ceil(N/2)),:);
        Off    = [OperatorGA(Problem,Parent(TournamentSelection(2,size(Parent,1),EI(index(1:ceil(N/2)))),:));OperatorGA(Problem,Parent,{0,0,1,20})];
        IFEs   = IFEs - size(Off,1);
    end
    PopDec = Best;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [top, top_confidences] = nms(boxes, confidences, overlap)

% top = nms(boxes, overlap)
% Non-maximum suppression.
% Greedily select high-scoring detections and skip detections
% that are significantly covered by a previously selected detection.

if isempty(boxes)
  top = [];
else
  x1 = boxes(:,1);
  y1 = boxes(:,2);
  x2 = boxes(:,3);
  y2 = boxes(:,4);
  s = confidences';
  area = (x2-x1+1) .* (y2-y1+1);

  [vals, I] = sort(s);
  pick = [];
  while ~isempty(I)
    last = length(I);
    i = I(last);
    pick = [pick; i];
    suppress = [last];
    for pos = 1:last-1
      j = I(pos);
      xx1 = max(x1(i), x1(j));
      yy1 = max(y1(i), y1(j));
      xx2 = min(x2(i), x2(j));
      yy2 = min(y2(i), y2(j));
      w = xx2-xx1+1;
      h = yy2-yy1+1;
      if w > 0 && h > 0
        % compute overlap
        o = w * h / area(j);
        if o > overlap
          suppress = [suppress; pos];
        end
      end
    end
    I(suppress) = [];
  end
  top = boxes(pick,:);
  top_confidences = confidences(pick);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [xtable, ytable, utable, vtable, typevector,correlation_map] = piv_FFTensemble (autolimit,filepath,video_frame_selection,bg_img_A,bg_img_B,clahe,highp,intenscap,clahesize,highpsize,wienerwurst,wienerwurstsize,roi_inpt,maskiererx,maskierery,interrogationarea,step,subpixfinder,passes,int2,int3,int4,mask_auto,imdeform,repeat,do_pad)
%this funtion performs the  PIV analysis. It is a modification of the
%pivFFTmulti, and will do ensemble correlation. That is a suitable
%algorithm for low seeding density as it happens in microPIV.
warning off %#ok<*WNOFF> %MATLAB:log:logOfZero
%% pre-processing is done in this function
result_conv_ensemble = zeros(interrogationarea,interrogationarea); % prepare empty result_conv
if isempty(video_frame_selection) %list with image files was passed
	amount_input_imgs=size(filepath,1);
else
	amount_input_imgs=numel(video_frame_selection);
end
total_analyses_amount=amount_input_imgs / 2 * passes;
from_total = 0;
tic
skippy=0;
for ensemble_i1=1:2:amount_input_imgs
	if isempty(video_frame_selection) %list with image files was passed
		%detect if it is b16 or standard pixel image
		[~,~,ext] = fileparts(filepath{1});
		if strcmp(ext,'.b16')
			image1=f_readB16(filepath{ensemble_i1});
			image2=f_readB16(filepath{ensemble_i1+1});
		else
			image1=imread(filepath{ensemble_i1});
			image2=imread(filepath{ensemble_i1+1});
		end
	else % video file was passed
		image1 = read(filepath,video_frame_selection(ensemble_i1));
		image2 = read(filepath,video_frame_selection(ensemble_i1+1));
	end
	if size(image1,3)>1
		image1=uint8(mean(image1,3));
		image2=uint8(mean(image2,3));
		%disp('Warning: To optimize speed, your images should be grayscale, 8 bit!')
	end
	%Subtract background (if existent)
	if ~isempty(bg_img_A)
		image1=image1-bg_img_A;
	end
	if ~isempty(bg_img_B)
		image2=image2-bg_img_B;
	end
	%if autolimit == 1 %if autolimit is desired: do autolimit for each image seperately
		if size(image1,3)>1
			stretcher = stretchlim(rgb2gray(image1));
		else
			stretcher = stretchlim(image1);
		end
		minintens1 = stretcher(1);
		maxintens1 = stretcher(2);
		if size(image2,3)>1
			stretcher = stretchlim(rgb2gray(image2));
		else
			stretcher = stretchlim(image2);
		end
		minintens2 = stretcher(1);
		maxintens2 = stretcher(2);
	%end
	image1 = PIVlab_preproc (image1,roi_inpt,clahe, clahesize,highp,highpsize,intenscap,wienerwurst,wienerwurstsize,minintens1,maxintens1);
	image2 = PIVlab_preproc (image2,roi_inpt,clahe, clahesize,highp,highpsize,intenscap,wienerwurst,wienerwurstsize,minintens2,maxintens2);
	if numel(roi_inpt)>0
		xroi=roi_inpt(1);
		yroi=roi_inpt(2);
		widthroi=roi_inpt(3);
		heightroi=roi_inpt(4);
		image1_roi=double(image1(yroi:yroi+heightroi,xroi:xroi+widthroi));
		image2_roi=double(image2(yroi:yroi+heightroi,xroi:xroi+widthroi));
	else
		xroi=0;
		yroi=0;
		image1_roi=double(image1);
		image2_roi=double(image2);
	end
	gen_image1_roi = image1_roi;
	gen_image2_roi = image2_roi;
	%prepare a matrix for calculating the average mask of all images
	if ensemble_i1==1
		average_mask=zeros(size(image1_roi));
	end
	%get mask from mask list
	ximask={};
	yimask={};
	if size(maskiererx,2)>=ensemble_i1
		for j=1:size(maskiererx,1)
			if isempty(maskiererx{j,ensemble_i1})==0
				ximask{j,1}=maskiererx{j,ensemble_i1}; %#ok<*AGROW>
				yimask{j,1}=maskierery{j,ensemble_i1};
			else
				break
			end
		end
		if size(ximask,1)>0
			mask_inpt=[ximask yimask];
		else
			mask_inpt=[];
		end
	else
		mask_inpt=[];
	end
	if numel(mask_inpt)>0
		cellmask=mask_inpt;
		mask=zeros(size(image1_roi));
		for i=1:size(cellmask,1)
			masklayerx=cellmask{i,1};
			masklayery=cellmask{i,2};
			mask = mask + poly2mask(masklayerx-xroi,masklayery-yroi,size(image1_roi,1),size(image1_roi,2)); %kleineres eingangsbild und maske geshiftet
		end
	else
		mask=zeros(size(image1_roi));
	end
	mask(mask>1)=1;
	gen_mask = mask;
	try
		average_mask=average_mask + mask; %will fail if images are not same dimensions.
	catch
		cancel = 1;
		hgui=getappdata(0,'hgui');
		setappdata(hgui, 'cancel', cancel);
		text(gca(getappdata(0,'hgui')),10,10,'Error: Image dimensions inconsistent!','color',[1 0 0],'fontsize',20)
		drawnow;
		break
	end
	miniy=1+(ceil(interrogationarea/2));
	minix=1+(ceil(interrogationarea/2));
	maxiy=step*(floor(size(image1_roi,1)/step))-(interrogationarea-1)+(ceil(interrogationarea/2)); %statt size deltax von ROI nehmen
	maxix=step*(floor(size(image1_roi,2)/step))-(interrogationarea-1)+(ceil(interrogationarea/2));

	numelementsy=floor((maxiy-miniy)/step+1);
	numelementsx=floor((maxix-minix)/step+1);

	LAy=miniy;
	LAx=minix;
	LUy=size(image1_roi,1)-maxiy;
	LUx=size(image1_roi,2)-maxix;
	shift4centery=round((LUy-LAy)/2);
	shift4centerx=round((LUx-LAx)/2);
	if shift4centery<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data....
		shift4centery=0;
	end
	if shift4centerx<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data....
		shift4centerx=0;
	end
	miniy=miniy+shift4centery;
	minix=minix+shift4centerx;
	maxix=maxix+shift4centerx;
	maxiy=maxiy+shift4centery;

	image1_roi=padarray(image1_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi)));
	image2_roi=padarray(image2_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi)));
	mask=padarray(mask,[ceil(interrogationarea/2) ceil(interrogationarea/2)],0);

	if (rem(interrogationarea,2) == 0) %for the subpixel displacement measurement
		SubPixOffset=1;
	else
		SubPixOffset=0.5;
	end
	xtable=zeros(numelementsy,numelementsx);
	ytable=xtable; %#ok<*NASGU>
	utable=xtable;
	vtable=xtable;
	typevector=ones(numelementsy,numelementsx);

	%% MAINLOOP
	try %check if used from GUI
		handles=guihandles(getappdata(0,'hgui'));
		GUI_avail=1;
		hgui=getappdata(0,'hgui');
		cancel=getappdata(hgui, 'cancel');
		if cancel == 1
			break
			%disp('user cancelled');
		end

	catch %#ok<CTCH>
		GUI_avail=0;
		disp('no GUI')
	end
	% divide images by small pictures
	% new index for image1_roi and image2_roi
	s0 = (repmat((miniy:step:maxiy)'-1, 1,numelementsx) + repmat(((minix:step:maxix)-1)*size(image1_roi, 1), numelementsy,1))';
	s0 = permute(s0(:), [2 3 1]);
	s1 = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
	ss1 = repmat(s1, [1, 1, size(s0,3)])+repmat(s0, [interrogationarea, interrogationarea, 1]);

	image1_cut = image1_roi(ss1);
	image2_cut = image2_roi(ss1);

	if do_pad==1 && passes == 1 %only on first pass
		%subtract mean to avoid high frequencies at border of correlation:
		image1_cut=image1_cut-mean(image1_cut,[1 2]);
		image2_cut=image2_cut-mean(image2_cut,[1 2]);
		% padding (faster than padarray) to get the linear correlation:
		image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
		image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];
	end
	%do fft2:

	result_conv = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);
	if do_pad==1 && passes == 1
		%cropping of correlation matrix:
		result_conv =result_conv((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
	end

	%% repeated  Correlation in the first pass (might make sense to repeat more often to make it even more robust...)
	if repeat == 1 && passes == 1
		ms=round(step/4); %multishift parameter so groß wie viertel int window
		%Shift left bot
		s0B = (repmat((miniy+ms:step:maxiy+ms)'-1, 1,numelementsx) + repmat(((minix-ms:step:maxix-ms)-1)*size(image1_roi, 1), numelementsy,1))';
		s0B = permute(s0B(:), [2 3 1]);
		s1B = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
		ss1B = repmat(s1B, [1, 1, size(s0B,3)])+repmat(s0B, [interrogationarea, interrogationarea, 1]);
		image1_cutB = image1_roi(ss1B);
		image2_cutB = image2_roi(ss1B);
		if do_pad==1 && passes == 1
			%subtract mean to avoid high frequencies at border of correlation:
			image1_cutB=image1_cutB-mean(image1_cutB,[1 2]);
			image2_cutB=image2_cutB-mean(image2_cutB,[1 2]);
			% padding (faster than padarray) to get the linear correlation:
			image1_cutB=[image1_cutB zeros(interrogationarea,interrogationarea-1,size(image1_cutB,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cutB,3))];
			image2_cutB=[image2_cutB zeros(interrogationarea,interrogationarea-1,size(image2_cutB,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cutB,3))];
		end
		result_convB = fftshift(fftshift(real(ifft2(conj(fft2(image1_cutB)).*fft2(image2_cutB))), 1), 2);
		if do_pad==1 && passes == 1
			%cropping of correlation matrix:
			result_convB =result_convB((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
		end

		%Shift right bot
		s0C = (repmat((miniy+ms:step:maxiy+ms)'-1, 1,numelementsx) + repmat(((minix+ms:step:maxix+ms)-1)*size(image1_roi, 1), numelementsy,1))';
		s0C = permute(s0C(:), [2 3 1]);
		s1C = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
		ss1C = repmat(s1C, [1, 1, size(s0C,3)])+repmat(s0C, [interrogationarea, interrogationarea, 1]);
		image1_cutC = image1_roi(ss1C);
		image2_cutC = image2_roi(ss1C);
		if do_pad==1 && passes == 1
			%subtract mean to avoid high frequencies at border of correlation:
			image1_cutC=image1_cutC-mean(image1_cutC,[1 2]);
			image2_cutC=image2_cutC-mean(image2_cutC,[1 2]);
			% padding (faster than padarray) to get the linear correlation:
			image1_cutC=[image1_cutC zeros(interrogationarea,interrogationarea-1,size(image1_cutC,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cutC,3))];
			image2_cutC=[image2_cutC zeros(interrogationarea,interrogationarea-1,size(image2_cutC,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cutC,3))];
		end
		result_convC = fftshift(fftshift(real(ifft2(conj(fft2(image1_cutC)).*fft2(image2_cutC))), 1), 2);
		if do_pad==1 && passes == 1
			%cropping of correlation matrix:
			result_convC =result_convC((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
		end

		%Shift left top
		s0D = (repmat((miniy-ms:step:maxiy-ms)'-1, 1,numelementsx) + repmat(((minix-ms:step:maxix-ms)-1)*size(image1_roi, 1), numelementsy,1))';
		s0D = permute(s0D(:), [2 3 1]);
		s1D = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
		ss1D = repmat(s1D, [1, 1, size(s0D,3)])+repmat(s0D, [interrogationarea, interrogationarea, 1]);
		image1_cutD = image1_roi(ss1D);
		image2_cutD = image2_roi(ss1D);

		if do_pad==1 && passes == 1
			%subtract mean to avoid high frequencies at border of correlation:
			image1_cutD=image1_cutD-mean(image1_cutD,[1 2]);
			image2_cutD=image2_cutD-mean(image2_cutD,[1 2]);
			% padding (faster than padarray) to get the linear correlation:
			image1_cutD=[image1_cutD zeros(interrogationarea,interrogationarea-1,size(image1_cutD,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cutD,3))];
			image2_cutD=[image2_cutD zeros(interrogationarea,interrogationarea-1,size(image2_cutD,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cutD,3))];
		end
		result_convD = fftshift(fftshift(real(ifft2(conj(fft2(image1_cutD)).*fft2(image2_cutD))), 1), 2);
		if do_pad==1 && passes == 1
			%cropping of correlation matrix:
			result_convD =result_convD((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
		end

		%Shift right top
		s0E = (repmat((miniy-ms:step:maxiy-ms)'-1, 1,numelementsx) + repmat(((minix+ms:step:maxix+ms)-1)*size(image1_roi, 1), numelementsy,1))';
		s0E = permute(s0E(:), [2 3 1]);
		s1E = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
		ss1E = repmat(s1E, [1, 1, size(s0E,3)])+repmat(s0E, [interrogationarea, interrogationarea, 1]);
		image1_cutE = image1_roi(ss1E);
		image2_cutE = image2_roi(ss1E);
		if do_pad==1 && passes == 1
			%subtract mean to avoid high frequencies at border of correlation:
			image1_cutE=image1_cutE-mean(image1_cutE,[1 2]);
			image2_cutE=image2_cutE-mean(image2_cutE,[1 2]);
			% padding (faster than padarray) to get the linear correlation:
			image1_cutE=[image1_cutE zeros(interrogationarea,interrogationarea-1,size(image1_cutE,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cutE,3))];
			image2_cutE=[image2_cutE zeros(interrogationarea,interrogationarea-1,size(image2_cutE,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cutE,3))];
		end
		result_convE = fftshift(fftshift(real(ifft2(conj(fft2(image1_cutE)).*fft2(image2_cutE))), 1), 2);
		if do_pad==1 && passes == 1
			%cropping of correlation matrix:
			result_convE =result_convE((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
		end
		result_conv=result_conv.*result_convB.*result_convC.*result_convD.*result_convE;
	end

	if mask_auto == 1
		%das zentrum der Matrize (3x3) mit dem mittelwert ersetzen = Keine Autokorrelation
		%MARKER
		h = fspecial('gaussian', 3, 1.5);
		h=h/h(2,2);
		h=1-h;
		%h=repmat(h,1,1,size(result_conv,3));
		h=repmat(h,[1,1,size(result_conv,3)]);
		h=h.*result_conv((interrogationarea/2)+SubPixOffset-1:(interrogationarea/2)+SubPixOffset+1,(interrogationarea/2)+SubPixOffset-1:(interrogationarea/2)+SubPixOffset+1,:);
		result_conv((interrogationarea/2)+SubPixOffset-1:(interrogationarea/2)+SubPixOffset+1,(interrogationarea/2)+SubPixOffset-1:(interrogationarea/2)+SubPixOffset+1,:)=h;
	end
	%apply mask
	ii = find(mask(ss1(round(interrogationarea/2+1), round(interrogationarea/2+1), :)));
	result_conv(:,:, ii) = 0;
	%average the correlation matrices
	try
		result_conv_ensemble=result_conv_ensemble+result_conv;
	catch % older matlab releases
		result_conv_ensemble = zeros(size(result_conv));
		result_conv_ensemble=result_conv_ensemble+result_conv;
	end

	if GUI_avail==1
		progri=ensemble_i1/(amount_input_imgs)*100;
		from_total=from_total+1;
		set(handles.progress, 'string' , ['Pass ' int2str(1) ' progress: ' int2str(progri) '%' ])
		set(handles.overall, 'string' , ['Total progress: ' int2str(from_total / total_analyses_amount * 100) '%'])
		zeit=toc;
		done=from_total;
		tocome=total_analyses_amount-done;
		zeit=zeit/done*tocome;
		hrs=zeit/60^2;
		mins=(hrs-floor(hrs))*60;
		secs=(mins-floor(mins))*60;
		hrs=floor(hrs);
		mins=floor(mins);
		secs=floor(secs);
		set(handles.totaltime,'string', ['Time left: ' sprintf('%2.2d', hrs) 'h ' sprintf('%2.2d', mins) 'm ' sprintf('%2.2d', secs) 's']);

		%xxx update display every 10 frames...?
		%aber wie, dann müsste man peakfinder machen
		if skippy ==0
			[xtable,ytable,utable, vtable] = peakfinding (result_conv_ensemble, mask, interrogationarea,minix,step,maxix,miniy,maxiy,SubPixOffset,ss1,subpixfinder);
			if verLessThan('matlab','8.4')
				delete (findobj(getappdata(0,'hgui'),'type', 'hggroup'))
			else
				delete (findobj(getappdata(0,'hgui'),'type', 'quiver'))
			end
			hold on;
			vecscale=str2double(get(handles.vectorscale,'string'));
			%Problem: wenn colorbar an, zï¿½hlt das auch als aexes...
			colorbar('off')

			%u_table original gibts nicjt, braichts auch nicht...
			quiver ((findobj(getappdata(0,'hgui'),'type', 'axes')),xtable(isnan(utable)==0)+xroi-interrogationarea/2,ytable(isnan(utable)==0)+yroi-interrogationarea/2,utable(isnan(utable)==0)*vecscale,vtable(isnan(utable)==0)*vecscale,'Color', [1-(from_total / total_analyses_amount) (from_total / total_analyses_amount) 0.15],'autoscale','off')
			%quiver ((findobj(getappdata(0,'hgui'),'type', 'axes')),xtable(isnan(utable)==1)+xroi-interrogationarea/2,ytable(isnan(utable)==1)+yroi-interrogationarea/2,utable(isnan(utable)==1)*vecscale,vtable(isnan(utable)==1)*vecscale,'Color',[0.7 0.15 0.15], 'autoscale','off')
			hold off
			drawnow;
		end
		if skippy <10
			skippy=skippy+1;
		else
			skippy=0;
		end
		try
			drawnow limitrate
		catch
			drawnow
		end
	else
		fprintf('.');
	end
	if passes==1 % only 1 pass selected, so correlation coefficient will be calculated in this (first & final) pass.
		if ensemble_i1==1 %first image pair
			correlation_map=zeros(size(typevector));
			corr_map_cnt=0;
		end
		for cor_i=1:size(image1_cut,3)
			correlation_map(cor_i)=correlation_map(cor_i) + corr2(image1_cut(:,:,cor_i),image2_cut(:,:,cor_i));
		end
		corr_map_cnt=corr_map_cnt+1;
	end
end
%correlation_map=[];
if cancel == 0
	%% Correlation matrix of pass 1 is done.
	[xtable,ytable,utable, vtable] = peakfinding (result_conv_ensemble, mask, interrogationarea,minix,step,maxix,miniy,maxiy,SubPixOffset,ss1,subpixfinder);

	for multipass=1:passes-1
		% unfortunately, preprocessing has to be done again for every pass, otherwise i would have to save the modified data somehow.
		if multipass==1
			interrogationarea=round(int2/2)*2;
		end
		if multipass==2
			interrogationarea=round(int3/2)*2;
		end
		if multipass==3
			interrogationarea=round(int4/2)*2;
		end
		result_conv_ensemble = zeros(interrogationarea,interrogationarea); % prepare empty result_conv
		skippy=0;
		for ensemble_i1=1:2:amount_input_imgs
			if skippy <10
				skippy=skippy+1;
			else
				skippy=0;
			end
			if isempty(video_frame_selection) %list with image files was passed
				if strcmp(ext,'.b16')
					image1=f_readB16(filepath{ensemble_i1});
					image2=f_readB16(filepath{ensemble_i1+1});
				else
					image1=imread(filepath{ensemble_i1});
					image2=imread(filepath{ensemble_i1+1});
				end
			else % video file was passed
				image1 = read(filepath,video_frame_selection(ensemble_i1));
				image2 = read(filepath,video_frame_selection(ensemble_i1+1));
			end
			if size(image1,3)>1
				image1=uint8(mean(image1,3));
				image2=uint8(mean(image2,3));
			end
			%subtract bg if present
			if ~isempty(bg_img_A)
				image1=image1-bg_img_A;
			end
			if ~isempty(bg_img_B)
				image2=image2-bg_img_B;
			end
			%if autolimit == 1 %if autolimit is desired: do autolimit for each image seperately
				if size(image1,3)>1
					stretcher = stretchlim(rgb2gray(image1));
				else
					stretcher = stretchlim(image1);
				end
				minintens1 = stretcher(1);
				maxintens1 = stretcher(2);
				if size(image2,3)>1
					stretcher = stretchlim(rgb2gray(image2));
				else
					stretcher = stretchlim(image2);
				end
				minintens2 = stretcher(1);
				maxintens2 = stretcher(2);
			%end
			image1 = PIVlab_preproc (image1,roi_inpt,clahe, clahesize,highp,highpsize,intenscap,wienerwurst,wienerwurstsize,minintens1,maxintens1);
			image2 = PIVlab_preproc (image2,roi_inpt,clahe, clahesize,highp,highpsize,intenscap,wienerwurst,wienerwurstsize,minintens2,maxintens2);
			if numel(roi_inpt)>0
				xroi=roi_inpt(1);
				yroi=roi_inpt(2);
				widthroi=roi_inpt(3);
				heightroi=roi_inpt(4);
				image1_roi=double(image1(yroi:yroi+heightroi,xroi:xroi+widthroi));
				image2_roi=double(image2(yroi:yroi+heightroi,xroi:xroi+widthroi));
			else
				xroi=0;
				yroi=0;
				image1_roi=double(image1);
				image2_roi=double(image2);
			end
			gen_image1_roi = image1_roi;
			gen_image2_roi = image2_roi;
			if GUI_avail==1
				progri=ensemble_i1/(amount_input_imgs)*100;
				from_total=from_total+1;
				set(handles.progress, 'string' , ['Pass ' int2str(multipass+1) ' progress: ' int2str(progri) '%' ])
				set(handles.overall, 'string' , ['Total progress: ' int2str(from_total / total_analyses_amount * 100) '%'])

				zeit=toc;
				done=from_total;
				tocome=total_analyses_amount-done;
				zeit=zeit/done*tocome;
				hrs=zeit/60^2;
				mins=(hrs-floor(hrs))*60;
				secs=(mins-floor(mins))*60;
				hrs=floor(hrs);
				mins=floor(mins);
				secs=floor(secs);
				set(handles.totaltime,'string', ['Time left: ' sprintf('%2.2d', hrs) 'h ' sprintf('%2.2d', mins) 'm ' sprintf('%2.2d', secs) 's']);
				try
					drawnow limitrate

				catch
					drawnow
				end
			else
				fprintf('.');
			end
			%multipass validation, smoothing
			utable_orig=utable;
			vtable_orig=vtable;
			[utable,vtable] = PIVlab_postproc (utable,vtable,[],[], [], 1,4, 1,1.5);
			if GUI_avail==1
				cancel=getappdata(hgui, 'cancel');
				if cancel == 1
					break
					%disp('user cancelled');
				end
				if skippy ==0
					if verLessThan('matlab','8.4')
						delete (findobj(getappdata(0,'hgui'),'type', 'hggroup'))
					else
						delete (findobj(getappdata(0,'hgui'),'type', 'quiver'))
					end
					hold on;
					vecscale=str2double(get(handles.vectorscale,'string'));
					%Problem: wenn colorbar an, zï¿½hlt das auch als aexes...
					colorbar('off')
					quiver ((findobj(getappdata(0,'hgui'),'type', 'axes')),xtable(isnan(utable)==0)+xroi-interrogationarea/2,ytable(isnan(utable)==0)+yroi-interrogationarea/2,utable(isnan(utable)==0)*vecscale,vtable(isnan(utable)==0)*vecscale,'Color', [1-(from_total / total_analyses_amount) (from_total / total_analyses_amount) 0.15],'autoscale','off')

					%                    quiver ((findobj(getappdata(0,'hgui'),'type', 'axes')),xtable(isnan(utable)==0)+xroi-interrogationarea/2,ytable(isnan(utable)==0)+yroi-interrogationarea/2,utable_orig(isnan(utable)==0)*vecscale,vtable_orig(isnan(utable)==0)*vecscale,'Color', [0.15 0.7 0.15],'autoscale','off')
					%quiver ((findobj(getappdata(0,'hgui'),'type', 'axes')),xtable(isnan(utable)==1)+xroi-interrogationarea/2,ytable(isnan(utable)==1)+yroi-interrogationarea/2,utable_orig(isnan(utable)==1)*vecscale,vtable_orig(isnan(utable)==1)*vecscale,'Color',[0.7 0.15 0.15], 'autoscale','off')
					drawnow
					hold off
				end
			end
			%replace nans
			utable=inpaint_nans(utable,4);
			vtable=inpaint_nans(vtable,4);
			%smooth predictor
			try
				if multipass<passes-1
					utable = smoothn(utable,0.6); %stronger smoothing for first passes
					vtable = smoothn(vtable,0.6);
				else
					utable = smoothn(utable); %weaker smoothing for last pass
					vtable = smoothn(vtable);
				end
			catch
				%old matlab versions: gaussian kernel
				h=fspecial('gaussian',5,1);
				utable=imfilter(utable,h,'replicate');
				vtable=imfilter(vtable,h,'replicate');
			end
			if multipass==1
				interrogationarea=round(int2/2)*2;
			end
			if multipass==2
				interrogationarea=round(int3/2)*2;
			end
			if multipass==3
				interrogationarea=round(int4/2)*2;
			end
			step=interrogationarea/2;

			%bildkoordinaten neu errechnen:
			image1_roi = gen_image1_roi;
			image2_roi = gen_image2_roi;
			%get mask from mask list
			ximask={};
			yimask={};
			if size(maskiererx,2)>=ensemble_i1
				for j=1:size(maskiererx,1)
					if isempty(maskiererx{j,ensemble_i1})==0
						ximask{j,1}=maskiererx{j,ensemble_i1}; %#ok<*AGROW>
						yimask{j,1}=maskierery{j,ensemble_i1};
					else
						break
					end
				end
				if size(ximask,1)>0
					mask_inpt=[ximask yimask];
				else
					mask_inpt=[];
				end
			else
				mask_inpt=[];
			end
			if numel(mask_inpt)>0
				cellmask=mask_inpt;
				mask=zeros(size(image1_roi));
				for i=1:size(cellmask,1)
					masklayerx=cellmask{i,1};
					masklayery=cellmask{i,2};
					mask = mask + poly2mask(masklayerx-xroi,masklayery-yroi,size(image1_roi,1),size(image1_roi,2)); %kleineres eingangsbild und maske geshiftet
				end
			else
				mask=zeros(size(image1_roi));
			end
			mask(mask>1)=1;
			gen_mask = mask;
			miniy=1+(ceil(interrogationarea/2));
			minix=1+(ceil(interrogationarea/2));
			maxiy=step*(floor(size(image1_roi,1)/step))-(interrogationarea-1)+(ceil(interrogationarea/2)); %statt size deltax von ROI nehmen
			maxix=step*(floor(size(image1_roi,2)/step))-(interrogationarea-1)+(ceil(interrogationarea/2));

			numelementsy=floor((maxiy-miniy)/step+1);
			numelementsx=floor((maxix-minix)/step+1);

			LAy=miniy;
			LAx=minix;
			LUy=size(image1_roi,1)-maxiy;
			LUx=size(image1_roi,2)-maxix;
			shift4centery=round((LUy-LAy)/2);
			shift4centerx=round((LUx-LAx)/2);
			if shift4centery<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data....
				shift4centery=0;
			end
			if shift4centerx<0 %shift4center will be negative if in the unshifted case the left border is bigger than the right border. the vectormatrix is hence not centered on the image. the matrix cannot be shifted more towards the left border because then image2_crop would have a negative index. The only way to center the matrix would be to remove a column of vectors on the right side. but then we weould have less data....
				shift4centerx=0;
			end
			miniy=miniy+shift4centery;
			minix=minix+shift4centerx;
			maxix=maxix+shift4centerx;
			maxiy=maxiy+shift4centery;

			image1_roi=padarray(image1_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi)));
			image2_roi=padarray(image2_roi,[ceil(interrogationarea/2) ceil(interrogationarea/2)], min(min(image1_roi)));
			mask=padarray(mask,[ceil(interrogationarea/2) ceil(interrogationarea/2)],0);
			if (rem(interrogationarea,2) == 0) %for the subpixel displacement measurement
				SubPixOffset=1;
			else
				SubPixOffset=0.5;
			end

			xtable_old=xtable;
			ytable_old=ytable;
			typevector=ones(numelementsy,numelementsx);
			xtable = repmat((minix:step:maxix), numelementsy, 1) + interrogationarea/2;
			ytable = repmat((miniy:step:maxiy)', 1, numelementsx) + interrogationarea/2;

			%xtable alt und neu geben koordinaten wo die vektoren herkommen.
			%d.h. u und v auf die gewï¿½nschte grï¿½ï¿½e bringen+interpolieren

			utable=interp2(xtable_old,ytable_old,utable,xtable,ytable,'*spline');
			vtable=interp2(xtable_old,ytable_old,vtable,xtable,ytable,'*spline');

			utable_1= padarray(utable, [1,1], 'replicate');
			vtable_1= padarray(vtable, [1,1], 'replicate');

			%add 1 line around image for border regions... linear extrap

			firstlinex=xtable(1,:);
			firstlinex_intp=interp1(1:1:size(firstlinex,2),firstlinex,0:1:size(firstlinex,2)+1,'linear','extrap');
			xtable_1=repmat(firstlinex_intp,size(xtable,1)+2,1);

			firstliney=ytable(:,1);
			firstliney_intp=interp1(1:1:size(firstliney,1),firstliney,0:1:size(firstliney,1)+1,'linear','extrap')';
			ytable_1=repmat(firstliney_intp,1,size(ytable,2)+2);

			X=xtable_1; %original locations of vectors in whole image
			Y=ytable_1;
			U=utable_1; %interesting portion of u
			V=vtable_1; % "" of v

			X1=X(1,1):1:X(1,end)-1;
			Y1=(Y(1,1):1:Y(end,1)-1)';
			X1=repmat(X1,size(Y1, 1),1);
			Y1=repmat(Y1,1,size(X1, 2));

			U1 = interp2(X,Y,U,X1,Y1,'*linear');
			V1 = interp2(X,Y,V,X1,Y1,'*linear');

			image2_crop_i1 = interp2(1:size(image2_roi,2),(1:size(image2_roi,1))',double(image2_roi),X1+U1,Y1+V1,imdeform); %linear is 3x faster and looks ok...

			xb = find(X1(1,:) == xtable_1(1,1));
			yb = find(Y1(:,1) == ytable_1(1,1));

			% divide images by small pictures
			% new index for image1_roi
			s0 = (repmat((miniy:step:maxiy)'-1, 1,numelementsx) + repmat(((minix:step:maxix)-1)*size(image1_roi, 1), numelementsy,1))';
			s0 = permute(s0(:), [2 3 1]);
			s1 = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
			ss1 = repmat(s1, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
			% new index for image2_crop_i1
			s0 = (repmat(yb-step+step*(1:numelementsy)'-1, 1,numelementsx) + repmat((xb-step+step*(1:numelementsx)-1)*size(image2_crop_i1, 1), numelementsy,1))';
			s0 = permute(s0(:), [2 3 1]) - s0(1);
			s2 = repmat((1:2*step)',1,2*step) + repmat(((1:2*step)-1)*size(image2_crop_i1, 1),2*step,1);
			ss2 = repmat(s2, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
			image1_cut = image1_roi(ss1);
			image2_cut = image2_crop_i1(ss2);
			if do_pad==1 && multipass==passes-1
				%subtract mean to avoid high frequencies at border of correlation:
				try
					image1_cut=image1_cut-mean(image1_cut,[1 2]);
					image2_cut=image2_cut-mean(image2_cut,[1 2]);
				catch
					for oldmatlab=1:size(image1_cut,3);
						image1_cut(:,:,oldmatlab)=image1_cut(:,:,oldmatlab)-mean(mean(image1_cut(:,:,oldmatlab)));
						image2_cut(:,:,oldmatlab)=image2_cut(:,:,oldmatlab)-mean(mean(image2_cut(:,:,oldmatlab)));
					end
				end
				% padding (faster than padarray) to get the linear correlation:
				image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
				image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];
			end
			%do fft2:
			result_conv = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);
			if do_pad==1 && multipass==passes-1
				%cropping of correlation matrix:
				result_conv =result_conv((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
			end

			%% repeated correlation
			if repeat == 1 && multipass==passes-1
				ms=round(step/4); %multishift parameter so groß wie viertel int window

				%Shift left bot
				image2_crop_i1 = interp2(1:size(image2_roi,2),(1:size(image2_roi,1))',double(image2_roi),X1+U1-ms,Y1+V1+ms,imdeform); %linear is 3x faster and looks ok...
				xb = find(X1(1,:) == xtable_1(1,1));
				yb = find(Y1(:,1) == ytable_1(1,1));
				s0 = (repmat((miniy+ms:step:maxiy+ms)'-1, 1,numelementsx) + repmat(((minix-ms:step:maxix-ms)-1)*size(image1_roi, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]);
				s1 = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
				ss1 = repmat(s1, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				s0 = (repmat(yb-step+step*(1:numelementsy)'-1, 1,numelementsx) + repmat((xb-step+step*(1:numelementsx)-1)*size(image2_crop_i1, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]) - s0(1);
				s2 = repmat((1:2*step)',1,2*step) + repmat(((1:2*step)-1)*size(image2_crop_i1, 1),2*step,1);
				ss2 = repmat(s2, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				image1_cut = image1_roi(ss1);
				image2_cut = image2_crop_i1(ss2);
				if do_pad==1 && multipass==passes-1
					%subtract mean to avoid high frequencies at border of correlation:
					try
						image1_cut=image1_cut-mean(image1_cut,[1 2]);
						image2_cut=image2_cut-mean(image2_cut,[1 2]);
					catch
						for oldmatlab=1:size(image1_cut,3);
							image1_cut(:,:,oldmatlab)=image1_cut(:,:,oldmatlab)-mean(mean(image1_cut(:,:,oldmatlab)));
							image2_cut(:,:,oldmatlab)=image2_cut(:,:,oldmatlab)-mean(mean(image2_cut(:,:,oldmatlab)));
						end
					end
					% padding (faster than padarray) to get the linear correlation:
					image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
					image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];
				end
				result_convB = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);
				if do_pad==1 && multipass==passes-1
					%cropping of correlation matrix:
					result_convB =result_convB((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
				end
				%Shift right bot
				image2_crop_i1 = interp2(1:size(image2_roi,2),(1:size(image2_roi,1))',double(image2_roi),X1+U1+ms,Y1+V1+ms,imdeform); %linear is 3x faster and looks ok...
				xb = find(X1(1,:) == xtable_1(1,1));
				yb = find(Y1(:,1) == ytable_1(1,1));
				s0 = (repmat((miniy+ms:step:maxiy+ms)'-1, 1,numelementsx) + repmat(((minix+ms:step:maxix+ms)-1)*size(image1_roi, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]);
				s1 = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
				ss1 = repmat(s1, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				s0 = (repmat(yb-step+step*(1:numelementsy)'-1, 1,numelementsx) + repmat((xb-step+step*(1:numelementsx)-1)*size(image2_crop_i1, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]) - s0(1);
				s2 = repmat((1:2*step)',1,2*step) + repmat(((1:2*step)-1)*size(image2_crop_i1, 1),2*step,1);
				ss2 = repmat(s2, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				image1_cut = image1_roi(ss1);
				image2_cut = image2_crop_i1(ss2);
				if do_pad==1 && multipass==passes-1
					%subtract mean to avoid high frequencies at border of correlation:
					try
						image1_cut=image1_cut-mean(image1_cut,[1 2]);
						image2_cut=image2_cut-mean(image2_cut,[1 2]);
					catch
						for oldmatlab=1:size(image1_cut,3);
							image1_cut(:,:,oldmatlab)=image1_cut(:,:,oldmatlab)-mean(mean(image1_cut(:,:,oldmatlab)));
							image2_cut(:,:,oldmatlab)=image2_cut(:,:,oldmatlab)-mean(mean(image2_cut(:,:,oldmatlab)));
						end
					end
					% padding (faster than padarray) to get the linear correlation:
					image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
					image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];
				end
				result_convC = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);
				if do_pad==1 && multipass==passes-1
					%cropping of correlation matrix:
					result_convC =result_convC((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
				end
				%Shift left top
				image2_crop_i1 = interp2(1:size(image2_roi,2),(1:size(image2_roi,1))',double(image2_roi),X1+U1-ms,Y1+V1-ms,imdeform); %linear is 3x faster and looks ok...
				xb = find(X1(1,:) == xtable_1(1,1));
				yb = find(Y1(:,1) == ytable_1(1,1));
				s0 = (repmat((miniy-ms:step:maxiy-ms)'-1, 1,numelementsx) + repmat(((minix-ms:step:maxix-ms)-1)*size(image1_roi, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]);
				s1 = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
				ss1 = repmat(s1, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				s0 = (repmat(yb-step+step*(1:numelementsy)'-1, 1,numelementsx) + repmat((xb-step+step*(1:numelementsx)-1)*size(image2_crop_i1, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]) - s0(1);
				s2 = repmat((1:2*step)',1,2*step) + repmat(((1:2*step)-1)*size(image2_crop_i1, 1),2*step,1);
				ss2 = repmat(s2, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				image1_cut = image1_roi(ss1);
				image2_cut = image2_crop_i1(ss2);
				if do_pad==1 && multipass==passes-1
					%subtract mean to avoid high frequencies at border of correlation:
					try
						image1_cut=image1_cut-mean(image1_cut,[1 2]);
						image2_cut=image2_cut-mean(image2_cut,[1 2]);
					catch
						for oldmatlab=1:size(image1_cut,3)
							image1_cut(:,:,oldmatlab)=image1_cut(:,:,oldmatlab)-mean(mean(image1_cut(:,:,oldmatlab)));
							image2_cut(:,:,oldmatlab)=image2_cut(:,:,oldmatlab)-mean(mean(image2_cut(:,:,oldmatlab)));
						end
					end
					% padding (faster than padarray) to get the linear correlation:
					image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
					image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];
				end
				result_convD = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);
				if do_pad==1 && multipass==passes-1
					%cropping of correlation matrix:
					result_convD =result_convD((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
				end
				%Shift right top
				image2_crop_i1 = interp2(1:size(image2_roi,2),(1:size(image2_roi,1))',double(image2_roi),X1+U1+ms,Y1+V1-ms,imdeform); %linear is 3x faster and looks ok...
				xb = find(X1(1,:) == xtable_1(1,1));
				yb = find(Y1(:,1) == ytable_1(1,1));
				s0 = (repmat((miniy-ms:step:maxiy-ms)'-1, 1,numelementsx) + repmat(((minix+ms:step:maxix+ms)-1)*size(image1_roi, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]);
				s1 = repmat((1:interrogationarea)',1,interrogationarea) + repmat(((1:interrogationarea)-1)*size(image1_roi, 1),interrogationarea,1);
				ss1 = repmat(s1, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				s0 = (repmat(yb-step+step*(1:numelementsy)'-1, 1,numelementsx) + repmat((xb-step+step*(1:numelementsx)-1)*size(image2_crop_i1, 1), numelementsy,1))';
				s0 = permute(s0(:), [2 3 1]) - s0(1);
				s2 = repmat((1:2*step)',1,2*step) + repmat(((1:2*step)-1)*size(image2_crop_i1, 1),2*step,1);
				ss2 = repmat(s2, [1, 1, size(s0,3)]) + repmat(s0, [interrogationarea, interrogationarea, 1]);
				image1_cut = image1_roi(ss1);
				image2_cut = image2_crop_i1(ss2);
				if do_pad==1 && multipass==passes-1
					%subtract mean to avoid high frequencies at border of correlation:
					try
						image1_cut=image1_cut-mean(image1_cut,[1 2]);
						image2_cut=image2_cut-mean(image2_cut,[1 2]);
					catch
						for oldmatlab=1:size(image1_cut,3)
							image1_cut(:,:,oldmatlab)=image1_cut(:,:,oldmatlab)-mean(mean(image1_cut(:,:,oldmatlab)));
							image2_cut(:,:,oldmatlab)=image2_cut(:,:,oldmatlab)-mean(mean(image2_cut(:,:,oldmatlab)));
						end
					end
					% padding (faster than padarray) to get the linear correlation:
					image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
					image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];
				end
				result_convE = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);
				if do_pad==1 && multipass==passes-1
					%cropping of correlation matrix:
					result_convE =result_convE((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);
				end
				result_conv=result_conv.*result_convB.*result_convC.*result_convD.*result_convE;
			end

			if mask_auto == 1
				%limit peak search arena....
				emptymatrix=zeros(size(result_conv,1),size(result_conv,2),size(result_conv,3));
				%emptymatrix=emptymatrix+0.1;
				if interrogationarea > 8 % masking central peak will not work for extrmely small interrogation areas. And it also doesn't make sense.
					sizeones=4;
					%h = fspecial('gaussian', sizeones*2+1,1);
					h=fspecial('disk',4);
					h=h/max(max(h));
					%h=repmat(h,1,1,size(result_conv,3));
					h=repmat(h,[1,1,size(result_conv,3)]);
					emptymatrix((interrogationarea/2)+SubPixOffset-sizeones:(interrogationarea/2)+SubPixOffset+sizeones,(interrogationarea/2)+SubPixOffset-sizeones:(interrogationarea/2)+SubPixOffset+sizeones,:)=h;
					result_conv = result_conv .* emptymatrix;
				else
					disp('All interrogation areas must be larger than 8 pixels for disabling auto correlation successfully.')
				end
			end

			%apply mask ---
			ii = find(mask(ss1(round(interrogationarea/2+1), round(interrogationarea/2+1), :)));
			result_conv(:,:, ii) = 0;
			%add alle result_conv
			try
				result_conv_ensemble=result_conv_ensemble+result_conv;
			catch % older matlab releases
				result_conv_ensemble = zeros(size(result_conv));
				result_conv_ensemble=result_conv_ensemble+result_conv;
			end

			if multipass==passes-1 %correlation strength only in last pass

				if ensemble_i1==1 %first image pair
					correlation_map=zeros(size(typevector));
					corr_map_cnt=0;
				end
				%Correlation strength
				for cor_i=1:size(image1_cut,3)
					correlation_map(cor_i)=correlation_map(cor_i)+corr2(image1_cut(:,:,cor_i),image2_cut(:,:,cor_i));
				end
				corr_map_cnt=corr_map_cnt+1;
			end
		end
		[xtable,ytable,utable2, vtable2] = peakfinding (result_conv_ensemble, [], interrogationarea,minix,step,maxix,miniy,maxiy,SubPixOffset,ss1,subpixfinder);
		utable = utable+utable2;
		vtable = vtable+vtable2;
	end
	if cancel == 0
		%mask only if all frames are masked
		%apply mask
		nrx=0;
		nrxreal=0;
		nry=0;
		average_mask=padarray(average_mask,[ceil(interrogationarea/2) ceil(interrogationarea/2)],0);
		for jmask = miniy:step:maxiy %vertical loop
			nry=nry+1;
			for imask = minix:step:maxix % horizontal loop
				nrx=nrx+1;%used to determine the pos of the vector in resulting matrix
				if nrxreal < numelementsx
					nrxreal=nrxreal+1;
				else
					nrxreal=1;
				end
				%fehlerzeile:
				if average_mask(round(jmask+interrogationarea/2),round(imask+interrogationarea/2)) >= amount_input_imgs/2
					typevector(nry,nrxreal)=0;
				end
			end
		end
		xtable=xtable-ceil(interrogationarea/2);
		ytable=ytable-ceil(interrogationarea/2);

		xtable=xtable+xroi;
		ytable=ytable+yroi;
	end
	%% Write correlation matrices to the workspace
	%{
try
    counter=evalin('base','counter');
    counter=counter+1;
    assignin('base','counter',counter);
    all_matrices=evalin('base','all_matrices');
    all_matrices{end+1}=result_conv_ensemble;
    assignin('base','all_matrices',all_matrices);
    disp('appended matrix')
catch
    assignin('base','counter',1);
    all_matrices{1}=result_conv_ensemble;
    assignin('base','all_matrices',all_matrices);
    disp('created new matrix')
end
	%}
	correlation_map = permute(reshape(correlation_map, [size(xtable')]), [2 1 3])/corr_map_cnt;
	%clear Correlation map in masked area
	correlation_map(typevector==0) = 0;
end

function [xtable,ytable,utable, vtable] = peakfinding (result_conv_ensemble, mask, interrogationarea,minix,step,maxix,miniy,maxiy,SubPixOffset,ss1,subpixfinder)
minres = permute(repmat(squeeze(min(min(result_conv_ensemble))), [1, size(result_conv_ensemble, 1), size(result_conv_ensemble, 2)]), [2 3 1]);
deltares = permute(repmat(squeeze(max(max(result_conv_ensemble))-min(min(result_conv_ensemble))),[ 1, size(result_conv_ensemble, 1), size(result_conv_ensemble, 2)]), [2 3 1]);
result_conv_ensemble = ((result_conv_ensemble-minres)./deltares)*255;

%apply mask ---
if isempty (mask)==0
	ii = find(mask(ss1(round(interrogationarea/2+1), round(interrogationarea/2+1), :)));
	result_conv_ensemble(:,:, ii) = 0;
end

[y, x, z] = ind2sub(size(result_conv_ensemble), find(result_conv_ensemble==255));

% we need only one peak from each couple pictures
[z1, zi] = sort(z);
dz1 = [z1(1); diff(z1)];
i0 = find(dz1~=0);
x1 = x(zi(i0));
y1 = y(zi(i0));
z1 = z(zi(i0));

xtable = repmat((minix:step:maxix)+interrogationarea/2, length(miniy:step:maxiy), 1);
ytable = repmat(((miniy:step:maxiy)+interrogationarea/2)', 1, length(minix:step:maxix));

if subpixfinder==1
	[vector] = SUBPIXGAUSS (result_conv_ensemble,interrogationarea, x1, y1, z1, SubPixOffset);
elseif subpixfinder==2
	[vector] = SUBPIX2DGAUSS (result_conv_ensemble,interrogationarea, x1, y1, z1, SubPixOffset);
end
vector = permute(reshape(vector, [size(xtable') 2]), [2 1 3]);

utable = vector(:,:,1);
vtable = vector(:,:,2);

function [vector] = SUBPIXGAUSS(result_conv, interrogationarea, x, y, z, SubPixOffset)
%was hat peak nr.1 für einen Durchmesser?
%figure;imagesc((1-im2bw(uint8(result_conv(:,:,155)),0.9)).*result_conv(:,:,101))
xi = find(~((x <= (size(result_conv,2)-1)) & (y <= (size(result_conv,1)-1)) & (x >= 2) & (y >= 2)));
x(xi) = [];
y(xi) = [];
z(xi) = [];
xmax = size(result_conv, 2);
vector = NaN(size(result_conv,3), 2);
if(numel(x)~=0)
	ip = sub2ind(size(result_conv), y, x, z);
	%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ï¿½ Israel Institute of Technology
	%http://urapiv.wordpress.com
	f0 = log(result_conv(ip));
	f1 = log(result_conv(ip-1));
	f2 = log(result_conv(ip+1));
	peaky = y + (f1-f2)./(2*f1-4*f0+2*f2);
	f0 = log(result_conv(ip));
	f1 = log(result_conv(ip-xmax));
	f2 = log(result_conv(ip+xmax));
	peakx = x + (f1-f2)./(2*f1-4*f0+2*f2);

	SubpixelX=peakx-(interrogationarea/2)-SubPixOffset;
	SubpixelY=peaky-(interrogationarea/2)-SubPixOffset;
	vector(z, :) = [SubpixelX, SubpixelY];
end

function [vector] = SUBPIX2DGAUSS(result_conv, interrogationarea, x, y, z, SubPixOffset)
xi = find(~((x <= (size(result_conv,2)-1)) & (y <= (size(result_conv,1)-1)) & (x >= 2) & (y >= 2)));
x(xi) = [];
y(xi) = [];
z(xi) = [];
xmax = size(result_conv, 2);
vector = NaN(size(result_conv,3), 2);
if(numel(x)~=0)
	c10 = zeros(3,3, length(z));
	c01 = c10;
	c11 = c10;
	c20 = c10;
	c02 = c10;
	ip = sub2ind(size(result_conv), y, x, z);

	for i = -1:1
		for j = -1:1
			%following 15 lines based on
			%H. Nobach ï¿½ M. Honkanen (2005)
			%Two-dimensional Gaussian regression for sub-pixel displacement
			%estimation in particle image velocimetry or particle position
			%estimation in particle tracking velocimetry
			%Experiments in Fluids (2005) 38: 511ï¿½515
			c10(j+2,i+2, :) = i*log(result_conv(ip+xmax*i+j));
			c01(j+2,i+2, :) = j*log(result_conv(ip+xmax*i+j));
			c11(j+2,i+2, :) = i*j*log(result_conv(ip+xmax*i+j));
			c20(j+2,i+2, :) = (3*i^2-2)*log(result_conv(ip+xmax*i+j));
			c02(j+2,i+2, :) = (3*j^2-2)*log(result_conv(ip+xmax*i+j));
			%c00(j+2,i+2)=(5-3*i^2-3*j^2)*log(result_conv_norm(maxY+j, maxX+i));
		end
	end
	c10 = (1/6)*sum(sum(c10));
	c01 = (1/6)*sum(sum(c01));
	c11 = (1/4)*sum(sum(c11));
	c20 = (1/6)*sum(sum(c20));
	c02 = (1/6)*sum(sum(c02));
	%c00=(1/9)*sum(sum(c00));

	deltax = squeeze((c11.*c01-2*c10.*c02)./(4*c20.*c02-c11.^2));
	deltay = squeeze((c11.*c10-2*c01.*c20)./(4*c20.*c02-c11.^2));
	peakx = x+deltax;
	peaky = y+deltay;

	SubpixelX = peakx-(interrogationarea/2)-SubPixOffset;
	SubpixelY = peaky-(interrogationarea/2)-SubPixOffset;

	vector(z, :) = [SubpixelX, SubpixelY];
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function varargout = expression( nm, varargin )

%EXPRESSION Declares a single CVX object for storing subexpressions.
%   EXPRESSION x
%   where x is a valid MATLAB variable nm, declares a scalar expression
%   holder for the current cvx model. Like a variable, an expression holder
%   can be used in constraints and objectives, according to the DCP ruleset.
%   However, unlike a variable, an expression holder is initialized to
%   zero, because the intent is for it to hold intermediate computations.
%
%   EXPRESSION x(n1,n2,...,nk)
%   declares a vector, matrix, or array expression holder with dimensions
%   n1, n2, ..., nk, each of which must be nonnegative integers. The value
%   of the expression holder is initialized to zero.
%
%   Examples:
%      variable x y
%      expression z
%      z = 2 * x - y;
%
%   See also EXPRESSIONS.

global cvx___
prob = evalin( 'caller', 'cvx_problem', '[]' );
if ~isa( prob, 'cvxprob' ),
    cvx_throw( 'No CVX model exists in this scope.' );
elseif isempty( cvx___.problems ) || cvx___.problems( end ).self ~= prob,
    cvx_throw( 'Internal CVX data corruption. Please CLEAR ALL and rebuild your model.' );
elseif nargin > 1,
    cvx_throw( 'Too many input arguments.\nTrying to declare multiple expression holders? Use the EXPRESSIONS keyword instead.', 1 ); %#ok
end

%
% Step 1: separate the name from the parenthetical, verify the name
%

toks = regexp( nm, '^\s*([a-zA-Z]\w*)\s*(\(.*\))?\s*$', 'tokens' );
if isempty( toks ),
    cvx_throw( 'Invalid variable specification: %s', nm );
end
toks = toks{1};
x.name = toks{1};
x.size = toks{2};
if x.name(end) == '_',
    cvx_throw( 'Invalid expression specification: %s\n   Names ending in underscores are reserved for internal use.', nm );
end

%
% Step 2: Parse the size. In effect, all we do here is surround what is
% replace the surrounding parentheses with square braces and evaluate. All
% that matters is the result is a valid size vector. In particular, it
% need to be a simple comma-delimited list.
%

if isempty( x.size ),
	x.size = [1,1];
else
    try
        x.size = evalin( 'caller', [ '[', x.size(2:end-1), '];' ] );
    catch exc
        throw( exc );
    end
    x.size = cvx_get_dimlist( x.size );
end

%
% Step 3. Initialize
%

v = cvx( x.size, [] );
if nargout > 0,
    varargout{1} = v;
else
    assignin( 'caller', x.name, v );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% VL_AIBHIST  Compute histogram over VL_AIB tree
%  H = VL_AIBHIST(PARENTS, DATA) computes the histogram of the data
%  points DATA on the VL_AIB tree defined by PARENTS. Each element of
%  DATA indexes one of the leaves of the VL_AIB tree.
%
%  H = VL_AIBHIST(PARENTS, DATA, 'HIST') treats DATA as an histograms.
%  In this case each compoment of DATA is the number of occurences of
%  the VL_AIB leaves corresponding to that component.
%
%  H has the same dimension of parents and counts how many data points
%  are descendent of the corresponding node of the VL_AIB tree.
%
%  See also: VL_HELP(), VL_AIB(), VL_AIBCUTPUSH().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid function
%   g = SIGMOID(z) computes the sigmoid of z.

% You need to return the following variables correctly
g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the sigmoid of each value of z (z can be a matrix,
%               vector or scalar).

g = 1 ./ (1 + exp(-z));

% =============================================================

end

repo: machine-learning


% Copy this template configuration file to your VOT workspace.
% Enter the full path to the ECO repository root folder.

ECO_repo_path = ########

tracker_label = 'ECO';
tracker_command = generate_matlab_command('benchmark_tracker_wrapper(''ECO'', ''VOT2016_DEEP_settings'', true)', {[ECO_repo_path '/VOT_integration/benchmark_wrapper']});
tracker_interpreter = 'matlab';
tracker_trax = false;


text:
function varargout=convfft(varargin)
%
% For calling details please see v_convfft.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_convfft.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: convfft.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_convfft(varargin{:});
else
    v_convfft(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [flowval cut R F] = max_flow(A,u,v,varargin)
% MAX_FLOW Compute the max flow on A from u to v.
%
% flowval=max_flow(A,u,v) computes the maximum flow on the network defined by
% the adjacency structure A, with source u and sink v.
%
% [flowval cut R F] = max_flow(A,u,v) returns the maximum flow in the
% network A with source u and sink v as well as additional information.
% For each vertex on the source side of the mincut, mincut(i) = 1,
% for each vertex on the sink side, mincut(i) = -1.
% R is the residual graph.  R(i,j) is the amount of unused capacity
% on edge (i,j).  F is the flow graph, F(i,j) is the amount of used
% capacity on edge (i,j).  F, A, and R satisfy the relationship A = F + R.
%
% The optional parameter algname specifies the algorithm used to compute
% the maximum flow.  For reference, the push relabel method is likely the
% best general purpose algorithm.  The Edmunds-Karp algorithm
%
% ... = max_flow(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   options.algname: the max flow/min cut algorithm
%     [{'push_relabel'} | 'edmunds_karp' | 'kolmogorov']
%   options.fix_diag: remove any diagonal entries [0 | {1}]
%
% Note: the values on A are interpreted as integers, please round them
% yourself to get the best interpretation.  The code uses the floor of
% the values in A.
%
% Example:
%    load('graphs/max_flow_example.mat')
%    max_flow(A,1,8)

% David Gleich
% Copyright, Stanford University, 2006-2008

%% History
%  2006-04-16: Initial version
%  2006-05-31: Added full2sparse check
%  2007-07-08: Added additional algname
%    Fixed transpose option to implement the pretranspose
%    Fixed documentation bug
%  2007-07-09: Added non-negative edge capacities check
%  2008-09-23: Fixed "check" changing the input (Bug #273796)
%  2008-10-07: Changed options parsing
%    Added fix_diag option
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

options = struct('algname', 'push_relabel','fix_diag',1);
options = merge_options(options, varargin{:});

% no negative capacities and no diagonal entries allowed
if options.fix_diag, A = A - diag(diag(A)); end
if check, check_matlab_bgl(A,struct('noneg',1,'nodiag',1)); end

% max_flow will transpose the data inside

% but ~trans means the input is already transposed, so pre-transpose
if ~trans, A = A'; end

n = size(A,1);

if nargout == 2
    [flowval cut] = max_flow_mex(A,u,v,lower(options.algname));
elseif nargout >= 3
    [flowval cut ri rj rv] = max_flow_mex(A,u,v,lower(options.algname));

    % If anyone needs this operation to be more efficient, send me email,
    % and I can make max_flow_mex return this more efficiently.
    R = sparse(ri,rj,rv,n,n);
    if ~trans
        R = R';
    end
else
    flowval = max_flow_mex(A,u,v,lower(options.algname));
end

if nargout >= 4
    F = A - R;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore