text:
function mimgR = red_channel_mean2(ops)

% numPlanes = length(ops.planesToProcess);

%% build file list with red channel

if (isfield(ops, 'SubDirsRed') && ~isempty(ops.SubDirsRed))
   subDirsRed = ops.SubDirsRed;
else
    if (isfield(ops, 'expred') && ~isempty(ops.expred))
        for i = 1:length(ops.expred)
            subDirsRed{i} = sprintf('%d', ops.expred(i));
        end
    else
        warning('could not find red channel info, returning...')
        return;
    end
end
%%
% build file list
fsrootRED = [];

for j = 1:length(subDirsRed)
    fsrootRED{j} = dir(fullfile(ops.RootDir, subDirsRed{j}, '*.tif'));

    for k = 1:length(fsrootRED{j})
        fsrootRED{j}(k).name = fullfile(ops.RootDir, subDirsRed{j}, fsrootRED{j}(k).name);
    end
end

%%
fsRED = [];
for j = 1:length(fsrootRED)
    fsRED = [fsRED {fsrootRED{j}.name}];
end

%%
D = [];
for k = 1:length(fsRED)
    data = loadFramesBuff(fsRED{k});

    data = cat(3, D, data);
end

%%
clear mimgR
for j = 1:ops.nplanes
    i0 = ops.nchannels_red + (j-1)*ops.nchannels_red;
    %mimgR(:,:,j) = mean(data(:,:,i0:ops.nplanes*ops.nchannels:end), 3);
    datj = data(:,:,i0:ops.nplanes*ops.nchannels_red:end);
    ops1  = AlignIterativeKriging(datj, ops);
    mimgR(:,:,j) = ops1.mimg;
end

%%

author: cortex-lab


text:
function linplus_test2655 ( )

%*****************************************************************************80
%
%% TEST2655 tests R8GB_TO_R8S3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 8;
  ml = 2;
  mu = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST2655\n' );
  fprintf ( 1, '  R8GB_TO_R8S3 copies a R8GB matrix to a R8S3 matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =      %d\n', m );
  fprintf ( 1, '  Matrix columns N =   %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );

  a = r8gb_indicator ( m, n, ml, mu );

  r8gb_print ( m, n, ml, mu, a, '  The R8GB matrix:' );

  nz_num = r8gb_nz_num ( m, n, ml, mu, a );

  fprintf ( 1, '  Nonzeros NZ_NUM =    %d\n', nz_num );

  [ isym, row, col, b ] = r8gb_to_r8s3 ( m, n, ml, mu, a, nz_num );

  r8s3_print ( m, n, nz_num, isym, row, col, b, '  The R8S3 matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function R = functionRlocalscattering(M,theta,ASDdeg,antennaSpacing,distribution)
%Generate the spatial correlation matrix for the local scattering model,
%defined in (2.23) for different angular distributions.
%
%INPUT:
%M              = Number of antennas
%theta          = Nominal angle
%ASDdeg         = Angular standard deviation around the nominal angle
%                 (measured in degrees)
%antennaSpacing = (Optional) Spacing between antennas (in wavelengths)
%distribution   = (Optional) Choose between 'Gaussian', 'Uniform', and
%                'Laplace' angular distribution. Gaussian is default
%
%OUTPUT:
%R              = M x M spatial correlation matrix
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.1 (Last edited: 2017-11-16)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Set the antenna spacing if not specified by input
if  nargin < 4

    %Half a wavelength distance
    antennaSpacing = 1/2;

end

%Set angular distribution to Gaussian if not specified by input
if nargin<5
    distribution = 'Gaussian';
end

%Compute the ASD in radians based on input
ASD = ASDdeg*pi/180;

%The correlation matrix has a Toeplitz structure, so we only need to
%compute the first row of the matrix
firstRow = zeros(M,1);

%Go through all the columns of the first row
for column = 1:M

    %Distance from the first antenna
    distance = antennaSpacing*(column-1);

    %For Gaussian angular distribution
    if strcmp(distribution,'Gaussian')

        %Define integrand of (2.23)
        F = @(Delta)exp(1i*2*pi*distance*sin(theta+Delta)).*exp(-Delta.^2/(2*ASD^2))/(sqrt(2*pi)*ASD);

        %Compute the integral in (2.23) by including 20 standard deviations
        firstRow(column) = integral(F,-20*ASD,20*ASD);

    %For uniform angular distribution
    elseif strcmp(distribution,'Uniform')

        %Set the upper and lower limit of the uniform distribution
        limits = sqrt(3)*ASD;

        %Define integrand of (2.23)
        F = @(Delta)exp(1i*2*pi*distance*sin(theta+Delta))/(2*limits);

        %Compute the integral in (2.23) over the entire interval
        firstRow(column) = integral(F,-limits,limits);

    %For Laplace angular distribution
    elseif strcmp(distribution,'Laplace')

        %Set the scale parameter of the Laplace distribution
        LaplaceScale = ASD/sqrt(2);

        %Define integrand of (2.23)
        F = @(Delta)exp(1i*2*pi*distance*sin(theta+Delta)).*exp(-abs(Delta)/LaplaceScale)/(2*LaplaceScale);

        %Compute the integral in (2.23) by including 20 standard deviations
        firstRow(column) = integral(F,-20*ASD,20*ASD);

    end

end

%Compute the spatial correlation matrix by utilizing the Toeplitz structure
R = toeplitz(firstRow);

author: emilbjornson

repo: massivemimobook


text:
function [K] = kv1u0(x, xp, hyp, ubarp, vbarp, dt, i)

logsigmau = hyp(1);
logthetau = hyp(2);
logsigmav = hyp(3);
logthetav = hyp(4);

a1 = hyp(5);
a2 = hyp(6);

n_x = size(x,1);
n_xp = size(xp,1);

x = repmat(x,1,n_xp);
xp = repmat(xp',n_x,1);

ubarp = repmat(ubarp',n_x,1);
vbarp = repmat(vbarp',n_x,1);

switch i

case 0

K=dt.*exp(1).^(logsigmav+(-2).*logthetav+(-1/2).*exp(1).^((-1).*logthetav) ...
  .*(x+(-1).*xp).^2).*(a2.*exp(1).^(2.*logthetav).*(ubarp.^2+vbarp.^2)+ ...
  a1.*((-1).*exp(1).^logthetav+(x+(-1).*xp).^2));

case 1 % logsigmau

K=0;

case 2 % logthetau

K=0;

case 3 % logsigmav

K=dt.*exp(1).^(logsigmav+(-2).*logthetav+(-1/2).*exp(1).^((-1).*logthetav) ...
  .*(x+(-1).*xp).^2).*(a2.*exp(1).^(2.*logthetav).*(ubarp.^2+vbarp.^2)+ ...
  a1.*((-1).*exp(1).^logthetav+(x+(-1).*xp).^2));

case 4 % logthetav

K=(1/2).*dt.*exp(1).^(logsigmav+(-3).*logthetav+(-1/2).*exp(1).^((-1).* ...
  logthetav).*(x+(-1).*xp).^2).*(a1.*(2.*exp(1).^(2.*logthetav)+(-5).*exp( ...
  1).^logthetav.*(x+(-1).*xp).^2+(x+(-1).*xp).^4)+a2.*exp(1).^(2.* ...
  logthetav).*(ubarp.^2+vbarp.^2).*(x+(-1).*xp).^2);

case 5 % a1

K=dt.*exp(1).^(logsigmav+(-2).*logthetav+(-1/2).*exp(1).^((-1).*logthetav) ...
  .*(x+(-1).*xp).^2).*((-1).*exp(1).^logthetav+(x+(-1).*xp).^2);

case 6 % a2

K=dt.*exp(1).^(logsigmav+(-1/2).*exp(1).^((-1).*logthetav).*(x+(-1).*xp) ...
  .^2).*(ubarp.^2+vbarp.^2);

otherwise

        K = zeros(n_x, n_xp);
end

if K == 0

    K = zeros(n_x, n_xp);

end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function tf = isfield2( S, fs, isinit )
% Similar to isfield but also test whether fields are initialized.
%
% A more comprehensive test of what fields are present [and optionally
% initialized] in a stuct S.  fs is either a single field name or a cell
% array of field names.  The presence of all fields in fs are tested for in
% S, tf is true iif all fs are present. Additionally, if isinit==1, then tf
% is true iff every field fs of every element of S is nonempty (test done
% using isempty).
%
% USAGE
%  tf = isfield2( S, fs, [isinit] )
%
% INPUTS
%  S        - struct array
%  fs       - cell of string name or string
%  isinit   - [0] if true additionally test if all fields are initialized
%
% OUTPUTS
%  tf      - true or false, depending on results of above tests
%
% EXAMPLE
%  isfield2( struct('a',1,'b',2), {'a','b'}, 1 )
%
% See also ISFIELD
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.10
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

tf = all( isfield(S,fs) );
if( ~tf || nargin<3 || ~isinit ); return; end

% now optionally check if fields are initialized
if( iscell(fs) )
  for i=1:length(fs)
    for j=1:numel(S)
      if( isempty(S(j).(fs{i})) ); tf=false; return; end;
    end;
  end;
else
  for j=1:numel(S)
    if( isempty(S(j).(fs)) ); tf=false; return; end;
  end;
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [bdEdge,bdNode,isBdEdge,isBdNode] = findboundaryedge3(edge,elem2edge,bdFlag)
%% FINDBOUNDARYEDGE3 find boundary edges of a three dimensional mesh
%
% [bdEdge,bdNode] = findboundaryedge3(edge,elem2edge,bdFlag) finds all boundary
% edges and boundary nodes of 3D mesh using bdFlag.
%
% For a 2D mesh, use findboundary
%
% See also findboundary3, findboundary
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

% Find boundary edges and nodes
NE = size(edge,1);
N = max(edge(:));
isBdEdge = false(NE,1);
isBdNode = false(N,1);
isBdEdge(elem2edge(bdFlag(:,1) == 1,[4,5,6])) = true;
isBdEdge(elem2edge(bdFlag(:,2) == 1,[2,3,6])) = true;
isBdEdge(elem2edge(bdFlag(:,3) == 1,[1,3,5])) = true;
isBdEdge(elem2edge(bdFlag(:,4) == 1,[1,2,4])) = true;
bdEdge = edge(isBdEdge,:);
isBdNode(bdEdge) = true;
bdNode = find(isBdNode);

repo: ifem


text:
% ========================================================================
% Classification
% USAGE: [prediction, accuracy, err] = classification(D, W, data, Hlabel,
%                                       sparsity)
% Inputs
%       D               -learned dictionary
%       W               -learned classifier parameters
%       data            -testing features
%       Hlabel          -labels matrix for testing feature
%       iterations      -iterations for KSVD
%       sparsity        -sparsity threshold
% outputs
%       prediction      -predicted labels for testing features
%       accuracy        -classification accuracy
%       err             -misclassfication information
%                       [errid featureid groundtruth-label predicted-label]
%
% Author: Zhuolin Jiang (zhuolin@umiacs.umd.edu)
% Date: 10-16-2011
% ========================================================================

function [prediction, accuracy, err] = classification(D, W, data, Hlabel, sparsity)

% sparse coding
G = D'*D;
% Gamma = omp(D'*data,G,sparsity);
Gamma = myOMP(data, D, sparsity);
% classify process
errnum = 0;
err = [];
prediction = [];
for featureid=1:size(data,2)
    spcode = Gamma(:,featureid);
    score_est =  W * spcode;
    score_gt = Hlabel(:,featureid);
    [maxv_est, maxind_est] = max(score_est);  % classifying
    [maxv_gt, maxind_gt] = max(score_gt);
    prediction = [prediction maxind_est];
    if(maxind_est~=maxind_gt)
        errnum = errnum + 1;
        err = [err;errnum featureid maxind_gt maxind_est];
    end
end
accuracy = (size(data,2)-errnum)/size(data,2);

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
% FPMat: FlexPDE Matlab Interface
%
%   Subject:    FlexPDE Matlab Interface
%   Gategory:   Utility
%   Author:     Mohammad Rahmani
%   Created:    Dec. 3, 2004
%   Version:    1.0
%
% Files
%   fpreade  - reads the unformatted data produced by FlexPDE ELEVATION plot command.
%   fpreadh  - reads the unformatted data exported by FlexPDE HISTORY plot command.
%   fpreadsc - reads the unformatted data EXPORTED by FlexPDE on a rectangular grid.
%
% Notes: Tested under Matlab R13, R14 and FlexPDE 4.x.
%
%
% Jan. 13, 2005
% Mohammad Rahmani
% Chemical Eng. Dep.
% Amir-Kabir Uni. of Technology
% Tehran, IRAN
% Mohammad.Rahmani@Gmail.com
%
%

author: Sable

repo: mcbench-benchmarks


text:
function A = suma_stitch(imaname)
% A function to put together a series of suma recorder images
% Very crude, only for the intrepid.
% Example: A = suma_stitch('imageseries_');
%The would be the images as spat out by the recorder window after
% a suma 'r' key stroke with SUMA_SnapshotOverSampling > 1

figure(1); clf
[e, em, lst] = zglobb({sprintf('%s*', imaname)});
N_im = length(lst)
N_1 = round(sqrt(N_im))
k=1
figure(1); clf
for (i=1:1:N_1),
   for (j=1:1:N_1),
      %k = ((i-1)+(j-1)*N_1)+1
      lst(k).name
      a = imread(lst(k).name); size(a)
      if (k==1), A = zeros(N_1.*size(a,1), N_1.*size(a,2), size(a,3), 'uint8'); end
      istrt = (N_1-i)*size(a,1) + 1;
      istp = istrt + size(a,1) - 1;
      jstrt = (j-1)*size(a,2) + 1;
      jstp = jstrt + size(a,2) - 1;
      istrt, istp-istrt, jstrt, jstp-jstrt,
      A(istrt:istp,jstrt:jstp,:) = a;
      subplot (N_1, N_1, k); image(a); title(lst(k).name); drawnow;
      k = k +1;
   end
end

figure(2); clf
image(A); axis square;  drawnow
info = imfinfo(lst(1).name);info.Format
imwrite(A, sprintf('%s_stitch.%s', imaname, info.Format), info.Format);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function mask = mask_create_from_image_set(imgs, outname, atleastN, varargin)
% Take a set of images and create a mask of voxels in which at least N
% subjects have valid (not exactly zero, non NaN) data.
%
% :Usage:
% ::
%
%     mask = mask_create_from_image_set(imgs, outname, atleastN, ['sum'])
%
% This makes a useful results mask for a set of images, i.e., in a
% group analysis.
%
% :Optional: 'sum' input writes the sum image instead of the mask image,
% so that the values in the image reflect the number of input images
% with valid values.
%
% compatible with SPM5 and above only!
%
% :Examples:
% ::
%
%    mask_create_from_image_set(EXPT.SNPM.P{1}, 'mask_all_valid.img');
%
%    imgs = filenames('vascular_mask_*img');
%    mask_create_from_image_set(imgs, 'vascular_group_sum.img', 6, 'sum');
%
% ..
%    Tor Wager, April 2, 2008
%    Edit Sept 13, 2008: Sum image output option
% ..

    write_sum = 0;

    for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                % reserved keywords
                case {'sum', 'sumimage'}, write_sum = 1;

                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
    end

    imgs = strvcat(imgs); % enforce string

    if nargin < 3 || isempty(atleastN), atleastN = size(imgs, 1); end

    spm_defaults

    disp('Mapping images.');
    V = spm_vol(imgs);

    disp('Checking dimensions and spaces.');
    anybad = iimg_check_volinfo(V(1), V);
    if anybad, error('Exiting.'); end

    disp('Reading images.');
    dat = spm_read_vols(V);

    % number of subjects for which each vox is valid
    sumimg = sum(abs(dat) > eps & ~isnan(dat), 4);

    % threshold
    mask = double(sumimg >= atleastN);

    % write output
    if ~isempty(outname)

        outV = struct( 'fname', outname, 'mat', V(1).mat, 'dim', V(1).dim, 'dt', V(1).dt, 'pinfo', [1 0 0]', 'n', [1 V(1).n(2)]  );
        outV.descrip = sprintf('At least %3.0f out of %3.0f images valid', atleastN, size(imgs, 1));

        fprintf('Writing: %s\n', outV.fname);

        if write_sum
            spm_write_vol(outV, sumimg);
        else
            spm_write_vol(outV, mask);
        end
    end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [F,J] = bananaobj(x)
% Evaluate the vector function and the Jacobian matrix for
% the system of nonlinear equations derived from the general
% n-dimensional Rosenbrock function.
% Get the problem size
n = length(x);
if n == 0, error('Input vector, x, is empty.'); end
if mod(n,2) ~= 0,
   error('Input vector, x, must have an even number of components.');
end
% Evaluate the vector function
odds  = 1:2:n;
evens = 2:2:n;
F=x(:);
% F = zeros(n,1);
F(odds,1)  = 1-x(odds);
F(evens,1) = 10.*(x(evens)-x(odds).^2);
% Evaluate the Jacobian matrix if nargout > 1
if nargout > 1
   c = -ones(n/2,1);    C = sparse(odds,odds,c,n,n);
   d = 10*ones(n/2,1);  D = sparse(evens,evens,d,n,n);
   e = -20.*x(odds);    E = sparse(evens,odds,e,n,n);
   J = C + D + E;
end


text:
function [glyphImage]=textImage(varargin)

% function [glyphImage]=textImage(textInput,FontName,FontSize,padAmount)
% ------------------------------------------------------------------------
% This function generates a 2D image containing text as provided in the
% input. The text is of high intensity while the background is of low
% intensity. The font type and font size can be specified allong with the
% amount of padding around the text.
%
%
% ------------------------------------------------------------------------

%%

switch nargin
    case 0
        textInput='GIBBON';
        FontName=vision.internal.getDefaultFont();
        FontSize=50;
        padAmount=0;
    case 1
        textInput=varargin{1};
        FontName=vision.internal.getDefaultFont();
        FontSize=50;
        padAmount=0;
    case 2
        textInput=varargin{1};
        FontName=varargin{2};
        FontSize=50;
        padAmount=0;
    case 3
        textInput=varargin{1};
        FontName=varargin{2};
        FontSize=varargin{3};
        padAmount=0;
    case 4
        textInput=varargin{1};
        FontName=varargin{2};
        FontSize=varargin{3};
        padAmount=varargin{4};
end

position=[0 0];

if isa(textInput,'cell')
    textString=char(textInput);
else
    textString=textInput;
end

if isempty(FontName)
    FontName=vision.internal.getDefaultFont();
end

%%
FONT=listTrueTypeFonts(FontName);

try
[glyphBitmapArray, ...
    glyphIdxFromCharcode, ...
    glyphBitmapStartIdx, ...
    glyphWidths, ...
    glyphHeights, ...
    glyphXAdvances, ...
    glyphLeftBearings, ...
    glyphTopBearings, ...
    fontAscend, ...
    fontDescend, ...
    fontLinespace, ...
    maxBitmapSize]=visionPopulateGlyphBuffer(FONT.fileName,FONT.faceIndex, ...
    FontSize,false);
catch ME
    warning(['Check if ',FontName,' is a valid font use listTrueTypeFonts to see available fonts']);
    rethrow(ME);
end

glyphStruct.glyphBitmapArray     = glyphBitmapArray;
glyphStruct.glyphIdxFromCharcode = glyphIdxFromCharcode;
glyphStruct.glyphBitmapStartIdx  = glyphBitmapStartIdx;
glyphStruct.glyphWidths          = glyphWidths;
glyphStruct.glyphHeights         = glyphHeights;
glyphStruct.glyphXAdvances       = glyphXAdvances;
glyphStruct.glyphLeftBearings    = glyphLeftBearings;
glyphStruct.glyphTopBearings     = glyphTopBearings;
%
glyphStruct.fontAscend            = fontAscend;
glyphStruct.fontDescend           = fontDescend;
glyphStruct.fontLinespace         = fontLinespace;

fontStruct.fontAscend            = fontAscend;
fontStruct.fontDescend           = fontDescend;
fontStruct.fontLinespace         = fontLinespace;

glyphStruct.maxBitmapSize         = maxBitmapSize;

%%

glyphIdxFromCharcode=glyphStruct.glyphIdxFromCharcode;
glyphXAdvances=glyphStruct.glyphXAdvances;
fontHeightWLinegap  = fontStruct.fontLinespace;
fontHeightWOLinegap = fontStruct.fontAscend - fontStruct.fontDescend;

b1 = 1; % for converting 0-based to 1-based
spaceGlyphIdx = glyphIdxFromCharcode(32+b1);
if (spaceGlyphIdx==0)
    spaceCharWidth = int32(fontHeightWOLinegap/4);
else
    spaceCharWidth = int32(glyphXAdvances(spaceGlyphIdx+b1));
end

MARGIN_LeftOrRight = spaceCharWidth; % used 3 before
MARGIN_TopOrBottom = spaceCharWidth; % used 3 before

%%

position = int32(position);

tbLocationXY = int32(position);

tbLocationX = tbLocationXY(1);
tbLocationY = tbLocationXY(2);

tbTopLeftY = tbLocationY;
tbTopLeftX = tbLocationX;

textLocationXY.x = tbTopLeftX + int32(MARGIN_LeftOrRight);
textLocationXY.y = tbTopLeftY + int32(MARGIN_TopOrBottom);

%%

textStringU16 = uint16(textString);

penX = int32(textLocationXY.x);

% go to reference baseline (near the middle of the glyph)
penY = int32(textLocationXY.y) + fontStruct.fontAscend;

oneI32 = int32(1);

isNewLineChar = (textStringU16 == uint16(10));

numChars = numel(textStringU16);

glyphImages=cell(1,numChars);
glyph_I=zeros(numChars,2);
glyph_J=zeros(numChars,2);
q=1;
for ii=1:size(textStringU16,1)
    %Set position
    if ii>1
        penY = penY + fontHeightWLinegap; % go to next line
    end
    penX = textLocationXY.x;
    for jj=1:size(textStringU16,2)

        %see logic in mdlOutputs of sviptextrender.cpp

        % reset x position to the beginning on a line

        thisCharcode = textStringU16(ii,jj);
        thisCharcode_1b = thisCharcode+b1;
        thisGlyphIdx = glyphStruct.glyphIdxFromCharcode(thisCharcode_1b);
        thisGlyphIdx_1b = thisGlyphIdx + b1;
        glyphExists = (thisGlyphIdx ~= 0);
        if ~glyphExists
            penX = penX + int32(spaceCharWidth);
        else
            thisGlyphW = glyphStruct.glyphWidths(thisGlyphIdx_1b);
            thisGlyphH = glyphStruct.glyphHeights(thisGlyphIdx_1b);

            xx=penX+int32(glyphStruct.glyphLeftBearings(thisGlyphIdx_1b));
            yy=penY-int32(glyphStruct.glyphTopBearings(thisGlyphIdx_1b));

            start_I = yy;
            end_I = yy+int32(thisGlyphH)- oneI32;
            start_J = xx;
            end_J = xx+int32(thisGlyphW)- oneI32;

            bitmapStartIdx_1b = glyphStruct.glyphBitmapStartIdx(thisGlyphIdx_1b) + b1;
            bitmapEndIdx_1b =  bitmapStartIdx_1b + uint32(thisGlyphW*thisGlyphH) - uint32(1);
            thisGlyphBitmap = glyphStruct.glyphBitmapArray(bitmapStartIdx_1b:bitmapEndIdx_1b);
            thisGlyphBitmap = reshape(thisGlyphBitmap,[thisGlyphW thisGlyphH])';

            glyphImages{q}=thisGlyphBitmap;
            glyph_I(q,:)=[start_I end_I];
            glyph_J(q,:)=[start_J end_J];

            % update X position for next character
            penX=penX+int32(glyphStruct.glyphXAdvances(thisGlyphIdx_1b));

        end
        q=q+1;
    end
end

%Crop around text
glyph_I=glyph_I-min(glyph_I(:))+1;
glyph_J=glyph_J-min(glyph_J(:))+1;

%% Create image

glyphImage=zeros(max(glyph_I(:))+min(glyph_I(:))-1,max(glyph_J(:))+min(glyph_J(:))-1);
for q=1:1:numChars
    glyphImage(glyph_I(q,1):glyph_I(q,2),glyph_J(q,1):glyph_J(q,2))=glyphImages{q};
end
glyphImage=double(glyphImage);

if padAmount>0
    siz=size(glyphImage)+2*padAmount;
    glyphImagePad=zeros(siz);
    glyphImagePad(padAmount+1:padAmount+size(glyphImage,1),padAmount+1:padAmount+size(glyphImage,2))=glyphImage;
    glyphImage=glyphImagePad;
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function yval = dif_val ( ntab, xtab, diftab, xval )

%*****************************************************************************80
%
%% DIF_VAL evaluates a divided difference polynomial at a point.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl de Boor,
%    A Practical Guide to Splines,
%    Springer Verlag, 1978.
%
%  Parameters:
%
%    Input, integer NTAB, the number of divided difference
%    coefficients in DIFTAB, and the number of points XTAB.
%
%    Input, real XTAB(NTAB), the X values upon which the
%    divided difference polynomial is based.
%
%    Input, real DIFTAB(NTAB), the divided difference
%    polynomial coefficients.
%
%    Input, real XVAL, the value where the polynomial
%    is to be evaluated.
%
%    Output, real YVAL, the value of the polynomial at XVAL.
%
  yval = diftab(ntab);
  for i = 1 : ntab-1
    yval = diftab(ntab-i) + ( xval - xtab(ntab-i) ) * yval;
  end

  return
end


text:
classdef GradientMethod < handle

    properties (Access = private)
        lineSearch
        differentiableFunction
        designVariable
    end

    methods (Access = public)

        function obj = GradientMethod(cParams)
            obj.differentiableFunction = cParams.differentiableFunction;
            obj.designVariable = cParams.designVariable;
            obj.createLineSearch()
        end

        function compute(obj)
            t = obj.lineSearch.value;
            obj.differentiableFunction.computeGradient();
            g = obj.differentiableFunction.gradient;
            x = obj.designVariable.value;
            x = x - t*g;
            obj.designVariable.value = x;
        end

    end

    methods (Access = private)

        function createLineSearch(obj)
            L = obj.differentiableFunction.lipschitzConstant;
            s.value = 1/L;
            obj.lineSearch = ConstantLineSearch(s);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function status = test_secondary_source_selection(modus)
%TEST_SECONDARY_SOURCE_SELECTION tests the correctness of
%secondary_source_selection()
%
%   Usage: status = test_secondary_source_selection(modus)
%
%   Input parameters:
%       modus   - 0: numerical (quiet)
%                 1: visual (not available)
%                 2: numerical verbose
%
%   Output parameters:
%       status - true or false

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

status = false;

%% ===== Checking of input  parameters ===================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Main ============================================================
conf = SFS_config;
conf.secondary_sources.number = 16;
% reference values
ref_selection_circular_pw = ...
    [0   1   1   1   1   1   1   1   0   0   0   0   0   0   0   0]';
ref_selection_circular_ps = ...
    [0   0   1   1   1   1   1   0   0   0   0   0   0   0   0   0]';
ref_selection_circular_fs = ...
    [1   1   1   1   1   0   0   0   0   0   0   0   0   0   0   0]';
ref_selection_linear_pw = ...
    [1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1]';
ref_selection_linear_ps = ...
    [1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1]';
ref_selection_linear_fs = ...
    [0   0   0   0   0   0   0   0   1   1   1   1   1   1   1   1]';
ref_selection_box_pw = ...
    [0   0   0   0   1   1   1   1   0   0   0   0   0   0   0   0]';
ref_selection_box_ps = ...
    [0   0   0   0   1   1   1   1   0   0   0   0   0   0   0   0]';
ref_selection_box_fs = ...
    [0   1   1   1   1   1   1   0   0   0   0   0   0   0   0   0]';
% Calculate current values
% circular array
conf.secondary_source.geometry = 'circular';
x0 = secondary_source_positions(conf);
[~,selection_circular_pw] = secondary_source_selection(x0,[0 -1 0],'pw');
[~,selection_circular_ps] = secondary_source_selection(x0,[0 2.5 0],'ps');
[~,selection_circular_fs] = secondary_source_selection(x0,[0.5 0.5 0 -1 -1 0],'fs');
% linear array
conf.secondary_sources.geometry = 'linear';
x0 = secondary_source_positions(conf);
[~,selection_linear_pw] = secondary_source_selection(x0,[0 -1 0],'pw');
[~,selection_linear_ps] = secondary_source_selection(x0,[0 1 0],'ps');
[~,selection_linear_fs] = secondary_source_selection(x0,[0.5 -0.5 0 -1 -1 0],'fs');
% box form array
conf.secondary_sources.geometry = 'box';
x0 = secondary_source_positions(conf);
[~,selection_box_pw] = secondary_source_selection(x0,[0 -1 0],'pw');
[~,selection_box_ps] = secondary_source_selection(x0,[0 3.5 0],'ps');
[~,selection_box_fs] = secondary_source_selection(x0,[0.5 0.5 0 -1 -1 0],'fs');

if modus==0
    % Numerical mode (quiet)
    if ~all(eq(ref_selection_circular_pw,selection_circular_pw)) || ...
       ~all(eq(ref_selection_circular_ps,selection_circular_ps)) || ...
       ~all(eq(ref_selection_circular_fs,selection_circular_fs)) || ...
       ~all(eq(ref_selection_linear_pw,selection_linear_pw)) || ...
       ~all(eq(ref_selection_linear_ps,selection_linear_ps)) || ...
       ~all(eq(ref_selection_linear_fs,selection_linear_fs)) || ...
       ~all(eq(ref_selection_box_pw,selection_box_pw)) || ...
       ~all(eq(ref_selection_box_ps,selection_box_ps)) || ...
       ~all(eq(ref_selection_box_fs,selection_box_fs))
        return;
    end
elseif modus==2
    message = 'wrong secondary source selection for a';
    if ~all(eq(ref_selection_circular_pw,selection_circular_pw))
        error('%s: %s circular array and a plane wave.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_circular_ps,selection_circular_ps))
        error('%s: %s circular array and a point source.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_circular_fs,selection_circular_fs))
        error('%s: %s circular array and a focused source.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_linear_pw,selection_linear_pw))
        error('%s: %s linear array and a plane wave.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_linear_ps,selection_linear_ps))
        error('%s: %s linear array and a point source.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_linear_fs,selection_linear_fs))
        error('%s: %s linear array anda focused source.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_box_pw,selection_box_pw))
        error('%s: %s box shaped array and a plane wave.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_box_ps,selection_box_ps))
        error('%s: %s box shaped array and a point source.', ...
            upper(mfilename),message);
    end
    if ~all(eq(ref_selection_box_fs,selection_box_fs))
        error('%s: %s box shaped array and a focused source.', ...
            upper(mfilename),message);
    end
end

status = true;

author: sfstoolbox

repo: sfs-matlab


text:
function euler = dcm2euler(DCMbn)
% dcm2euler: converts from body-to-nav DCM to Euler angles.
%
% INPUT
%   DCMbn, 3x3 body-to-nav DCM.
%
% OUTPUT
%   euler, 3x1 Euler angles [roll pitch yaw] (rad, rad, rad).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
% 	Titterton, D.H. and Weston, J.L. (2004). Strapdown
% Inertial Navigation Technology (2nd Ed.). Institution
% of Engineering and Technology, USA. Eq. 11.4. Eq. 3.66, p. 46.
%
%   R. Gonzalez, J. Giribet, and H.D. Patiño,. An approach to
% benchmarking of loosely coupled low-cost navigation systems,
% Mathematical and Computer Modelling of Dynamical Systems, vol. 21,
% issue 2, pp. 272-287. Eq. 15.
%
% Version: 002
% Date:    2016/11/26
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

phi   =  atan( DCMbn(3,2) ./ DCMbn(3,3) );  % roll
theta = -asin( DCMbn(3,1) );                % pitch
psi   =  atan2( DCMbn(2,1), DCMbn(1,1) );   % yaw

euler = [phi theta psi];

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function c_hist = my_color_hist_weight(img, weight)
% img in format (H,W,3)
% weight in format (H,W)

    num_bins = 10;
    c_hist = zeros(num_bins,num_bins,num_bins);
    img = single(img)/256;

    for row = 1:size(img,1)
        for col = 1:size(img,2)
            val = img(row,col,:);
            hist_bin = floor(val*num_bins)+1;
            c_hist(hist_bin(1),hist_bin(2),hist_bin(3)) = c_hist(hist_bin(1),hist_bin(2),hist_bin(3)) + weight(row,col);
        end
    end

%{
    weight = weight/max(weight(:));
    img(:,:,1) = floor(img(:,:,1) .* weight);
    img(:,:,2) = floor(img(:,:,2) .* weight);
    img(:,:,3) = floor(img(:,:,3) .* weight);
    r = histogram(img(:,:,1),100);
    g = histogram(img(:,:,2),100);
    b = histogram(img(:,:,3),100);
  %}
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function json = spm_jsonread(filename, opts)
% JSON (JavaScript Object Notation) parser - a compiled routine
% FORMAT json = spm_jsonread(filename, opts)
% filename - name of a JSON file or JSON string
% json     - JSON structure
% opts     - structure of optional parameters:
%              replacementStyle: string to control how non-alphanumeric
%                characters are replaced {'underscore','hex','delete','nop'}
%                [Default: 'underscore']
%
% References:
%   JSON Standard: http://www.json.org/
%   JSMN C parser: http://zserge.com/jsmn.html
%   jsondecode: http://www.mathworks.com/help/matlab/ref/jsondecode.html
%__________________________________________________________________________
% Copyright (C) 2015-2017 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_jsonread.m 7045 2017-03-17 10:41:12Z guillaume $

%-This is merely the help file for the compiled routine
error('spm_jsonread.c not compiled - see Makefile')


text:
function corrected = apply_cmatrix(im,cmatrix)
% CORRECTED = apply_cmatrix(IM,CMATRIX)
%
% Applies CMATRIX to RGB input IM. Finds the appropriate weighting of the
% old color planes to form the new color planes, equivalent to but much
% more efficient than applying a matrix transformation to each pixel.
if size(im,3)~=3
error('Apply cmatrix to RGB image only.')
end
r = cmatrix(1,1)*im(:,:,1)+cmatrix(1,2)*im(:,:,2)+cmatrix(1,3)*im(:,:,3);
g = cmatrix(2,1)*im(:,:,1)+cmatrix(2,2)*im(:,:,2)+cmatrix(2,3)*im(:,:,3);
b = cmatrix(3,1)*im(:,:,1)+cmatrix(3,2)*im(:,:,2)+cmatrix(3,3)*im(:,:,3);
corrected = cat(3,r,g,b);
end

author: AbdoKamel

repo: sidd-ground-truth-image-estimation


text:
function determ = vand2_determinant ( n, x )

%*****************************************************************************80
%
%% VAND2_DETERMINANT returns the determinant of the VAND2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix desired.
%
%    Input, real X(N), the values that define A.
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  for i = 1 : n
    for j = 1 : i - 1
      determ = determ * ( x(i) - x(j) );
    end
  end

  return
end


text:
function test_prf_full()
%Validate calculation of prf model, including gFit, sFit, hrfFit, final fit
%
%  test_prf()
%
% Tests: ip2volTSeries, rmMain, rmGridFit, rmSearchFit, rmHrfSearchFit
%
% INPUTS
%  No inputs
%
% RETURNS
%  No returns
%
% Example: test_prf_full()
%
% See also MRVTEST
%
% Copyright Stanford team, mrVista, 2015

%% Initialize the key variables and data path
% Data directory (where the mrSession file is located)
dataDir = mrtInstallSampleData('functional','prfInplane');

% This is the validation file
storedPRF = mrtGetValididationData('prfFull');

% These are the items we stored in the validation file
%
% val.roiname    = rmGet(m.model{1}, 'roiname');
% val.eccmn      = nanmean(ecc);
% val.sigmn      = nanmean(sig);
% val.xmax       = max(x);
% val.vemax      = max(ve);
% val.coordsSz   = length(coords);
% save(vFile, '-struct', 'val')

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

% There can be several data types - we will compute PRF model from Averages
dataType = 'Averages';

%% Transfer time series from inplane to gray
ip = initHiddenInplane(); % Foregoes interface - loads data silently

% Set dataTYPE:
ip = viewSet(ip, 'Current DataType', dataType);

% Same for Gray view - initialize hidden view and set dataTYPE
vw = initHiddenGray();
vw = viewSet(vw, 'Current DataType', dataType);

% Transfer time series from inplane to gray
vw = ip2volTSeries(ip,vw,1,'linear'); clear ip;

%% calculate the pRF model
vw = loadROI(vw, 'RV1.mat', 3, [], 0, 1);
vw = rmMain(vw,'RV1.mat' ,5,'min pRF size', 0.5, 'max pRF size', 3,...
    'number of sigmas', 6, 'outerlimit', 0, ...
    'coarse sample', false, 'decimate', 0);
mGauss = load(viewGet(vw, 'rm file'));

coords = rmGet(mGauss.model{1}, 'roiindices');
eccGauss = rmCoordsGet('gray', mGauss.model{1}, 'ecc', coords);
sigGauss = rmCoordsGet('gray', mGauss.model{1}, 'sigma', coords);
xGauss   = rmCoordsGet('gray', mGauss.model{1}, 'x', coords);
veGauss  = rmCoordsGet('gray', mGauss.model{1}, 've', coords);
%% Return to original directory
cd(curDir)

%% Validate the results
tol = 1e-4;

assertEqual(storedPRF.roiname,rmGet(mGauss.model{1}, 'roiname'));

assertEqual(storedPRF.coordsSz, length(coords));

assertElementsAlmostEqual(storedPRF.eccmn, nanmean(eccGauss), 'relative', tol);

assertElementsAlmostEqual(storedPRF.sigmn, nanmean(sigGauss), 'relative', tol);

assertElementsAlmostEqual(storedPRF.xmax, max(xGauss), 'relative', tol);

assertElementsAlmostEqual(storedPRF.vemax, max(veGauss), 'relative', tol);

mrvCleanWorkspace;

%% Test CSS pRF model separately
% [ERK 11/17/21]: While one can request multiple models within one command
% using the variable input 'pRFModel',{'onegaussian','css'}, here we test
% the css model separately as it avoids getting exponents fixed to 1. This
% is because the rmDefineParameters.m function only defines parameters
% based on the first model ('one gaussian', not the second 'css').
vw = rmMain(vw,'RV1.mat' ,5,'min pRF size', 0.5, 'max pRF size', 3,...
    'number of sigmas', 6, 'outerlimit', 0, ...
    'coarse sample', false, 'decimate', 0, 'pRFModel',{'css'});
mCSS = load(viewGet(vw, 'rm file'));

coords = rmGet(mCSS.model{1}, 'roiindices');

% Check if css exponent ranges between 0-1.
expCSS  = rmCoordsGet('gray', mCSS.model{1}, 'exponent', coords);
assert(min(expCSS)>=0);
assert(max(expCSS)<=1);

%% TO DO: Store PRF results and compare against current output %%
eccCSS = rmCoordsGet('gray', mCSS.model{1}, 'ecc', coords);
sigCSS = rmCoordsGet('gray', mCSS.model{1}, 'sigma', coords);
xCSS   = rmCoordsGet('gray', mCSS.model{1}, 'x', coords);
veCSS  = rmCoordsGet('gray', mCSS.model{1}, 've', coords);

%% End Script

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Auto2Fit Test Problems
%http://www.7d-soft.com/en/AF_Examples.htm#5:%20Chart
clear all
clc
load auto2fitData

%% Problem 1.1 - Mount Curve
clc
%Fitting Function
fun = @(m,xdata) m(1)/m(2) * exp(-0.5 * ((xdata-m(3))./m(2)).^2);
%Data
xdata = data11(:,2);
ydata = data11(:,3);

%Problem Setup
opts = optiset('display','iter','solver','nomad');
Opt = opti('fun',fun,'data',xdata,ydata,'x0',[1;3;460],'bounds',[0;0;100],[5;10;1000],'options',opts)

%Solve & Plot
[x,f,e,i] = solve(Opt)
plot(Opt)

%% Problem 1.1 - Mount Curve [multi-solve]
clc
%Fitting Function
fun = @(m,xdata) m(1)/m(2) * exp(-0.5 * ((xdata-m(3))./m(2)).^2);
%Data
xdata = data11(:,2);
ydata = data11(:,3);

%Problem Setup
opts = optiset('display','final','solver','nl2sol');
Opt = opti('fun',fun,'data',xdata,ydata,'bounds',[0;0;100],[5;10;1000],'options',opts)

%Solve & Plot
[x,f,e,i] = multisolve(Opt,[],12)
plot(Opt)

%% Problem 1.2 - Cos Function
clc
%Fitting Function
fun = @(m,xdata) (m(1)./xdata(:,1) - cos(m(2).*xdata(:,2))) .* xdata(:,3)./xdata(:,1) .* m(3);
%Data
xdata = data12(:,2:4);
ydata = data12(:,end);

%Problem Setup
opts = optiset('display','iter','solver','nomad');
Opt = opti('fun',fun,'data',xdata,ydata,'x0',[1;40;1],'bounds',[0;0;0],[5;100;5],'options',opts)

%Solve & Plot
[x,f,e,i] = solve(Opt)
plot(Opt)

%% Problem 1.2 - Cos Function [multi-solve]
clc
%Fitting Function
fun = @(m,xdata) (m(1)./xdata(:,1) - cos(m(2).*xdata(:,2))) .* xdata(:,3)./xdata(:,1) .* m(3);
%Data
xdata = data12(:,2:4);
ydata = data12(:,end);

%Problem Setup
opts = optiset('display','final','solver','nl2sol');
Opt = opti('fun',fun,'data',xdata,ydata,'x0',[1;40;1],'bounds',[0;0;0],[5;100;5],'options',opts)

%Solve & Plot
[x,f,e,i] = multisolve(Opt)
plot(Opt)

%% Problem 1.2 - Cos Function [multi-solve, 2 variable version]
clc
%Fitting Function
fun = @(m,xdata) (2.4923./xdata(:,1) - cos(m(1)*xdata(:,2))) .* xdata(:,3)./xdata(:,1) .* m(2);
%Data
xdata = data12(:,2:4);
ydata = data12(:,end);

%Problem Setup
opts = optiset('display','final','solver','auto');
Opt = opti('fun',fun,'data',xdata,ydata,'bounds',[0;0],[100;5],'options',opts)

%Solve & Plot
[x,f,e,i] = multisolve(Opt,[],[5 3])

figure(1);
plot(Opt)
figure(2);
multiplot(Opt,1,75)

%% Problem 1.3 - SinCos Function (data appears wrong)
% clc
% %Fitting Function
% fun = @(a,xdata) (a(1) + xdata(:,1)./a(2) + cos(a(3) * xdata(:,2)./xdata(:,3)))./(a(4) .* sin(xdata(:,1) + xdata(:,2) + xdata(:,3)));
% %Data
% xdata = data13(:,2:4);
% ydata = data13(:,end);
%
% %Problem Setup
% opts = optiset('display','iter','solver','nomad');
% Opt = opti('fun',fun,'data',xdata,ydata,'x0',[1 1 10 1],'options',opts)
%
% %Solve & Plot
% [x,f,e,i] = solve(Opt)
% plot(Opt)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [L,S,obj,err,iter] = rpca(X,lambda,opts)

% Solve the Robust Principal Component Analysis minimization problem by M-ADMM
%
% min_{L,S} ||L||_*+lambda*loss(S), s.t. X=L+S
% loss(S) = ||S||_1 or ||S||_{2,1}
%
% ---------------------------------------------
% Input:
%       X       -    d*n matrix
%       lambda  -    >0, parameter
%       opts    -    Structure value in Matlab. The fields are
%           opts.loss       -   'l1' (default): loss(S) = ||S||_1
%                               'l21': loss(S) = ||S||_{2,1}
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       L       -    d*n matrix
%       S       -    d*n matrix
%       obj     -    objective function value
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 19/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;
loss = 'l1';

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'loss');        loss = opts.loss;            end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

[d,n] = size(X);

L = zeros(d,n);
S = L;
Y = L;

iter = 0;
for iter = 1 : max_iter
    Lk = L;
    Sk = S;
    % update L
    [L,nuclearnormL] = prox_nuclear(-S+X-Y/mu,1/mu);
    % update S
    if strcmp(loss,'l1')
        S = prox_l1(-L+X-Y/mu,lambda/mu);
    elseif strcmp(loss,'l21')
        S = prox_l21(-L+X-Y/mu,lambda/mu);
    else
        error('not supported loss function');
    end

    dY = L+S-X;
    chgL = max(max(abs(Lk-L)));
    chgS = max(max(abs(Sk-S)));
    chg = max([chgL chgS max(abs(dY(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = nuclearnormL+lambda*comp_loss(S,loss);
            err = norm(dY,'fro');
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y = Y + mu*dY;
    mu = min(rho*mu,max_mu);
end
obj = nuclearnormL+lambda*comp_loss(S,loss);
err = norm(dY,'fro');

function out = comp_loss(E,loss)

switch loss
    case 'l1'
        out = norm(E(:),1);
    case 'l21'
        out = 0;
        for i = 1 : size(E,2)
            out = out + norm(E(:,i));
        end
end

author: canyilu

repo: LibADMM-toolbox


text:
% Test file for @chebfun/realsqrt.m.

function pass = test_realsqrt(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

%% These should fail:

x = chebfun('x', pref);
try
    realsqrt(x);
    pass(1) = false;
catch ME
    pass(1) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:realsqrt:complexRes');
end

try
    realsqrt(1i*abs(x));
    pass(2) = false;
catch ME
    pass(2) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:realsqrt:complexRes');
end

%% These should pass:

realsqrt(x.^2);
pass(3) = true;

realsqrt(abs(x));
pass(4) = true;

%% Test array-valued:

X = [x x];
try
    realsqrt(1i*X);
    pass(5) = false;
catch ME
    pass(5) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:realsqrt:complexRes');
end

try
    realsqrt(X);
    pass(6) = false;
catch ME
    pass(6) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:realsqrt:complexRes');
end

realsqrt(X.^2);
pass(7) = true;

realsqrt(abs(X));
pass(8) = true;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ret = syntheticCalib(noiseType, noiseSigma, calibInfo)
%% Compare four calibration methods using synthetic data.
% Notations
% gt:
%     ground truth
% wb:
%     white board, also called calibration board where checkerboard is
%     attached to the center of it.
% cb: checkerboard
% sl:
%     structured light that is projected by the projector,
%     e.g, slPtsCamImg: structured light points in camera image space
% PR: proposed method
% DG: degraded proposed method (w/o BA)
% GH: generalized global homography method
% MT: Moreno & Taubin method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

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

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%
disp('Calibrating...');

% for debug
verbose = 0;

format short

%% Noise description (noiseType)
% 0. No noise
% 1. Camera image noise
% 2. Projector image noise (before projection).
% 3. Calibration board noise.
% 4. All noises above

%% 1. Load ground truth parameters
% Load calibration data as synthetic params ground truth
gt = Calibration.loadGroundTruth('simulation_ground_truth.yml', calibInfo);

%% 2. Generate synthetic data using gt params
gt = Calibration.generateSyntheticData(gt, noiseType, noiseSigma, verbose);
synthData = gt.synthData;
PR = synthData.PR;
MT = synthData.MT;

%% 3. Calibrate camera using checkerboard points

% camParams is global homography and Moreno & Taubin method's camera
% paramter.
camParams = Calibration.calibrateInitGuess(synthData.cbPtsModel, synthData.cbPtsCamImg, calibInfo);

%% 4. Use calibrated cam params to undistort camera checkerboard and SL nodes

% # of calibration board poses
numSets = length(synthData.cbPtsCamImg);

% checkerboard points in projector image space
cbPtsPrjImgGH = cell(1, numSets);

% undistorted checkerboard and sl grid nodes in camera image space
cbPtsCamImgUndistort = cell(1, numSets);
slPtsCamImgUndistortPR = cell(1, numSets);

for i = 1:numSets
    % undistort checkerboard points in camera image
    cbPtsCamImgUndistort{i} = ImgProc.cvUndistortPoints(synthData.cbPtsCamImg{i}, camParams.camK, camParams.camKc);

    % undistort grid points in camera image
    slPtsCamImgUndistortPR{i} = ImgProc.cvUndistortPoints(PR.slPtsCamImg{i}, camParams.camK, camParams.camKc);

    % find a global homography between camera image space and projector
    % image space
    [Hcp, ~] = cv.findHomography(slPtsCamImgUndistortPR{i}, PR.slPtsPrjImg{i});

    % transform camera checkerboard points to projector image space using
    % global Hcp
    % NOTE: the cb corners have to be undistorted
    cbPtsPrjImgGH{i} = ImgProc.applyHomography(cbPtsCamImgUndistort{i}, Hcp);
end

%% 5. Compute SL nodes in world space and model space for PR and DG

% estimate grid points in world space using camera calibration tvecs and
% rvecs
slPtsModelPR = cell(1, numSets); % estimated grid points in model space
slPtsWorldPR = cell(1, numSets); % estimated grid points in world space

for i = 1:numSets
    % according to Zhang's method the homoraphy from cb model space is:
    % Hmc = lambda*K*[r1, r2, t], where r1, r2 are 1st and 2nd column of R
    R = cv.Rodrigues(camParams.rVecs{i});
    Hmc = camParams.camK * [R(:, 1:2), camParams.tVecs{i}];
    lambda = 1 / norm(inv(camParams.camK) * Hmc(:, 1));
    Hmc = lambda * Hmc;

    % transform camera image grid points to white board model space using
    % Hmc, then use calibrated tvecs and rvecs to transform grid points from
    % model space to world space.
    % NOTE: the cb corners have to be undistorted
    curSlPtsModelPR = ImgProc.applyHomography(slPtsCamImgUndistortPR{i}, inv(Hmc));
    slPtsModelPR{i} = [curSlPtsModelPR, zeros(length(curSlPtsModelPR), 1)];

    % transform grid points from model space to world space.
    slPtsWorldPR{i} = Calibration.rotateTranslatePoints(slPtsModelPR{i}, camParams.tVecs{i}', camParams.rVecs{i}, true);
end

%% 6. [Global Homography] Calibrate using checkerboard points and global homography

prjParamsGH = Calibration.calibrateInitGuess(synthData.cbPtsModel, cbPtsPrjImgGH, calibInfo);

% 2. Stereo calibration using OpenCV
stereoParamsGH = Calibration.calibrateStereoInitGuess(synthData.cbPtsModel, synthData.cbPtsCamImg, cbPtsPrjImgGH, camParams, prjParamsGH, calibInfo);
stereoParamsGH.sets = calibInfo.sets;
stereoParamsGH.dataName = calibInfo.dataName;

% 3. Calculate reprojection errors
stereoParamsGH = Calibration.calcReprojectionError(stereoParamsGH, synthData.cbPtsModel, synthData.cbPtsCamImg, cbPtsPrjImgGH);

if (verbose)
    Reconstruct.visualizePts3d(cell2mat(stereoParamsGH.worldPts'), stereoParamsGH.R, stereoParamsGH.T, 'Global homography calibration: Xm ');
end

%% 7. Subsample SL nodes to save computation time and RAM

for i = 1:length(slPtsModelPR)
    idx = round(linspace(1, length(slPtsModelPR{i}), 500));
    slPtsModelPR{i} = slPtsModelPR{i}(idx, :);
    PR.slPtsCamImg{i} = PR.slPtsCamImg{i}(idx, :);
    PR.slPtsPrjImg{i} = PR.slPtsPrjImg{i}(idx, :);

    % also do the same to slPtsWorldPR
    slPtsWorldPR{i} = slPtsWorldPR{i}(idx, :);
end

%% 8. [Degraded proposed] w/o BA, use DG for short
% 1. Use 'withoutBA' option to specify no bundle adjustment on Xm
stereoParamsDG = Calibration.stereoCalibrate(slPtsModelPR, PR.slPtsCamImg, PR.slPtsPrjImg, calibInfo.camImgSize, calibInfo.prjImgSize, 'withoutBA');
stereoParamsDG.sets = calibInfo.sets;
stereoParamsDG.dataName = calibInfo.dataName;

% 2. Calculate reprojection errors
stereoParamsDG = Calibration.calcReprojectionError(stereoParamsDG, stereoParamsDG.modelPts, stereoParamsDG.camImgPts, stereoParamsDG.prjImgPts);

% debug
if (verbose)
    Reconstruct.visualizePts3d(cell2mat(stereoParamsDG.worldPts'), stereoParamsDG.R, stereoParamsDG.T, 'Propsoed w/o BA calibration: Xm ');
end

%% 9. [proposed] with BA, use PR for short
% 1. Use 'BA' option to specify bundle adjustment on Xm
stereoParamsPR = Calibration.stereoCalibrate(slPtsModelPR, PR.slPtsCamImg, PR.slPtsPrjImg, calibInfo.camImgSize, calibInfo.prjImgSize, 'BA');
stereoParamsPR.sets = calibInfo.sets;
stereoParamsPR.dataName = calibInfo.dataName;

% 2. Calculate reprojection errors
stereoParamsPR = Calibration.calcReprojectionError(stereoParamsPR, stereoParamsPR.modelPts, stereoParamsPR.camImgPts, stereoParamsPR.prjImgPts);

% debug
if (verbose)
    Reconstruct.visualizePts3d(cell2mat(stereoParamsDG.worldPts'), stereoParamsDG.R, stereoParamsDG.T, 'Propsoed calibration: Xm ');
end

%% 10. Compute local homographies and warp corners for Moreno & Taubin method

% same cb corners as GH method
cbPtsModelMT = synthData.cbPtsModel;
cbPtsCamImgMT = synthData.cbPtsCamImg;

% warped checkerboard corners in projector image space (use local homographies)
cbPtsPrjImgMT = [];

% local homography windows size (default parameter)
winSize = 40/2;

% number of cb corners
numCorners = calibInfo.numCorners;

% for each sets (calibration board pose)
for i = 1:numSets
    cbPtsPrjImgMT{i} = zeros(numCorners, 2);

    % undistort grid points in camera image
    curSlPtsCamImgUndistortMT = ImgProc.cvUndistortPoints(MT.slPtsCamImg{i}, camParams.camK, camParams.camKc);
    curSlPtsPrjImgMT = MT.slPtsPrjImg{i};

    % for each cb corner
    for j = 1:numCorners

        % The current checkerboard point
        curCbPtsCamImgRoundMT = round(cbPtsCamImgUndistort{i}(j, :));

        % only use the SL points that are within winSize range of cb point
        % to estimate local homographies
        i1 = curSlPtsCamImgUndistortMT(:, 1) < curCbPtsCamImgRoundMT(1) + winSize;
        i2 = curSlPtsCamImgUndistortMT(:, 1) > curCbPtsCamImgRoundMT(1) - winSize;
        i3 = curSlPtsCamImgUndistortMT(:, 2) < curCbPtsCamImgRoundMT(2) + winSize;
        i4 = curSlPtsCamImgUndistortMT(:, 2) > curCbPtsCamImgRoundMT(2) - winSize;
        inlierIdx = i1 & i2 & i3 & i4;

        if (nnz(inlierIdx) < 4)
            warning('Not enough points to fit a homography, the camera calibration is not very accurate');
            cbPtsPrjImgMT{i}(j, :) = nan;
        else
            % 2d grid points in camera and projector image that are within the window
            curSlPtsCamImgWin = curSlPtsCamImgUndistortMT(inlierIdx, :);
            curSlPtsPrjImgWin = curSlPtsPrjImgMT(inlierIdx, :);

            % calculate local homography
            Hcp = cv.findHomography(curSlPtsCamImgWin, curSlPtsPrjImgWin);

            % transform camera checkerboard corners to projector image
            % space using Hcp
            cbPtsPrjImgMT{i}(j, :) = ImgProc.applyHomography(cbPtsCamImgUndistort{i}(j, :), Hcp);
        end
    end
end

% get rid of nan if fail to estimate a local homogrphy
for i = 1:numSets
    nanIdx = sum(isnan(cbPtsPrjImgMT{i}), 2) > 0;

    cbPtsModelMT{i}(nanIdx, :) = [];
    cbPtsCamImgMT{i}(nanIdx, :) = [];
    cbPtsPrjImgMT{i}(nanIdx, :) = [];
end

%% 11. [Moreno & Taubin method] use MT for short

% 1. Calibrate projector using local homography warped checkerboard corners
prjParamsMT = Calibration.calibrateInitGuess(cbPtsModelMT, cbPtsPrjImgMT, calibInfo);

% 2. Stereo calibration using OpenCV
stereoParamsMT = Calibration.calibrateStereoInitGuess(cbPtsModelMT, cbPtsCamImgMT, cbPtsPrjImgMT, camParams, prjParamsMT, calibInfo);
stereoParamsMT.sets = calibInfo.sets;
stereoParamsMT.dataName = calibInfo.dataName;

% 3. Calculate reprojection errors
stereoParamsMT = Calibration.calcReprojectionError(stereoParamsMT, cbPtsModelMT, cbPtsCamImgMT, cbPtsPrjImgMT);

% debug
if (verbose)
    Reconstruct.visualizePts3d(cell2mat(stereoParamsMT.worldPts'), stereoParamsMT.R, stereoParamsMT.T, 'Moreno & Taubin calibration: Xm ');
end

%% 12. [debug] Compare GH and MT warped checkerboard points in projector image space
if (verbose)
    % MT is + and GH is o
    figure;
    title('warped checkerboard corners in projector image, MT is + and GH is o');
    hold on;
    planeColors = hsv(numSets);

    for i = 1:numSets
        plot(cbPtsPrjImgGH{i}(:, 1), cbPtsPrjImgGH{i}(:, 2), 'o', 'Color', planeColors(i, :));
        plot(cbPtsPrjImgMT{i}(:, 1), cbPtsPrjImgMT{i}(:, 2), '+', 'Color', planeColors(i, :));
    end
end

%% 9. Visualize all error metrics
[errMT, stereoParamsMT] = Calibration.showErrors(stereoParamsMT, gt, 'Moreno & Taubin');
[errGH, stereoParamsGH] = Calibration.showErrors(stereoParamsGH, gt, 'Global homography');
[errDG, stereoParamsDG] = Calibration.showErrors(stereoParamsDG, gt, 'Proposed w/o BA');
[errPR, stereoParamsPR] = Calibration.showErrors(stereoParamsPR, gt, 'Proposed');

%% Convert err to table
ret = [errMT; errGH; errDG; errPR];
resultsTable = struct2table(ret);
disp(resultsTable);

end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function [positions, time] = tracker_affine(video_path, img_files, pos, target_sz, ...
	padding, kernel, lambda, output_sigma_factor, interp_factor, cell_size, ...
	features, show_visualization)
%TRACKER Kernelized/Dual Correlation Filter (KCF/DCF) tracking.
%   This function implements the pipeline for tracking with the KCF (by
%   choosing a non-linear kernel) and DCF (by choosing a linear kernel).
%
%   It is meant to be called by the interface function RUN_TRACKER, which
%   sets up the parameters and loads the video information.
%
%   Parameters:
%     VIDEO_PATH is the location of the image files (must end with a slash
%      '/' or '\').
%     IMG_FILES is a cell array of image file names.
%     POS and TARGET_SZ are the initial position and size of the target
%      (both in format [rows, columns]).
%     PADDING is the additional tracked region, for context, relative to
%      the target size.
%     KERNEL is a struct describing the kernel. The field TYPE must be one
%      of 'gaussian', 'polynomial' or 'linear'. The optional fields SIGMA,
%      POLY_A and POLY_B are the parameters for the Gaussian and Polynomial
%      kernels.
%     OUTPUT_SIGMA_FACTOR is the spatial bandwidth of the regression
%      target, relative to the target size.
%     INTERP_FACTOR is the adaptation rate of the tracker.
%     CELL_SIZE is the number of pixels per cell (must be 1 if using raw
%      pixels).
%     FEATURES is a struct describing the used features (see GET_FEATURES).
%     SHOW_VISUALIZATION will show an interactive video if set to true.
%
%   Outputs:
%    POSITIONS is an Nx2 matrix of target positions over time (in the
%     format [rows, columns]).
%    TIME is the tracker execution time, without video loading/rendering.
%
%   Joao F. Henriques, 2014

	%if the target is large, lower the resolution, we don't need that much
	%detail
	resize_image = (sqrt(prod(target_sz)) >= 100);  %diagonal size >= threshold
	if resize_image,
		pos = floor(pos / 2);
		target_sz = floor(target_sz / 2);
	end

	%window size, taking padding into account
	window_sz = floor(target_sz * (1 + padding));

% 	%we could choose a size that is a power of two, for better FFT
% 	%performance. in practice it is slower, due to the larger window size.
% 	window_sz = 2 .^ nextpow2(window_sz);

	%create regression labels, gaussian shaped, with a bandwidth
	%proportional to target size
	output_sigma = sqrt(prod(target_sz)) * output_sigma_factor / cell_size;
	yf = fft2(gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));

	%store pre-computed cosine window
	cos_window = hann(size(yf,1)) * hann(size(yf,2))';

	if show_visualization,  %create video interface
		update_visualization = show_video(img_files, video_path, resize_image);
	end

	%note: variables ending with 'f' are in the Fourier domain.

	time = 0;  %to calculate FPS
	positions = zeros(numel(img_files), 2);  %to calculate precision

	for frame = 1:numel(img_files),
		%load image
		im = imread([video_path img_files{frame}]);
		if size(im,3) > 1,
			im = rgb2gray(im);
		end
		if resize_image,
			im = imresize(im, 0.5);
		end

		tic()

		if frame == 1
			camera_feat.detector = cv.FeatureDetector('SURF');
			camera_feat.extractor = cv.DescriptorExtractor('SURF');
			camera_feat.matcher = cv.DescriptorMatcher('FlannBased');
			camera_feat.last_keypoints = camera_feat.detector.detect(im);
			camera_feat.last_descriptors = camera_feat.extractor.compute(im, camera_feat.last_keypoints);
		end

		if frame > 1,
			[camera_feat, camera_H] = find_affine(camera_feat, im, pos, target_sz);
			pos = project_t(pos([2, 1]), camera_H);
            pos = pos([2, 1]);
            pos = min([size(im,1),size(im,2)], pos);
            pos = max([0,0], pos);
			%obtain a subwindow for detection at the position from last
			%frame, and convert to Fourier domain (its size is unchanged)
			patch = get_subwindow(im, pos, window_sz);
			zf = fft2(get_features(patch, features, cell_size, cos_window));

			%calculate response of the classifier at all shifts
			switch kernel.type
			case 'gaussian',
				kzf = gaussian_correlation(zf, model_xf, kernel.sigma);
			case 'polynomial',
				kzf = polynomial_correlation(zf, model_xf, kernel.poly_a, kernel.poly_b);
			case 'linear',
				kzf = linear_correlation(zf, model_xf);
			end
			response = real(ifft2(model_alphaf .* kzf));  %equation for fast detection

			%target location is at the maximum response. we must take into
			%account the fact that, if the target doesn't move, the peak
			%will appear at the top-left corner, not at the center (this is
			%discussed in the paper). the responses wrap around cyclically.
			[vert_delta, horiz_delta] = find(response == max(response(:)), 1);
			if vert_delta > size(zf,1) / 2,  %wrap around to negative half-space of vertical axis
				vert_delta = vert_delta - size(zf,1);
			end
			if horiz_delta > size(zf,2) / 2,  %same for horizontal axis
				horiz_delta = horiz_delta - size(zf,2);
			end
			pos = pos + cell_size * [vert_delta - 1, horiz_delta - 1];
		end

		%obtain a subwindow for training at newly estimated target position
		patch = get_subwindow(im, pos, window_sz);
		xf = fft2(get_features(patch, features, cell_size, cos_window));

		%Kernel Ridge Regression, calculate alphas (in Fourier domain)
		switch kernel.type
		case 'gaussian',
			kf = gaussian_correlation(xf, xf, kernel.sigma);
		case 'polynomial',
			kf = polynomial_correlation(xf, xf, kernel.poly_a, kernel.poly_b);
		case 'linear',
			kf = linear_correlation(xf, xf);
		end
		alphaf = yf ./ (kf + lambda);   %equation for fast training

		if frame == 1,  %first frame, train with a single image
			model_alphaf = alphaf;
			model_xf = xf;
		else
			%subsequent frames, interpolate model
			model_alphaf = (1 - interp_factor) * model_alphaf + interp_factor * alphaf;
			model_xf = (1 - interp_factor) * model_xf + interp_factor * xf;
		end

		%save position and timing
		positions(frame,:) = pos;
		time = time + toc();

		%visualization
		if show_visualization,
			box = [pos([2,1]) - target_sz([2,1])/2, target_sz([2,1])];
			stop = update_visualization(frame, box);
			if stop, break, end  %user pressed Esc, stop early

			drawnow
% 			pause(0.05)  %uncomment to run slower
		end

	end

	if resize_image,
		positions = positions * 2;
	end
end

author: flyers

repo: drone-tracking


text:
%% sphere_index
% Below is a demonstration of the features of the |sphere_index| function

%%
clear; close all; clc;

%% Syntax
% |[I J K]=sphere_index(r,Ic,Jc,Kc);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function sF = cos(sF, varargin)
% cost of a function
% Syntax
%   sF = cos(sF)
%   sF = cos(sF, 'bandwidth', bandwidth)
%
% Input
%  sF - @S2FunHarmonic
%
% Output
%  sF - @S2FunHarmonic
%
% Options
%  bandwidth - minimal degree of the spherical harmonic
%

sF = sF.quadrature(@(v) cos(sF.eval(v)),varargin{:});

end

author: mtex-toolbox

repo: mtex

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


text:
function [modeshape,k_term] = TangentialMode(n1,n2,l1,l2,rA,rB,r0A,r0B,k)
%-----------------------------------------------------
% Compute the Tangential Mode in a room.
% Syntax: [modeshape,k_term] =
%          TangentialMode(n1,n2,l1,l2,rA,rB,r0A,r0B,k)
%-----------------------------------------------------

parameter; % Define parameters
tau_m = (3*V)/(5*c*S*bta); % time constant of mth mode
An = sqrt(4); %normalized constant

% Initialisation of variables
modeshape = [];
k_m       = [];
k_term    = [];

% Calculate tangential modes shapes and wavenumber of modes
for n = 1:length(n2)
    modeshape(n,:) = An.*cos(n1.*pi*rA/l1).*cos(n*pi*rB/l2).*...
                     An.*cos(n1.*pi*r0A/l1).*cos(n*pi*r0B/l2);
    k_m(n,:) = sqrt((n1.*pi/l1).^2 + (n*pi/l2).^2);
end;

% Denominator term for Green function
k_term = k^2 - k_m.^2 - i*k/(tau_m*c);


text:
classdef SimpleClient < WebSocketClient
    %CLIENT Summary of this class goes here
    %   Detailed explanation goes here

    properties
    end

    methods
        function obj = SimpleClient(varargin)
            %Constructor
            obj@WebSocketClient(varargin{:});
        end
    end

    methods (Access = protected)
        function onOpen(~,message,~)
            % This function simply displays the message received
            fprintf('%s\n',message);
        end

        function onTextMessage(obj,message)
            % This function simply displays the message received
            fprintf('Message received:\n%s\n',message);
        end

        function onBinaryMessage(obj,bytearray)
            % This function simply displays the message received
            fprintf('Binary message received:\n');
            display(bytearray);
        end

        function onError(obj,message)
            % This function simply displays the message received
            fprintf('Error: %s\n',message);
        end

        function onClose(obj,message)
            % This function simply displays the message received
            fprintf('%s\n',message);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function nu = transport( X, Y, xi )
%TRANSPORT Calculates vector transport
%   NU = TRANSPORT( X, Y, XI ) calculates the vector transport from
%
%         T_X --> X_Y
%
%   where X and Y are given in Tucker format (ttensor) and applies it
%   to the tangent vector XI. XI is decomposed into
%     XI = Y_tilde x U1 x U2 x U3
%          + S x U1_tilde x U2 x U3
%          + S x U1 x U2_tilde x U3
%          + S x U1 x U2 x U3_tilde
%
%   where only Y_tilde, U1_tilde, U2_tilde and U3_tilde are stored
%   as a struct in factorized form.
%
%   See also calcProjection, addFactorized
%

%   GeomCG Tensor Completion. Copyright 2013 by
%   Michael Steinlechner
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

    V1 = Y.U{1}'*X.U{1};
    V2 = Y.U{2}'*X.U{2};
    V3 = Y.U{3}'*X.U{3};

    V1_tilde = Y.U{1}'*xi.U1_tilde;
    V2_tilde = Y.U{2}'*xi.U2_tilde;
    V3_tilde = Y.U{3}'*xi.U3_tilde;

    M1 = ttm( xi.Y_tilde, {V1, V2, V3} );
    M2 = ttm( X.core, {V1_tilde, V2, V3 } );
    M3 = ttm( X.core, {V1, V2_tilde, V3 } );
    M4 = ttm( X.core, {V1, V2, V3_tilde } );

    %first part:
    nu.Y_tilde = M1 + M2 + M3 + M4;

    %second part;

    Y1 = ttm(xi.Y_tilde, {X.U{1}, V2, V3} );
    Y2 = ttm(xi.Y_tilde, {V1, X.U{2}, V3} );
    Y3 = ttm(xi.Y_tilde, {V1, V2, X.U{3}} );

    G1_1 = ttm( X.core, {xi.U1_tilde,   V2,         V3} );
    G1_2 = ttm( X.core, {X.U{1},        V2_tilde,   V3} );
    G1_3 = ttm( X.core, {X.U{1},        V2,         V3_tilde} );

    G2_1 = ttm( X.core, {V1_tilde,  X.U{2},         V3} );
    G2_2 = ttm( X.core, {V1,        xi.U2_tilde,    V3} );
    G2_3 = ttm( X.core, {V1,        X.U{2},         V3_tilde} );

    G3_1 = ttm( X.core, {V1_tilde,  V2,         X.U{3}} );
    G3_2 = ttm( X.core, {V1,        V2_tilde,   X.U{3}} );
    G3_3 = ttm( X.core, {V1,        V2,         xi.U3_tilde} );

    S1_inv = pinv( double( tenmat( Y.core, 1 ) ));
    S2_inv = pinv( double( tenmat( Y.core, 2 ) ));
    S3_inv = pinv( double( tenmat( Y.core, 3 ) ));

    U1_tilde = double( tenmat(Y1 + G1_1 + G1_2 + G1_3, 1) ) * S1_inv;
    U2_tilde = double( tenmat(Y2 + G2_1 + G2_2 + G2_3, 2) ) * S2_inv;
    U3_tilde = double( tenmat(Y3 + G3_1 + G3_2 + G3_3, 3) ) * S3_inv;

    nu.U1_tilde = U1_tilde - Y.U{1} * ( Y.U{1}' * U1_tilde );
    nu.U2_tilde = U2_tilde - Y.U{2} * ( Y.U{2}' * U2_tilde );
    nu.U3_tilde = U3_tilde - Y.U{3} * ( Y.U{3}' * U3_tilde );

end

author: andrewssobral


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function options = get_training_options(epochs,miniBatch,lR,...
    checkpoint_dir,validation_data, validation_frequency, ...
    GPUDevice, L2Reg, checkpoint_period)

if nargin == 6
    L2Reg = 1.000e-7;
    checkpoint_period = 5;
elseif nargin == 7
    checkpoint_period = 5;
end
gpuDevice(GPUDevice);

if exist(checkpoint_dir,'dir') == 0
    mkdir(checkpoint_dir);
end

options = trainingOptions('adam', ...
    'InitialLearnRate',lR, ...
    'L2Regularization',L2Reg, ...
    'MaxEpochs',epochs, ...
    'MiniBatchSize',miniBatch, ...
    'CheckpointPath',checkpoint_dir, ...
    'Shuffle','once',...%'every-epoch', ...
    'VerboseFrequency',1,...
    'LearnRateSchedule','piecewise',...
    'LearnRateDropPeriod',75,... %drop learning rate every 250 epochs
    'LearnRateDropFactor',0.5,... %drop by 0.5
    'Plots','none',...%'training-progress',...
    'ExecutionEnvironment','gpu',...
    'GradientDecayFactor',0.9,...
    'SquaredGradientDecayFactor' ,0.999,...
    'ValidationData',validation_data,...
    'ValidationFrequency',validation_frequency, ...
    'OutputFcn',@(info)plotTrainingAccuracy(info,checkpoint_dir,...
    checkpoint_period));

author: mahmoudnafifi


text:
function [ y, m ] = month_carry_julian ( y, m )

%*****************************************************************************80
%
%% MONTH_CARRY_JULIAN carries a year of months on the Julian calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, the year and month.
%
%    Output, integer Y, integer M, the year and month.
%    On output, M is no greater than 12.
%
  while ( 1 )

    months = year_length_months_julian ( y );

    if ( m <= months )
      break
    end

    m = m - months;
    y = y + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = prior_unif(varargin)
%PRIOR_UNIF  Uniform prior structure
%
%  Description
%    P = PRIOR_UNIF creates uniform prior structure.
%
%  See also
%    PRIOR_*
%
% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_UNIFORM';
  ip.addOptional('p', [], @isstruct);
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Uniform';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Uniform')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  if init
    % set functions
    p.fh.pak = @prior_unif_pak;
    p.fh.unpak = @prior_unif_unpak;
    p.fh.lp = @prior_unif_lp;
    p.fh.lpg = @prior_unif_lpg;
    p.fh.recappend = @prior_unif_recappend;
  end

end

function [w, s, h] = prior_unif_pak(p, w)
  w=[];
  s={};
  h=[];
end

function [p, w] = prior_unif_unpak(p, w)
  w = w;
  p = p;
end

function lp = prior_unif_lp(x, p)
  lp = 0;
end

function lpg = prior_unif_lpg(x, p)
  lpg = zeros(size(x));
end

function rec = prior_unif_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function y = contrast_enhance ( x )

%*****************************************************************************80
%
%% CONTRAST_ENHANCE applies contrast enhancement to a pixel.
%
%  Discussion:
%
%    This function accepts an array representing a pixel
%    neighborhood, and computes a value for the center pixel according to
%    a contrast enhancement algorithm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, uint8 X(N,N), the pixel values.
%
%    Output, double Y, the new value to be assigned to the pixel
%    at X((N+1)/2,(N+1)/2) for enhanced contrast.
%

%
%  Convert image data from UINT8 to DOUBLE.
%
  x = double ( x );
%
%  Get the size of the array.
%
  n = size ( x, 1 );
%
%  Calculate the average value.
%
  x_average = sum ( sum ( x(:,:) ) ) / n / n;
%
%  The constant S should be chosen to be greater than 1;
%
  s = 3.0;
%
%  Compute the new value for the center pixel.
%
  y = ( 1.0 - s ) * x_average + s * x((n+1)/2,(n+1)/2);

  return
end


text:
% lw=load('wandell_fixed_19991208_V3A-V7_dti_leftOccFG+Splenium_wedge_dat');
% lr=load('wandell_fixed_19991208_V3A-V7_dti_leftOccFG+Splenium_ring_dat');
% rw=load('wandell_fixed_19991208_V3A-V7_dti_rightOccFG+Splenium_wedge_dat');
% rr=load('wandell_fixed_19991208_V3A-V7_dti_rightOccFG+Splenium_ring_dat');
mkdir(tempdir, 'matlabFigs');
baseDir = fullfile(tempdir, 'matlabFigs');
disp(['figures will be saved in ' baseDir]);
scanNames = {'ring','wedge','ring'};
scanNum = 2;

coThresh = 0.2;
fontName = 'Helvetica'; %'Comic Sans MS';
fontSize = 18;
leftPhaseLabels  = ['R-HM'; ' UVM'; 'L-HM'; ' LVM'; 'R-HM'];
rightPhaseLabels = ['L-HM'; ' LVM'; 'R-HM'; ' UVM'; 'L-HM'];

volView = getSelectedVolume;
d.roiCoords = getCurROIcoords(volView);
% We get the last coord, which is the left-right axis (lower = left)
% The following will swap coordinate pairs so that the left-most coord
% is always first.
fiberEndptCoordX = [d.roiCoords(3,1:2:end); d.roiCoords(3,2:2:end)];
swapThese = diff(fiberEndptCoordX)<0;
tmp = fiberEndptCoordX(1,swapThese);
fiberEndptCoordX(1,swapThese) = fiberEndptCoordX(2,swapThese);
fiberEndptCoordX(2,swapThese) = tmp;
subName = mrSESSION.subject;
[junk,sessionCode] = fileparts(pwd);
roiName = volView.ROIs(viewGet(volView, 'selectedROI')).name;

fname = fullfile(baseDir, [subName '_' sessionCode '_' roiName '_' scanNames{scanNum}]);
d.co = getCurDataROI(volView,'co',scanNum,d.roiCoords);
d.ph = getCurDataROI(volView,'ph',scanNum,d.roiCoords);

fiberEndptPh = [d.ph(1:2:end); d.ph(2:2:end)];
tmp = fiberEndptPh(1,swapThese);
fiberEndptPh(1,swapThese) = fiberEndptPh(2,swapThese);
fiberEndptPh(2,swapThese) = tmp;
goodVals = isfinite(fiberEndptPh(1,:)) & isfinite(fiberEndptPh(2,:)) & d.co(1:2:end)>coThresh & d.co(2:2:end)>coThresh;
fiberEndptPh = fiberEndptPh(:,goodVals);
% figure;
% subplot(1,2,1); scatter(fiberEndptPh(1,:), fiberEndptPh(2,:));
% title(roiName);
% subplot(1,2,2); hist(fiberEndptPh');

% Get data from entire occ lobe (eg. create a 60mm radius disk ROI
% centered in the calcarine- call the left 'allLeft' and the right
% 'allRight'.)

%Process left hemisphere
volView = getSelectedVolume;
roiNum = strmatch('allleft',lower({volView.ROIs(:).name}));
l.roiCoords = volView.ROIs(roiNum).coords;
l.co = getCurDataROI(volView,'co',scanNum,l.roiCoords);
l.ph = getCurDataROI(volView,'ph',scanNum,l.roiCoords);
allLeftPh = l.ph;
goodVals = isfinite(allLeftPh) & l.co>coThresh;
allLeftPh = allLeftPh(goodVals);
%figure; hist(allLeftPh',20);
cxPhLeft = exp(sqrt(-1)*allLeftPh);
centerPh = mean(cxPhLeft);
mnLeftPh = complexPh2PositiveRad(mean(cxPhLeft));
allLeftPhCentered = angle(cxPhLeft/centerPh);
fiberLeftPhCentered = angle(exp(sqrt(-1)*fiberEndptPh(1,:))/centerPh);

[p,t,df] = statTest(allLeftPhCentered, fiberLeftPhCentered, 't');
figure; set(gcf,'Position', [7, 50, 500, 700]);
subplot(2,1,1); hist(allLeftPhCentered,20);
set(gca, 'fontName', fontName, 'fontSize', fontSize, 'XLim', [-pi, pi], 'XTickLabel','');
pos = get(gca,'YLim');
text(-pi*.9,pos(2)*.9,['Mean phase = ' num2str(mnLeftPh,'%0.1f rad')]);
title(sprintf('t=%0.2f (p=%0.4f, df=%d)', t, p, df));
subplot(2,1,2); hist(fiberLeftPhCentered,20);
set(gca, 'fontName', fontName, 'fontSize', fontSize, 'XLim', [-pi, pi], 'XTickLabel','');
set(gca, 'XLim', [-pi, pi], 'XTick', [-pi, -pi/2, 0, pi/2, pi], 'XTickLabel', leftPhaseLabels);
set(gcf, 'PaperPositionMode', 'auto');
print(gcf, '-dpng', '-r120', [fname '_left.png']);

% Process the right
volView = getSelectedVolume;
roiNum = strmatch('allright',lower({volView.ROIs(:).name}));
r.roiCoords = volView.ROIs(roiNum).coords;
r.co = getCurDataROI(volView,'co',scanNum, r.roiCoords);
r.ph = getCurDataROI(volView,'ph',scanNum, r.roiCoords);
allRightPh = r.ph;
goodVals = isfinite(allRightPh) & r.co>coThresh;
allRightPh = allRightPh(goodVals);
%figure; hist(allRightPh',20);
cxPhRight = exp(sqrt(-1)*allRightPh);
centerPh = mean(cxPhRight);
mnRightPh = complexPh2PositiveRad(mean(cxPhRight));
allRightPhCentered = angle(cxPhRight/centerPh);
fiberRightPhCentered = angle(exp(sqrt(-1)*fiberEndptPh(2,:))/centerPh);

[p,t,df] = statTest(allRightPhCentered, fiberRightPhCentered, 't');
figure; set(gcf,'Position', [7, 50, 500, 700]);
subplot(2,1,1); hist(allRightPhCentered,20);
set(gca, 'fontName', fontName, 'fontSize', fontSize, 'XLim', [-pi, pi], 'XTickLabel','');
pos = get(gca,'YLim');
text(-pi*.9,pos(2)*.9,['Mean phase = ' num2str(mnRightPh,'%0.1f rad')]);
title(sprintf('t=%0.2f (p=%0.4f, df=%d)', t, p, df));
subplot(2,1,2); hist(fiberRightPhCentered,20);
set(gca, 'fontName', fontName, 'fontSize', fontSize, 'XLim', [-pi, pi], 'XTickLabel','');
set(gca, 'XLim', [-pi, pi], 'XTick', [-pi, -pi/2, 0, pi/2, pi], 'XTickLabel', rightPhaseLabels);
set(gcf, 'PaperPositionMode', 'auto');
print(gcf, '-dpng', '-r120', [fname '_right.png']);

save([fname '.mat'], 'd', 'l', 'r', 'coThresh');

% set(gcf, 'PaperPositionMode', 'auto');
% print(gcf, '-dpng', '-r120', '/tmp/bw_right_wedge.png']);
%print(gcf, '-depsc', '-tiff', [fileName '_scatter.eps']);


text:
function geometry_test033 ( )

%*****************************************************************************80
%
%% TEST033 tests LINE_EXP_PERP_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  ntest = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST033\n' );
  fprintf ( 1, '  LINE_EXP_PERP_2D is given an explicit line (P1,P2),\n' );
  fprintf ( 1, '  and another point P3.  It then finds a point\n' );
  fprintf ( 1, '  P4 on (P1,P2) so that (P1,P2) is perpendicular\n' );
  fprintf ( 1, '  to (P3,P4).\n' );

  p1(1:dim_num) = [ 1.0, 3.0 ];
  p2(1:dim_num) = [ 4.0, 0.0 ];

  p3test(1:dim_num,1:ntest) = [ ...
    0.0,  0.0; ...
    5.0, -1.0; ...
    5.0,  3.0 ]';

  r8vec_print ( dim_num, p1, '  Point P1:' );
  r8vec_print ( dim_num, p2, '  Point P2:' );

  for j = 1 : ntest

    p3(1:dim_num) = p3test(1:dim_num,j);

    r8vec_print ( dim_num, p3, '  Point P3:' );

    [ p4, flag ] = line_exp_perp_2d ( p1, p2, p3 );

    r8vec_print ( dim_num, p4, '  Point P4:' );

  end

  return
end


text:
% VL_DEMO_FRAME  Demonstrates vl_plotframe and vl_frame2oell

%% A standard frame is a unit circle centered at the origin with a
% radius pointing in the positive Y direction. Since images use a
% left-handed coordinate system (X axis pointing right, Y axis
% pointing downwards), the radius will point downwards.

A = eye(2) ;
T = [0;0] ;

figure(1) ; clf ;
vl_plotframe([T ; A(:)], 'b', 'linewidth', 5) ; hold on ;

% The same frame can be specified as 4 numbers: translation, scale
% (unit), and rotation (zero).
vl_plotframe([T ; 1 ; 0]) ;

axis equal ; set(gca,'ydir','reverse') ; grid on ;
vl_figaspect(1) ;
vl_demo_print(1, 'frame_standard') ;

%% A positive rotation is appears clockwise in the standard image
% reference frame.

angle = pi/4 ;
R = [cos(angle), -sin(angle) ; sin(angle) cos(angle)] ;
A = R ;

figure(2) ; clf ;
vl_plotframe([T ; A(:)], 'b', 'linewidth', 5) ; hold on ;

% The same frame can be specified as 4 numbers: translation, scale and
% rotation.
vl_plotframe([T ; 1 ; angle]) ;

axis equal ; set(gca,'ydir','reverse') ; grid on ;
vl_figaspect(1) ;
vl_demo_print(2, 'frame_rotated') ;

%% A more complex example, involving a generic affine
% frame. The frame is the result of applying a given affine
% transformation to the standard frame.

angle = pi/4 ;
R = [cos(angle), -sin(angle) ; sin(angle) cos(angle)] ;
S = [1 0 ; 0 .5] ;
A = R*S*R ;

figure(3) ; clf ;
vl_plotframe([T ; A(:)], 'b', 'linewidth', 3) ; hold on ;

axis equal ; set(gca,'ydir','reverse') ; grid on ;
vl_figaspect(1) ;
vl_demo_print(3, 'frame_affine') ;

%% An oriented ellipse can be easily transformed in an unoriented
% one by removing the orientation. The opposite transformation is
% possible as well, provided that a conventional orientation is
% assigned to it.

clear frames1 frames2 ;
randn('state',0) ;
rand('state',0) ;
for i = 1 : 5
  % sample a random ellipse
  angle1 = rand * 2*pi;
  angle2 = rand * 2*pi;
  R1 = [cos(angle1), -sin(angle1) ; sin(angle1) cos(angle1)] ;
  R2 = [cos(angle2), -sin(angle2) ; sin(angle2) cos(angle2)] ;
  A = 0.3*R1*diag(rand(2,1)+0.5)*R2 ;
  frames1(:,i) = [i ; 0 ; A(:)] ;

  % get the unoriented variant of the same ellipse
  S = A*A' ;
  frames2(:,i) = [i ; 0 ; S(1,1) ; S(1,2) ; S(2,2)] ;
end

% Recover oriented ellipses from the unoriented ones; note that
% vl_frame2oell assigns all frames an upright orientation.

frames3 = vl_frame2oell(frames2) ;

frames1(2,:) = 0 ;
frames2(2,:) = 1 ;
frames3(2,:) = 2 ;

figure(4) ; clf ;
vl_plotframe(frames1,'r') ;
vl_plotframe(frames2,'g') ;
vl_plotframe(frames3,'b') ;

axis equal tight ; set(gca,'ydir','reverse') ; grid on ;
vl_figaspect(5/3) ;
vl_demo_print(4, 'frame_frame2oell') ;

repo: StructCompletion


text:
function ph=phase4(cp,mode)
% PHASE4 4 Quadrant arctangent.
%        PHASE4 provides a 4 quadrant arctangent of a complex number such
%        that -360 < PH <= 0.  It also takes care of wrapping when
%        necessary.

% Author: Craig Borghesani
% Date: 8/17/94
% Revised: 10/24/94
% Copyright (c) 1999, Prentice-Hall

ph = atan2(imag(cp),real(cp));

if nargin == 1, % normal phase computation for frequency response plots

 ph = ph - 2*pi*(ph>1e-5);

else % unwrap phase values for gain (phase) plot

 [r,c] = size(ph);
 for k = 1:c,
  dph = diff(ph(:,k));
  loc_brk = find(abs(dph)>pi);
  if length(loc_brk),
   brks = loc_brk;
   if dph(loc_brk(1)) > 0, brks = [0;brks]; end
   if dph(loc_brk(length(loc_brk))) < 0, brks = [brks;length(ph)]; end
   for k2 = 1:2:length(brks),
    ph((brks(k2)+1):(brks(k2+1)),k) = ph((brks(k2)+1):(brks(k2+1)),k) + 2*pi;
   end
  end
 end

end

author: Sable

repo: mcbench-benchmarks


text:
classdef FemTestsSuite < handle

    methods

        function obj = FemTestsSuite()
            path = './tests/Source/FemTests/FemTests.m';
            suite = matlab.unittest.TestSuite.fromFile(path, 'Tag','FEM');
            results = suite.run;
            table(results)
        end

    end
end

author: SwanLab

repo: Swan


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@sym le {(@var{a}, @var{b})}
%% @defopx Operator @@sym {@var{a} <= @var{b}} {}
%% Test/define symbolic inequality, less than or equal to.
%%
%% Examples:
%% @example
%% @group
%% sym(1) <= sym(pi)
%%   @result{} (sym) True
%%
%% syms x
%% x <= 10
%%   @result{} (sym) x ≤ 10
%% @end group
%% @end example
%%
%% @seealso{@@sym/lt, @@sym/gt, @@sym/ge, @@sym/eq, @@sym/ne,
%%          @@sym/logical, @@sym/isAlways}
%% @end defop

function t = le(x, y)

  if (nargin ~= 2)
    print_usage ();
  end

  t = ineq_helper('<=', 'Le', sym(x), sym(y));

end

%!test
%! % simple
%! x = sym(1); y = sym(1); e = x <= y;
%! assert (logical (e))
%! x = sym(1); y = sym(2); e = x <= y;
%! assert (logical (e))

%!test
%! % array -- array
%! syms x
%! a = sym([1 3 3 2*x]);
%! b = sym([2 x 3 10]);
%! e = a <= b;
%! assert (isa (e, 'sym'))
%! assert (logical (e(1)))
%! assert (isa (e(2), 'sym'))
%! assert (isequal (e(2), 3 <= x))
%! assert (logical (e(3)))
%! assert (isa (e(4), 'sym'))
%! assert (isequal (e(4), 2*x <= 10))

repo: octsympy


text:
% Run example 3 to test multi-scale meshing, bound_courant_number, and interp
clearvars; clc

addpath('..')
addpath(genpath('../utilities/'))
addpath(genpath('../datasets/'))
addpath(genpath('../m_map/'))

ERR_TOL = 0.05;
ERR2_TOL = 0.01;
QUAL_TOL = 0.25;
Re2 = 111^2;

bbox = [-71.6 42.7; -64 30; -80 24; -85 38; -71.6 42.7]; %polygon boubox
min_el    = 1e3;  		        % minimum resolution in meters.
max_el    = 50e3; 		        % maximum resolution in meters.
wl        = 30;                 % 60 elements resolve M2 wavelength.
dt        = [50,0.5,6.0];      % Automatically set timestep based on nearshore res
grade     = 0.15;               % mesh grade in decimal percent.
R         = 3; 			        % Number of elements to resolve feature.

dem       = 'SRTM15+.nc';
coastline = 'GSHHS_f_L1';
gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,...
                'bbox',bbox);

fh1 = edgefx('geodata',gdat1,...
             'fs',R,'wl',wl,'max_el',max_el,...
             'dt',dt,'g',grade);

mshopts = meshgen('ef',fh1,'bou',gdat1,...
                  'plot_on',1,'proj','lam');
mshopts = mshopts.build;

m = mshopts.grd;

m = interp(m,gdat1);

%  Set desired courant criteria and time step parameter
dt     = 50; %[s]
max_cr = 4.0;
min_cr = 0.25;
max_it = 20;
m = bound_courant_number(m,dt,max_cr,min_cr,max_it);

Crnew = CalcCFL(m,dt);

if sum(Crnew > max_cr) > 0
    error(['Unable to bound maximum Courant number. Got ',...
        num2str(max(Crnew)),' expecting < ' num2str(max_cr)]);
else
    disp(['Bounded maximum Courant number. Got ',...
        num2str(max(Crnew)),' expecting < ' num2str(max_cr)]);
end
if sum(Crnew < min_cr) > 0
    error(['Unable to bound minimum Courant number. Got ',...
        num2str(min(Crnew)), ' expecting > ' num2str(min_cr)]);
else
    disp(['Bounded minimum Courant number. Got ',...
        num2str(min(Crnew)), ' expecting > ' num2str(min_cr)]);
end

disp('Passed: ECGC');

author: CHLNDDEV


text:
function imgOut = ConvertIPTtoICh(img, inverse)
%
%       imgOut = ConvertIPTtoICh(img, inverse)
%
%
%        Input:
%           -img: image to convert from IPT to ICh or from ICh to IPT.
%           -inverse: takes as values 0 or 1. If it is set to 1 the
%                     transformation from IPT to ICh is applied, otherwise
%                     the transformation from IChto IPT.
%
%        Output:
%           -imgOut: converted image in IPT or ICh (see inverse).
%
%     Copyright (C) 2016  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

check3Color(img);

imgOut = img;

if(inverse == 0)%ICh color space
    imgOut(:,:,2) = sqrt(img(:,:,2).^2 + img(:,:,3).^2);
    imgOut(:,:,3) = atan2(img(:,:,2), img(:,:,3));
else
    imgOut(:,:,2) = sin(img(:,:,3)) .* img(:,:,2);
    imgOut(:,:,3) = cos(img(:,:,3)) .* img(:,:,2);
end

end

author: banterle


text:
function mb_bfgs

% 2010 m.bangert@dkfz.de
% fitting a gaussian distribution to data with a quasi newton method using
% a bfgs or sr1 update for the hessian matrix
%
% example : mb_bfgs (function does not take any argument)

% set options
mode          = 'BFGS'; % update rule for hessian: either 'SR1' or 'BFGS'
visLineSearch = 1;
visBool       = 1;

fprintf(['\nFitting a gaussian distribution to data ...\n' ...
           'Using ' mode ' Hessian update\n']);

% sample data from a gaussian distribution
randn('seed',12345); % set the seed of the random number generator - if you do not want to do that comment...
numOfSamples = 250;
data(:,1)    = linspace(-2,4,numOfSamples) + 6/numOfSamples*randn(1,numOfSamples);
data(:,2)    = exp(-(linspace(-2,4,250)-1.5).^2) + .05*randn(1,250);

% plot data
if visBool
    close all
    figure
    hold on
    plot(data(:,1),data(:,2),'b+')
    xlabel('X')
    ylabel('Y')
    visHandle = [];
end

% specify input
x       = [1 1 1]'; % starting value for parameters
H       = eye(3);   % initial hessian approximation

% 1st calculation of objective function and gradient
objFunc         = @(x) sum((x(3).*exp(-(data(:,1)-x(1)).^2./x(2))-data(:,2)).^2);
objFuncValue    = objFunc(x);
oldObjFuncValue = objFuncValue + 1;
dx              = mb_numDiff(objFunc,x);

% iterate
iter      = 0;
numOfIter = 100;
prec      = 1e-5;

% convergence if gradient smaller than prec, change in objective function
% smaller than prec or maximum number of iteration reached...
while iter < numOfIter && abs((oldObjFuncValue-objFuncValue)/objFuncValue)>prec && norm(dx)>prec

    % increment iteration counter
    iter = iter + 1;

    % remember odl objective function value
    oldObjFuncValue = objFuncValue;

    % implementation of BFGS according to
    % http://en.wikipedia.org/wiki/BFGS_method
    % numbers corresponds

    % 1 obtain search direction
    dir = -(H\dx); % this is the efficient matlab notation for -inv(H)*dx

    % 2.1 linesearch to to find acceptable stepsize alpha
    %alpha = mb_backtrackingLineSearch(objFunc,objFuncValue,x,dx,dir);
    alpha = mb_nocLineSearch(objFunc,@(x) mb_numDiff(objFunc,x),x,dir,dir'*dx,objFuncValue);

    % 3
    p = alpha*dir;

        % visualize the line search (optional) used for debugging
        if visLineSearch
            figure
            hold on
            alphas   = linspace(0,2*alpha,100);
            alphaVis = NaN*alphas;
            for i = 1:100
                alphaVis(i) = objFunc(x+alphas(i)*dir);
            end
            plot(alphas,alphaVis,'b') % function in direction of line search
            plot(alpha,objFunc(x+alpha*dir),'r.') % the steplength found
            plot(alphas,objFuncValue + 1e-1*alphas*(dir'*dx),'g') % constraint
            pause(.5)
            close
        end

    % 2.2 update x
    x = x + p;

    % update objective function (and remember old objective function to
    % check for convergence)
    objFuncValue = objFunc(x);

    % 4 calculate difference of gradients
    dx_old = dx;
    dx     = mb_numDiff(objFunc,x);
    q      = dx-dx_old;

    % update hessian
    if strcmp(mode,'BFGS')
        H = H + (q*q')/(q'*p) - ((H*p)*(H*p)')/(p'*H*p);
    elseif strcmp(mode,'SR1')
        H = H + ((q-H*p)*(q-H*p)')/((q-H*p)'*p);
    end

    % update function with new x and plot
    if visBool
        pause(0.5)
        fitFunc = @(t) x(3).*exp(-(t-x(1)).^2./x(2));
        delete(visHandle);
        visHandle = plot(data(:,1),fitFunc(data(:,1)),'r','LineWidth',4);
        drawnow;
    end

    fprintf(1,'Iteration %d: alpha_min=%f, OF=%f\n',iter,alpha,objFuncValue);

end

fprintf(['\n' num2str(iter) ' iteration(s) performed to converge\n'])

fprintf(1,'Final solution: x=(%f,%f,%f)\n',x(1),x(2),x(3));

end


text:
%function s= multivariate_gauss(x,P,n,randstream)
function s= multivariate_gauss(x,P,n)%, randstream
%
% INPUTS:
%   (x, P) mean vector and covariance matrix
%   obtain n samples
% OUTPUT:
%   sample set
%
% Random sample from multivariate Gaussian distribution.
% Adapted from MVNORMRND (c) 1998, Harvard University.
%
% Tim Bailey 2004.

len = length(x);
S   = chol(P)';
%X   = randn(randstream, len,n);
X   = randn( len,n);
s   = S*X + x*ones(1,n);% Generate values from a normal distribution with mean x
%       and standard deviation S.

repo: Wheel-SLAM


text:
% The mutation function:
function nnum = SMutateBit(num)

% Selecting a position:
pos = ceil(rand*32);

% Flipping that bit:
change = bitset(0,pos);
nnum = bitxor(num,change);


text:
% Group (split) position data by time
%
% This function divides position data into groups based on value of groupParam parameter.
% You can use it to split position data into halves, or into groups of equal time.
%
%  USAGE
%   [groupIndices, edges] = general.groupByTime(pos, groupParam)
%   pos             Position samples. Matrix of size at least Nx3.
%   groupParam      Parameter that defines how groups are formed. Can be string
%                   or digit. The only possible value for string is 'half', which
%                   means that data will be divided into two groups.
%                   Numerical value defines group interval in seconds.
%   groupIndices    Vector of length equal to number of position samples. Each value
%                   indicates group belongingness. If groupIndices(3) == 1, means that
%                   position sample number 3 belongs to group number 1.
%   edges           Vector of edges that define groups. For example, if data is split
%                   in halves, then edges = [0 <middle_point> Inf].
%
%  EXAMPLE
%
%   pos = data.getPositions();
%   groupIndices = general.groupByTime(pos, 'half'); % divide data into halves
%   positionsOfFirstGroup = pos(groupIndices == 1, :);
%
%   general.groupByTime(pos, 1); % divide data into groups each duration of 1 second
%
function [groupIndices, edges] = groupByTime(pos, groupParam)
    if nargin < 2
        error('BNT:numArgs', 'Incorrect number of parameters (type ''help general.groupByTime'' for details).');
    end

    if size(pos, 2) < 3
        error('BNT:arg', 'Incorrect argument ''pos'' (type ''help general.groupByTime'' for details).');
    end

    duration = pos(end, 1) - pos(1, 1);

    if ischar(groupParam)
        if ~strcmpi(groupParam, 'half')
            error('Incorrect argument ''groupParam'' (type ''help general.groupByTime'' for details).');
        end
        midPoint = pos(1, 1) + round(duration / 2);
        edges = [0 midPoint inf];
    else
        grpInterval = groupParam;
        if ~helpers.isdscalar(grpInterval, '>0')
            error('Incorrect argument ''groupParam'' (type ''help general.groupByTime'' for details).');
        end
        edges = [0:grpInterval:duration];
        edges(end+1) = inf;
    end

    [~, groupIndices] = histc(pos(:, 1), edges);
end

author: kavli-ntnu


text:
function so = filterStream(si,filter)
% FILTERSTREAM  Filter (matimatically) a stream
so = {feval(filter,streamCar(si)), ...
        delayEval(@filterStream,{streamCdr(si),filter})};