text:
classdef ShapeContextDistanceExtractor < handle
    %SHAPECONTEXTDISTANCEEXTRACTOR  Implementation of the Shape Context descriptor and matching algorithm
    %
    % Proposed by [Belongie2002]. This implementation is packaged in a generic
    % scheme, in order to allow you the implementation of the common variations
    % of the original pipeline.
    %
    % ## References
    % [Belongie2002]:
    % > Belongie et al. "Shape Matching and Object Recognition using Shape
    % > Contexts" (PAMI 2002)
    %
    % See also: cv.ShapeContextDistanceExtractor.ShapeContextDistanceExtractor,
    %  cv.HausdorffDistanceExtractor, cv.matchShapes
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % The number of angular bins in the shape context descriptor.
        %
        % The number of angular bins for the Shape Context Descriptor used in
        % the shape matching pipeline. default 12
        AngularBins
        % The number of radial bins in the shape context descriptor.
        %
        % The number of radial bins for the Shape Context Descriptor used in
        % the shape matching pipeline. default 4
        RadialBins
        % The inner radius of the shape context descriptor.
        %
        % default 0.2
        InnerRadius
        % The outer radius of the shape context descriptor.
        %
        % default 2
        OuterRadius
        % default false
        RotationInvariant
        % The weight of the shape context distance in the final value of the
        % shape distance.
        %
        % The shape context distance between two shapes is defined as the
        % symmetric sum of shape context matching costs over best matching
        % points. The final value of the shape distance is a user-defined
        % linear combination of the shape context distance, an image
        % appearance distance, and a bending energy. default 1.0
        ShapeContextWeight
        % The weight of the image appearance cost in the final value of the
        % shape distance.
        %
        % The image appearance cost is defined as the sum of squared brightness
        % differences in Gaussian windows around corresponding image points.
        % The final value of the shape distance is a user-defined linear
        % combination of the shape context distance, an image appearance
        % distance, and a bending energy. If this value is set to a number
        % different from 0, is mandatory to set the images that correspond to
        % each shape. default 0.0
        ImageAppearanceWeight
        % The weight of the Bending Energy in the final distance value.
        %
        % The bending energy definition depends on what transformation is being
        % used to align the shapes. The final value of the shape distance is a
        % user-defined linear combination of the shape context distance, an
        % image appearance distance, and a bending energy. default 0.3
        BendingEnergyWeight
        % default 3
        Iterations
        % The standard deviation for the Gaussian window for the image
        % appearance cost.
        %
        % default 10.0
        StdDev
    end

    %% ShapeContextDistanceExtractor
    methods
        function this = ShapeContextDistanceExtractor(varargin)
            %SHAPECONTEXTDISTANCEEXTRACTOR  Constructor
            %
            %     obj = cv.ShapeContextDistanceExtractor()
            %     obj = cv.ShapeContextDistanceExtractor('OptionName',optionValue, ...)
            %
            % ## Options
            % * __AngularBins__ see
            %   cv.ShapeContextDistanceExtractor.AngularBins, default 12
            % * __RadialBins__ see
            %   cv.ShapeContextDistanceExtractor.RadialBins, default 4
            % * __InnerRadius__ see
            %   cv.ShapeContextDistanceExtractor.InnerRadius, default 0.2
            % * __OuterRadius__ see
            %   cv.ShapeContextDistanceExtractor.OuterRadius, default 2
            % * __Iterations__ see
            %   cv.ShapeContextDistanceExtractor.Iterations, default 3
            % * __CostExtractor__ an algorithm that defines the cost matrix
            %   between descriptors, specified as
            %   `{comparerType, 'OptionName',optionValue, ...}`. See
            %   cv.ShapeContextDistanceExtractor.setCostExtractor, where
            %   `comparerType` is one of:
            %   * __NormHistogramCostExtractor__
            %   * __EMDHistogramCostExtractor__
            %   * __ChiHistogramCostExtractor__ (default)
            %   * __EMDL1HistogramCostExtractor__
            % * __TransformAlgorithm__ an algorithm that defines the aligning
            %   transformation, specified as
            %   `{transformerType, 'OptionName',optionValue, ...}`. See
            %   cv.ShapeContextDistanceExtractor.setTransformAlgorithm, where
            %   `transformerType` is one of:
            %   * __ThinPlateSplineShapeTransformer__ (default)
            %   * __AffineTransformer__
            %
            % See also: cv.ShapeContextDistanceExtractor.computeDistance
            %
            this.id = ShapeContextDistanceExtractor_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.ShapeContextDistanceExtractor
            %
            if isempty(this.id), return; end
            ShapeContextDistanceExtractor_(this.id, 'delete');
        end

        function setImages(this, image1, image2)
            %SETIMAGES  Set the images that correspond to each shape, used in the calculation of the Image Appearance cost
            %
            %     obj.setImages(image1, image2)
            %
            % ## Input
            % * __image1__ Image corresponding to the shape defined by
            %   `contours1`.
            % * __image2__ Image corresponding to the shape defined by
            %   `contours2`.
            %
            % See also: cv.ShapeContextDistanceExtractor.getImages
            %
            ShapeContextDistanceExtractor_(this.id, 'setImages', image1, image2);
        end

        function setCostExtractor(this, comparerType, varargin)
            %SETCOSTEXTRACTOR  Set the algorithm used for building the shape context descriptor cost matrix
            %
            %     obj.setCostExtractor(comparerType)
            %     obj.setCostExtractor(comparerType, 'OptionName',optionValue,...)
            %
            % ## Input
            % * __comparerType__ an algorithm that defines the cost matrix
            %   between descriptors. One of:
            %   * __NormHistogramCostExtractor__ A norm based cost extraction.
            %     See cv.norm
            %   * __EMDHistogramCostExtractor__ An EMD based cost extraction.
            %     See cv.EMD
            %   * __ChiHistogramCostExtractor__ An Chi based cost extraction.
            %   * __EMDL1HistogramCostExtractor__ An EMD-L1 based cost
            %     extraction. See cv.EMDL1
            %
            % ## Options
            % The following are options accepted by all algorithms:
            %
            % * __NDummies__ default 25
            % * __DefaultCost__ default 0.2
            %
            % The following are options for the various algorithms:
            %
            % ### `NormHistogramCostExtractor`, `EMDHistogramCostExtractor`
            % * __NormFlag__ default 'L2'. This parameter matches the
            %   `NormType` and `DistType` flags of cv.norm and cv.EMD
            %   respectively.
            %
            % See also: cv.ShapeContextDistanceExtractor.getCostExtractor
            %
            ShapeContextDistanceExtractor_(this.id, 'setCostExtractor', comparerType, varargin{:});
        end

        function setTransformAlgorithm(this, transformerType, varargin)
            %SETTRANSFORMALGORITHM  Set the algorithm used for aligning the shapes
            %
            %     obj.setTransformAlgorithm(transformerType)
            %     obj.setTransformAlgorithm(transformerType, 'OptionName',optionValue,...)
            %
            % ## Input
            % * __transformerType__ an algorithm that defines the aligning
            %   transformation. One of:
            %   * __ThinPlateSplineShapeTransformer__ Definition of the
            %     transformation occupied in the paper [Bookstein89].
            %   * __AffineTransformer__ Wrapper class for the OpenCV Affine
            %     Transformation algorithm. See cv.estimateRigidTransform
            %
            % ## Options
            % The following are options for the various algorithms:
            %
            % ### `ThinPlateSplineShapeTransformer`
            % * __RegularizationParameter__ The regularization parameter for
            %   relaxing the exact interpolation requirements of the TPS
            %   algorithm. default 0
            %
            % ### `AffineTransformer`
            % * __FullAffine__ see cv.estimateRigidTransform, default true
            %
            % ## References
            % [Bookstein89]:
            % > "Principal Warps: Thin-Plate Splines and Decomposition of
            % > Deformations", by F.L. Bookstein (PAMI 1989)
            %
            % See also: cv.ShapeContextDistanceExtractor.getTransformAlgorithm
            %
            ShapeContextDistanceExtractor_(this.id, 'setTransformAlgorithm', transformerType, varargin{:});
        end

        function [image1, image2] = getImages(this)
            %GETIMAGES  Get the images that correspond to each shape, used in the calculation of the Image Appearance cost
            %
            %     [image1, image2] = obj.getImages()
            %
            % ## Output
            % * __image1__ Image corresponding to the shape defined by
            %   `contours1`.
            % * __image2__ Image corresponding to the shape defined by
            %   `contours2`.
            %
            % See also: cv.ShapeContextDistanceExtractor.setImages
            %
            [image1, image2] = ShapeContextDistanceExtractor_(this.id, 'getImages');
        end

        function value = getCostExtractor(this)
            %GETCOSTEXTRACTOR  Get the current algorithm used for building the shape context descriptor cost matrix
            %
            %     value = obj.getCostExtractor()
            %
            % ## Output
            % * __value__ output scalar struct
            %
            % See also: cv.ShapeContextDistanceExtractor.setCostExtractor
            %
            value = ShapeContextDistanceExtractor_(this.id, 'getCostExtractor');
        end

        function value = getTransformAlgorithm(this)
            %GETTRANSFORMALGORITHM  Get the current algorithm used for aligning the shapes
            %
            %     value = obj.getTransformAlgorithm()
            %
            % ## Output
            % * __value__ output scalar struct
            %
            % See also: cv.ShapeContextDistanceExtractor.setTransformAlgorithm
            %
            value = ShapeContextDistanceExtractor_(this.id, 'getTransformAlgorithm');
        end
    end

    %% ShapeDistanceExtractor
    methods
        function dist = computeDistance(this, contour1, contour2)
            %COMPUTEDISTANCE  Compute the shape distance between two shapes defined by its contours
            %
            %     dist = obj.computeDistance(contour1, contour2)
            %
            % ## Options
            % * __contour1__ Contour defining first shape. A numeric
            %   Nx2/Nx1x2/1xNx2 array or a cell-array of 2D points
            %   `{[x,y], ...}`
            % * __contour2__ Contour defining second shape. Same format as
            %   `contours1`.
            %
            % ## Output
            % * __dist__ output distance.
            %
            % See also: cv.ShapeContextDistanceExtractor.ShapeContextDistanceExtractor
            %
            dist = ShapeContextDistanceExtractor_(this.id, 'computeDistance', contour1, contour2);
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.ShapeContextDistanceExtractor.empty,
            %  cv.ShapeContextDistanceExtractor.load
            %
            ShapeContextDistanceExtractor_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.ShapeContextDistanceExtractor.clear,
            %  cv.ShapeContextDistanceExtractor.load
            %
            b = ShapeContextDistanceExtractor_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.ShapeContextDistanceExtractor.load
            %
            ShapeContextDistanceExtractor_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.ShapeContextDistanceExtractor.save
            %
            ShapeContextDistanceExtractor_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.ShapeContextDistanceExtractor.save,
            %  cv.ShapeContextDistanceExtractor.load
            %
            name = ShapeContextDistanceExtractor_(this.id, 'getDefaultName');
        end
    end

    %% Getters/Setters
    methods
        function value = get.AngularBins(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'AngularBins');
        end
        function set.AngularBins(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'AngularBins', value);
        end

        function value = get.RadialBins(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'RadialBins');
        end
        function set.RadialBins(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'RadialBins', value);
        end

        function value = get.InnerRadius(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'InnerRadius');
        end
        function set.InnerRadius(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'InnerRadius', value);
        end

        function value = get.OuterRadius(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'OuterRadius');
        end
        function set.OuterRadius(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'OuterRadius', value);
        end

        function value = get.RotationInvariant(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'RotationInvariant');
        end
        function set.RotationInvariant(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'RotationInvariant', value);
        end

        function value = get.ShapeContextWeight(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'ShapeContextWeight');
        end
        function set.ShapeContextWeight(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'ShapeContextWeight', value);
        end

        function value = get.ImageAppearanceWeight(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'ImageAppearanceWeight');
        end
        function set.ImageAppearanceWeight(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'ImageAppearanceWeight', value);
        end

        function value = get.BendingEnergyWeight(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'BendingEnergyWeight');
        end
        function set.BendingEnergyWeight(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'BendingEnergyWeight', value);
        end

        function value = get.Iterations(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'Iterations');
        end
        function set.Iterations(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'Iterations', value);
        end

        function value = get.StdDev(this)
            value = ShapeContextDistanceExtractor_(this.id, 'get', 'StdDev');
        end
        function set.StdDev(this, value)
            ShapeContextDistanceExtractor_(this.id, 'set', 'StdDev', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [X,nuclearnorm] = prox_nuclear(B,lambda)

% The proximal operator of the nuclear norm of a matrix
%
% min_X lambda*||X||_*+0.5*||X-B||_F^2
%
% version 1.0 - 18/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

[U,S,V] = svd(B,'econ');
S = diag(S);
svp = length(find(S>lambda));
if svp>=1
    S = S(1:svp)-lambda;
    X = U(:,1:svp)*diag(S)*V(:,1:svp)';
    nuclearnorm = sum(S);
else
    X = zeros(size(B));
    nuclearnorm = 0;
end

author: canyilu

repo: LibADMM-toolbox


text:
function mesh2off(OFFfilename, faces,vertices)
% output a mesh into an OFF file
% input:
% OFFfilename is the file name for the off file
% faces is Nx3 matrix for vertex index
% vertices is KxD matrix. D is the number of polygon size. =3 if triangles

file = fopen(OFFfilename,'w');
fprintf (file, 'OFF %d %d 0\n', size(vertices,1), size(faces,1));
for v=1:size(vertices,1)
    fprintf(file, '%f %f %f\n', vertices(v,1), vertices(v,2), vertices(v,3));
end

faces = faces -1; % matlab starts from 1, office starts from 0

for f=1:size(faces,1)
    fprintf(file, '%d', size(faces,2));
    for i=1:size(faces,2)
        fprintf(file, ' %d', faces(f,i));
    end
    fprintf(file, '\n');
end

fclose(file);

author: jianxiongxiao

repo: ProfXkit


text:
function [z,p,sig] = stouffer(p,varargin)
% :Usage:
% ::
%
%     [z,p,sig] = stouffer(p,[alph])
%
% :Inputs:
%
%   **p:**
%        p values in 4-D array
%
%        1st 3 dims are within images, dim4 = image
%
%        optional: alpha value for thresholding
%
% :Outputs:
%
%   **z:**
%        stouffer's combined test statistic, compare to normal
%
%   **p:**
%        p-values for combined test
%
%   **sig:**
%        signficance 1 / 0 binary mask, if alpha is specified
%
% :Described in:
%
% Lazar, N. A., Luna, B., Sweeney, J. A., & Eddy, W. F. (2002).
% Combining brains: a survey of methods for statistical pooling
% of information. Neuroimage, 16(2), 538-550.
%
% Stouffer, S. A., Suchman, E. A., DeVinney, L. C., Star, S. A., and
% Williams, R. M. 1949. The American Soldier: Vol. I. Adjustment
% During Army Life. Princeton University Press, Princeton.
%
% ..
%    tor wager
% ..

if length(varargin) > 0, alph = varargin{1};,else,alph = 0;,end

lastdim = length(size(p));
k = size(p,lastdim);

z = sum(norminv(1 - p),lastdim) ./ sqrt(k);
p = normcdf(1-z);

if alph, sig = p <= alph;,end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function F = acot(F, varargin)
%ACOT   Inverse cotangent of a CHEBFUN.
%   ACOT(F) computes the inverse cotangent of the CHEBFUN F.
%
%   ACOT(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also COT, ACOTD.

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

% Call the compose method:
F = compose(F, @acot, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function runTrackerDets(isSeqDisplay, gtPath, datasetPath, detPath, resPath, seqPath, trackerName)

createPath(resPath); % create the saving path of the tracking results
seqList = load(seqPath); % find the sequence list

% evaluate each sequence
for idSeq = 1:length(seqList)
    if(mod(idSeq,10)==0)
        disp(['tracking the sequence ' num2str(idSeq) '/' num2str(length(seqList)) '...']);
    end
    % load detections and sequence
    seqID = seqList(idSeq);
    seqName = sprintf('%05d',seqID);
    sequence.dataset = dir(fullfile(datasetPath, [seqName '/*.jpg']));
    sequence.seqPath = fullfile(datasetPath, seqName);
    sequence.seqName = seqName;
    img = imread(fullfile(datasetPath, seqName, sequence.dataset(1).name));
    [sequence.imgHeight, sequence.imgWidth, ~] = size(img);
    % nms processing
    detections = [];
    for i = 1:length(sequence.dataset)
        det = load(fullfile(detPath, sprintf('img%03d%03d_loc.txt',seqID,i)));
        if(size(det,2) ~= 3)
            det = det';
        end
        if(~isempty(det))
            idx = det(:,3) > 0;
            det = det(idx,:);
        end
        numdet = size(det,1);
        if(numdet>0)
            curdet = [repmat([i, -1], [numdet, 1]), det(:,1)-10, det(:,2)-10, repmat([20, 20], [numdet, 1]), det(:,3), repmat([-1, -1, -1], [numdet, 1])];
            detections = cat(1, detections, curdet);
        end
    end
    % show the tracking results
    dlmwrite(fullfile(resPath, [seqName '_det.txt']), detections);

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function point = circle3dPoint(circle, pos)
%CIRCLE3DPOINT Coordinates of a point on a 3D circle from its position.
%
%   output = circle3dPoint(input)
%
%   Example
%   % Draw some points on a 3D circle
%     figure('color','w'); hold on; view(130,-10);
%     circle = [10 20 30 50 90 45 0];
%     drawCircle3d(circle)
%     % origin point
%     pos1 = 0;
%     drawPoint3d(circle3dPoint(circle, pos1), 'ro')
%     % few points regularly spaced
%     drawPoint3d(circle3dPoint(circle, 10:10:40), '+')
%     % Draw point opposite to origin
%     drawPoint3d(circle3dPoint(circle, 180), 'k*')
%
%
%   See also
%   circles3d, circle3dPosition
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-21,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

pos=pos(:);

% extract circle coordinates
xc  = circle(1);
yc  = circle(2);
zc  = circle(3);
r   = circle(4);

theta   = circle(5);
phi     = circle(6);
psi     = circle(7);

% convert position to angle
t = pos * pi / 180;

% compute position on base circle
x   = r * cos(t);
y   = r * sin(t);
z   = zeros(length(pos),1);
pt0 = [x y z];

% compute transformation from local basis to world basis
trans   = localToGlobal3d(xc, yc, zc, theta, phi, psi);

% compute points of transformed circle
point   = transformPoint3d(pt0, trans);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ y2, m2, d2, f2 ] = ymdf_inc_roman ( y1, m1, d1, f1, days )

%*****************************************************************************80
%
%% YMDF_INC_ROMAN increments a Roman YMDF date by a DAYS days.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2103
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Input, real DAYS, the number of days to advance the date.
%
%    Output, integer Y2, M2, D2, real F2,
%    the incremented YMDF date.
%

%
%  Copy the parameters.
%
  y2 = y1;
  m2 = m1;
  d2 = d1;
  f2 = f1 + days;
%
%  Check the parameters.
%
  [ y2, m2, d2, f2, ierror ] = ymdf_check_roman ( y2, m2, d2, f2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function handle = inventory(game, player)

%% DISPLAY CARDS
res = get(0, 'ScreenSize');
handle = ...
    figure('Name', 'Inventory', ...
           'Units', 'pixels', ...
           'MenuBar', 'none', ...
           'NumberTitle', 'off', ...
           'Position', [res(3:4)/3, 500, 200]);

Infantry = ...
    axes  ('Parent', handle, ...
           'Position', [0.075, 0.1, 0.25, 0.8]);

Cavallery = ...
    axes  ('Parent', handle, ...
           'Position', [0.375, 0.1, 0.25, 0.8]);

Artillery = ...
    axes  ('Parent', handle, ...
           'Position', [0.675, 0.1, 0.25, 0.8]);

imshow('cards/infantry.jpg', 'Parent', Infantry)
imshow('cards/cavallery.jpg', 'Parent', Cavallery)
imshow('cards/artillery.jpg', 'Parent', Artillery)

%% DISPLAY NUMBERS
ud = get(game.Players(player), 'UserData');
n  = ud(3:7);
N  = [sum(n == 1), sum(n == 2), sum(n == 3)];

for i = 1:3
    uicontrol('style', 'text', ...
              'Parent', handle, ...
              'units', 'normalized', ...
              'FontSize', 18, ...
              'FontWeight', 'bold', ...
              'Position', [0.15 + (i-1) * 0.3, 0.2, 0.1, 0.12], ...
              'BackgroundColor', [245 228 156]/255, ...
              'String', num2str(N(i)));
end

author: Sable

repo: mcbench-benchmarks


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

%%
clear; close all; clc;

%% Syntax
% |L=isclose(A,B,toleranceLevel);|

%% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function region = p00_region ( problem )

%*****************************************************************************80
%
%% P00_REGION returns the name of the integration region for any problem.
%
%  Discussion:
%
%    I thought I was going to use this idea a lot, but most of my test
%    regions are boxes.
%
%    BALL
%      the interior of a 2D circle,
%      the interior of a 3D sphere,
%      the interior of an ND sphere.
%
%    BOX
%      a 1D finite line segment,
%      a 2D finite rectangle,
%      a 3D box,
%      an ND box.
%
%    SIMPLEX
%      a 2D triangle,
%      a 3D tetrahedron,
%      an ND simplex.
%      The "unit simplex" in ND is the set of nonnegative points X
%      such that sum ( X(1:N) ) <= 1.
%
%    SPACE
%      a 1D infinite line,
%      a 2D infinite place,
%      a 3D space,
%      an ND space.
%
%    SPHERE
%      the circumference of a 2D circle,
%      the surface of a 3D sphere,
%      the surface of an ND sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the number of the desired test problem.
%
%    Output, string REGION, the name of the integration region.
%
  if ( problem == 1 )
    region = p01_region ( );
  elseif ( problem == 2 )
    region = p02_region ( );
  elseif ( problem == 3 )
    region = p03_region ( );
  elseif ( problem == 4 )
    region = p04_region ( );
  elseif ( problem == 5 )
    region = p05_region ( );
  elseif ( problem == 6 )
    region = p06_region ( );
  elseif ( problem == 7 )
    region = p07_region ( );
  elseif ( problem == 8 )
    region = p08_region ( );
  elseif ( problem == 9 )
    region = p09_region ( );
  elseif ( problem == 10 )
    region = p10_region ( );
  elseif ( problem == 11 )
    region = p11_region ( );
  elseif ( problem == 12 )
    region = p12_region ( );
  elseif ( problem == 13 )
    region = p13_region ( );
  elseif ( problem == 14 )
    region = p14_region ( );
  elseif ( problem == 15 )
    region = p15_region ( );
  elseif ( problem == 16 )
    region = p16_region ( );
  elseif ( problem == 17 )
    region = p17_region ( );
  elseif ( problem == 18 )
    region = p18_region ( );
  elseif ( problem == 19 )
    region = p19_region ( );
  elseif ( problem == 20 )
    region = p20_region ( );
  elseif ( problem == 21 )
    region = p21_region ( );
  elseif ( problem == 22 )
    region = p22_region ( );
  elseif ( problem == 23 )
    region = p23_region ( );
  elseif ( problem == 24 )
    region = p24_region ( );
  elseif ( problem == 25 )
    region = p25_region ( );
  elseif ( problem == 26 )
    region = p26_region ( );
  elseif ( problem == 27 )
    region = p27_region ( );
  elseif ( problem == 28 )
    region = p28_region ( );
  elseif ( problem == 29 )
    region = p29_region ( );
  elseif ( problem == 30 )
    region = p30_region ( );
  elseif ( problem == 31 )
    region = p31_region ( );
  elseif ( problem == 32 )
    region = p32_region ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_REGION - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', problem );
    error ( 'P00_REGION - Fatal error!' );
  end

  return
end


text:
close all;
clear all;
clc;
rng('default');
N = 256;
% Number of signals
S = 2;
% Choosing non-zero values uniformly between (-b, -a) and (a, b)
a = 1;
b = 2;
% unsigned magnitudes of non-zero entries
XM = a + (b-a).*rand(N, S);
% Generate sign for non-zero entries randomly
sgn = sign(randn(N, S));
% Combine sign and magnitude
X = sgn .* XM;
% Creating noise using helper function
SNR = 15;
Noise = spx.data.noise.Basic.createNoise(X, SNR);
Y = X + Noise;

cs = spx.commons.SignalsComparison(X, Y);
cs.difference_norms()
cs.reference_norms()
cs.estimate_norms()
cs.error_to_signal_norms()
cs.signal_to_noise_ratios()
cs.summarize()

author: indigits

repo: sparse-plex


text:
classdef ComplianceVoigt2TensorConverterPS < FourthOrderVoigt2TensorConverterPS

    properties
    end

    methods (Access = public)

        function obj = ComplianceVoigt2TensorConverterPS(tensor)
           obj.computeConversion(tensor)
        end

    end

    methods (Access = protected,Static)

        function f = getVoigtFactor(iv,jv)
            if ((iv > 2 && jv <= 2) )  || ((iv <= 2 && jv > 2) )
                f = 2;
            elseif (iv>2 && jv>2)
                f = 4;
            else
                f = 1;
            end
        end

    end

    methods (Access = protected)

        function selectTensorClass(obj)
            obj.tensor = CompliancePlaneStressTensor();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function res = indtrial(this, label, flag)
% Method for getting trial indices based on condition labels
% FORMAT res = indtrial(this, label)
% this       - MEEG object
% label      - string or cell array of labels, 'GOOD' and 'BAD'
%              can be added to list of labels to select only
%              good or bad trials respectively
% flag       - 'GOOD' or 'BAD' to include only good or bad trials
%              respectively (all are selected by default)
%
% res        - vector of trial indices matching condition labels
%__________________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: indtrial.m 5059 2012-11-15 13:48:35Z vladimir $

if ischar(label)
    label = {label};
end

[dum, res] = match_str(label, conditions(this));

if nargin > 2
    if strcmpi(flag, 'GOOD')
        res = setdiff(res, badtrials(this));
    elseif strcmpi(flag, 'BAD')
        res = intersect(res, badtrials(this));
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear all; close all; clc

trainingData = csvread('cs-training.csv' , 1 , 1);
training_X = trainingData(:, 2:11); training_y = trainingData(:, 1);

sixCol = find(training_X(:, 5) ~= training_X(:, 5));
NsixCol = find(training_X(: , 5) == training_X(:, 5));

training_X(sixCol , 5) = 1.0 * sum(training_X(NsixCol , 5)) / size(NsixCol , 1);

eleCol = find(training_X(: , 10) ~= training_X(: , 10));
NeleCol = find(training_X(: , 10) == training_X(: , 10));

training_X(eleCol , 10) = 1.0 * sum(training_X(NeleCol , 10)) / size(NeleCol , 1);

[training_X , mu , sigma] = featureNormalize(training_X);

[m, n] = size(training_X);

training_X = [ones(m, 1) training_X];
initial_theta = zeros(n + 1, 1);

lambda = 1;

options = optimset('GradObj', 'on', 'MaxIter', 500);

[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, training_X, training_y, lambda)), initial_theta, options);

p = predict(theta, training_X);

fprintf('Train Accuracy: %f\n', mean(double(p == training_y)) * 100);

%"--------------------Test for TestingData-------------------------------"

testData = csvread('cs-test.csv' , 1 , 2);
test_X = testData;
[tm , tn] = size(test_X);

sixCol = find(test_X(:, 5) ~= test_X(: , 5));
NsixCol = find(test_X(: , 5) == test_X(: , 5));

test_X(sixCol , 5) = 1.0 * sum(test_X(NsixCol , 5)) / size(NsixCol , 1);

eleCol = find(test_X(: , 10) ~= test_X(: , 10));
NeleCol = find(test_X(: , 10) == test_X(: , 10));

test_X(eleCol , 10) = 1.0 * sum(test_X(NeleCol , 10)) / size(NeleCol , 1);

[test_X , mu , sigma] = featureNormalize(test_X);

test_X = [ones(tm , 1) test_X];

prediction = double(sigmoid(test_X * theta));

index = [1:tm]';

dlmwrite('result.csv' , [index , prediction]);

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function shape_test ( code )

%*****************************************************************************80
%
%% SHAPE_TEST verifies the shape function values at the basis nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character ( len = * ) CODE, identifies the element to be used.
%    Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL',
%    'T3', 'T4', 'T6' and 'T10'.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  SHAPE_TEST: Verify shape functions of type "%s"\n', code );

  element_order = order_code ( code );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Element order = %d\n', element_order );

  [ r, s, area ] = node_reference ( code );

  fprintf ( 1, '  Basis function values at basis nodes\n' );
  fprintf ( 1, '  should form the identity matrix.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : element_order
    [ t, dtdr, dtds ] = shape ( code, r(i), s(i) );
    for j = 1 : element_order
      fprintf ( 1, '  %7f', t(j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The R and S derivatives should sum to 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        dTdR sum        dTdS sum\n' );
  fprintf ( 1, '\n' );
  for i = 1 : element_order
    [ t, dtdr, dtds ] = shape ( code, r(i), s(i) );
    rsum = sum ( dtdr(1:element_order) );
    ssum = sum ( dtds(1:element_order) );
    fprintf ( 1, '  %14f  %14f\n', rsum, ssum );
  end

  return
end


text:
function [ G ] = gsp_bunny()
%GSP_BUNNY Create a graph of the stanford bunny
%   Usage :  G = gsp_bunny();
%
%   Output parameters:
%       G           : Resulting graph
%
%   'gsp_bunny()' creates a graph from the pointcloud of the Stanford Bunny model.
%
%   Example:::
%
%           G = gsp_bunny();
%           gsp_plot_graph(G);
%
%   References: turk1994zippered

    %Load the point cloud
    P = gsp_pointcloud('bunny');

    %Create the graph from the point cloud using an epsilon-neighborhood
    %connectivity
    param.type = 'knn';
    param.rescale = 1;
    param.center = 1;

    %Compute it
    G = gsp_nn_graph(double(P), param);
    %Reduce vertex size for plotting
    G.plotting.vertex_size = 10;
end

repo: gspbox


text:
function [net_conv, net_fc, opts] = mdnet_init(image, net)
% MDNET_INIT
% Initialize MDNet tracker.
%
% Hyeonseob Nam, 2015
%

%% set opts
% use gpu
opts.useGpu = true;

% model def
opts.net_file = net;

% test policy
opts.batchSize_test = 256; % <- reduce it in case of out of gpu memory

% bounding box regression
opts.bbreg = true;
opts.bbreg_nSamples = 1000;

% learning policy
opts.batchSize = 128;
opts.batch_pos = 32;
opts.batch_neg = 96;

% initial training policy
opts.learningRate_init = 0.0001; % x10 for fc6
opts.maxiter_init = 30;

opts.nPos_init = 500;
opts.nNeg_init = 5000;
opts.posThr_init = 0.7;
opts.negThr_init = 0.5;

% update policy
opts.learningRate_update = 0.0003; % x10 for fc6
opts.maxiter_update = 10;

opts.nPos_update = 50;
opts.nNeg_update = 200;
opts.posThr_update = 0.7;
opts.negThr_update = 0.3;

opts.update_interval = 10; % interval for long-term update

% data gathering policy
opts.nFrames_long = 100; % long-term period
opts.nFrames_short = 20; % short-term period

% cropping policy
opts.input_size = 107;
opts.crop_mode = 'wrap';
opts.crop_padding = 16;

% scaling policy
opts.scale_factor = 1.05;

% sampling policy
opts.nSamples = 256;
opts.trans_f = 0.6; % translation std: mean(width,height)*trans_f/2
opts.scale_f = 1; % scaling std: scale_factor^(scale_f/2)

% set image size
opts.imgSize = size(image);

%% load net
net = load(opts.net_file);
if isfield(net,'net'), net = net.net; end
net_conv.layers = net.layers(1:10);
net_fc.layers = net.layers(11:end);
clear net;

for i=1:numel(net_fc.layers)
    switch (net_fc.layers{i}.name)
        case {'fc4','fc5'}
            net_fc.layers{i}.filtersLearningRate = 1;
            net_fc.layers{i}.biasesLearningRate = 2;
        case {'fc6'}
            net_fc.layers{i}.filtersLearningRate = 10;
            net_fc.layers{i}.biasesLearningRate = 20;
    end
end

if opts.useGpu
    net_conv = vl_simplenn_move(net_conv, 'gpu') ;
    net_fc = vl_simplenn_move(net_fc, 'gpu') ;
else
    net_conv = vl_simplenn_move(net_conv, 'cpu') ;
    net_fc = vl_simplenn_move(net_fc, 'cpu') ;
end

end

author: flyers

repo: drone-tracking


text:
function amin_index = i4vec_amin_index ( n, a )

%*****************************************************************************80
%
%% I4VEC_AMIN_INDEX returns the index of the smallest magnitude in an I4VEC.
%
%  Discussion:
%
%    An I4VEC is a vector of integer values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries to be checked.
%
%    Input, integer A(N), the vector to be checked.
%
%    Output, integer AMIN_INDEX, the entry of the smallest magnitude.
%
  if ( n <= 0 )

    amin_index = 0;

  else

    aamin = a(1);
    amin_index = 1;

    for i = 2 : n

      if ( abs ( a(i) ) < aamin )
        aamin = abs ( a(i) );
        amin_index = i;
      end

    end

  end

  return
end


text:
 function [I] = wwwimread(url)

 %create a temporary file name with legal characters
 if 0
 tmpdir = '/tmp/';
 tmpfile = url;
 tmpfile = strrep(tmpfile,'//','');
 tmpfile = strrep(tmpfile,'/','-');
 tmpfile = strrep(tmpfile,'.','-');
 tmpfile = strrep(tmpfile,'http:',tmpdir);%run wget
 command = sprintf('wget %s -O %s',url,tmpfile);
 [s,w] = unix(command);
 end

 url
 s = urlread(url);
 size(s)
 fid = fopen(['./tmp.' url(end-2:end)], 'w');
 fwrite(fid, s);
 fclose(fid);

 I = imread(['./tmp.' url(end-2:end)]);
 delete(['./tmp.' url(end-2:end)]);

 if 0
 if (s == 0)
   %read in the file then delete the tmpfile
   I = imread(tmpfile);
   command = sprintf('rm -f %s',tmpfile);
   [s,w] = unix(command);
 else
   disp('error!');
   I = [];
 end;

end

author: zouchuhang

repo: LayoutNet


text:
function [pnt, dhk] = read_off(fn);

% READ_OFF reads vertices and triangles from a OFF format triangulation file
%
% [pnt, dhk] = read_off(filename)
%
% See also READ_TRI, READ_BND

% Copyright (C) 1998, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_off.m 945 2010-04-21 17:41:20Z roboos $

fid = fopen(fn, 'rt');
if fid~=-1

  % scan the file type
  [s, count] = fscanf(fid, '%s\n', 1);
  if ~strcmp(s,'OFF')
    msg = sprintf('wrong file type %s', s);
    error(msg);
  end

  % read the number of vertex points and triangles
  [val, count] = fscanf(fid, '%d', 3);
  Npnt = val(1)
  Ndhk = val(2)

  % read the vertex points
  pnt  = fscanf(fid, '%f', [3, Npnt]);
  pnt  = pnt(1:3,:)';

  % read the triangles
  dhk = fscanf(fid, '%d', [4, Ndhk]);
  dhk = (dhk(2:4,:)+1)';
  fclose(fid);

else
  error('unable to open file');
end

author: PatternRecognition

repo: OpenBMI


text:
%% FUNCTION LSSMTC
%   compute the mutual information
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Quanquan Gu, Jiayu Zhou, and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 17, 2012.
%

function MIhat = MutualInfo(L1,L2)
%   mutual information

%===========
L1 = L1(:);
L2 = L2(:);
if size(L1) ~= size(L2)
    error('size(L1) must == size(L2)');
end
L1 = L1 - min(L1) + 1;      %   min (L1) <- 1;
L2 = L2 - min(L2) + 1;      %   min (L2) <- 1;
%===========    make bipartition graph  ============
nClass = max(max(L1), max(L2));
G = zeros(nClass);
for i=1:nClass
    for j=1:nClass
        G(i,j) = length(find(L1 == i & L2 == j));
    end
end
sumG = sum(G(:));
%===========    calculate MIhat
P1 = sum(G,2);  P1 = P1/sumG;
P2 = sum(G,1);  P2 = P2/sumG;
H1 = sum(-P1.*log2(P1));
H2 = sum(-P2.*log2(P2));
P12 = G/sumG;
PPP = P12./repmat(P2,nClass,1)./repmat(P1,1,nClass);
PPP(abs(PPP) < 1e-12) = 1;
MI = sum(P12(:) .* log2(PPP(:)));
MIhat = MI / max(H1,H2);

MIhat = real(MIhat);

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [predict_label,accuracy] = jdknn(x_train,y_train,x_test,y_test,k,distance,rule)
%%jdknn:packaged knnclassify from Matlab,just for easier use.
%%==============================================================================
%%input:
%%------x_train,...,y_test  :   training and testing sets.                               [required]
%%------k                   :   number of nearest neighbors,default is 5.    [not required]
%%------distance            :   distance metric,choices are:                 [not required]
%%%%-----------euclidean,cityblock,cosine,correlation,Hamming(default is euclidean)
%%------rule                :rules used,choices are:                         [not required]
%%%%-----------nearest,random,consensus(default is nearest)

%%output:
%%------predict_label       :   predicted label vector for test case.
%%------accuracy            :   accuracy
%%==============================================================================
	predict_label = [];
    accuracy = 0;
    k_default = 5;
	distance_default = 'euclidean';
	rule_default = 'nearest';

	k_para = k_default;
	distance_para = distance_default;
	rule_para = rule_default;
    switch nargin
        case 4
            k_para = k_default;
            distance_para = distance_default;
            rule_para = rule_default;
        case 5
            k_para = k;
            distance_para = distance_default;
            rule_para = rule_default;
        case 6
            k_para = k;
            distance_para = distance;
            rule_para = rule_default;
        case 7
            k_para = k;
            distance_para = distance;
            rule_para = rule;
        otherwise
            fprintf('++++++Fatal error!Please check the input!');
            return
    end
    knn_model = fitcknn(x_train,y_train,'NumNeighbors',k_para);

	predict_label = knn_model.predict(x_test);
	accuracy = length(find(predict_label == y_test)) / length(y_test);
end

author: jindongwang

repo: activityrecognition


text:
%[ANR,SNR,STR]	=  ASORT(INP,'OPT',...);
% S		=  ASORT(INP,'OPT',...);
%		   to sort alphanumeric strings numerically if
%		   they contain one properly formatted number
%		   otherwise, ascii dictionary sorting is applied
%
% INP	unsorted input:
%	- a char array
%	- a cell array of strings
% OPT	options
%  -s	- sorting option
%	  '-s','ascend'					[def]
%	  '-s','descend'
%  -st	- force output form S				[def: nargout dependent]
%  -t	- replace matching template(s) with one space
%	  prior to sorting
%	  '-t','template'
%	  '-t',{'template1','template2',...}
%  -w	- remove space(s) prior to sorting
%
%	  NOTE	-t/-w options are processed in the
%		      order that they appear in
%		      the command line
%
%  -v	- verbose output				[def: quiet]
%  -d	- debug mode
%	  save additional output in S
%	  .c:	lex parser input
%	  .t:	lex parser table
%	  .n:	lex parser output
%	  .d:	numbers read from .n
%
% ANR	numerically sorted alphanumeric strings		[eg, 'f.-1.5e+2x.x']
%	- contain one number that can be read by
%	  <strread> | <sscanf>
% SNR	ascii dict  sorted alphanumeric strings
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=7212#
%
%	- contain more than one number			[eg, 'f.-1.5e +2.x']
%	- contain incomplete|ambiguous numbers		[eg, 'f.-1.5e+2.x']
% STR	ascii dict  sorted strings
%	- contain no numbers				[eg, 'a test']
%
% S	structure with fields
%	.anr
%	.srn
%	.str

% created:
%	us	03-Mar-2002
% modified:
%	us	30-Mar-2005 11:57:07 	/ TMW R14.sp2

%--------------------------------------------------------------------------------
function	varargout=asort(inp,varargin)

varargout(1:nargout)={[]};
if	~nargin
	help(mfilename);
	return;
end

% - common parameters/options
n=[];
ds=[];
anr={};
snr={};
str={};
smod='ascend';	% sorting option
tmpl={};	% template(s)
sflg=false;	% output  mode: structure
tflg=false;	% remove  template(s)
dflg=false;	% debug   mode
vflg=false;	% verbose output
wflg=false;	% remove  spaces

if	nargin > 1
	ix=find(strcmp('-s',varargin));
	if	~isempty(ix) && nargin > ix(end)+1
		smod=varargin{ix(end)+1};
	end
	ix=find(strcmp('-t',varargin));
	if	~isempty(ix) && nargin > ix(end)+1
		tflg=ix(end);
		tmpl=varargin{ix(end)+1};
	end
	if	find(strcmp('-d',varargin));
		dflg=true;
	end
	if	find(strcmp('-st',varargin));
		sflg=true;
	end
	if	find(strcmp('-v',varargin));
		vflg=true;
	end
	ix=find(strcmp('-w',varargin));
	if	~isempty(ix)
		wflg=ix(end);
	end
end
%   spec numbers
ntmpl={
	' inf '
	'+inf '
	'-inf '
	' nan '
	'+nan '
	'-nan '
	};
%   spec chars
ctmpl={
	'.'	% decimal point
	'd'	% exponent
	'e'	% exponent
	};

if	nargout <= 3
	varargout{1}=inp;
else
	disp(sprintf('ASORT> too many output args [%-1d/%-1d]\n',nargout,3));
	help(mfilename);
	return;
end
if	isempty(inp)
	disp(sprintf('ASORT> input is empty'));
	return;
end

ti=clock;
winp=whos('inp');
switch	winp.class
	case	'cell'
		if	~iscellstr(inp)
			disp(sprintf('ASORT> cell is not an array of strings'));
			return;
		end
		inp=inp(:);
		[ins,inx]=sort(inp);
	case	'char'
		%		[ins,inx]=sortrows(inp);
		inp=cstr(inp);
	otherwise
		disp(sprintf('ASORT> does not sort input of class <%s>',winp.class));
		return;
end

inp=inp(:);
inp=setinp(inp,tmpl,[tflg wflg]);
[ins,inx]=sort(inp);
if	strcmp(smod,'descend')
	ins=ins(end:-1:1,:);
	inx=inx(end:-1:1);
end
ins=inp(inx);
c=lower(char(ins));
wins=whos('c');
[cr,cc]=size(c);

% - LEXICAL PARSER
%--------------------------------------------------------------------------------
% - extend input on either side for search
c=[' '*ones(cr,2) c ' '*ones(cr,2)];

% - search for valid alphanumeric items in strings
%   numbers/signs
t=(c>='0'&c<='9');
t=t|c=='-';
t=t|c=='+';
[tr,tc]=size(t);
%   decimal points
%   note: valid numbers with dec points must follow these templates
%         nr.nr
%	  sign.nr
%         nr.<SPACE>
%         <SPACE>.nr
ix1=	 t(:,1:end-2) & ...
	~isletter(c(:,1:end-2)) & ...
	c(:,2:end-1)=='.';
t(:,2:end-1)=t(:,2:end-1)|ix1;
ix1=	(t(:,3:end) & ...
	(~isletter(c(:,3:end)) & ...
	~isletter(c(:,1:end-2))) | ...
	(c(:,3:end)=='e' | ...
	c(:,3:end)=='d')) & ...
	c(:,2:end-1)=='.';
t(:,2:end-1)=t(:,2:end-1)|ix1;
%		t(:,3:end)=t(:,3:end)|ix1;
%   signs
t(c=='-')=false;
t(c=='+')=false;
ix1=	 t(:,3:end) & ...
	(c(:,2:end-1)=='-' | ...
	c(:,2:end-1)=='+');
t(:,2:end-1)=t(:,2:end-1)|ix1;
%   exponents
ix1=	 t(:,1:end-2) & ...
	(c(:,2:end-1)=='e' | ...
	c(:,2:end-1)=='d');
t(:,2:end-1)=t(:,2:end-1)|ix1;
%   spec numbers
c=reshape(c.',1,[]);
t=t';
ic=[];
for	j=1:numel(ntmpl)
	ic=[ic,strfind(c,ntmpl{j})];
end
ic=sort(ic);
for	i=1:numel(ic)
	ix=ic(i)+0:ic(i)+4;
	t(ix)=true;
end
t=t';
c=reshape(c.',[tc,tr]).';
t(c==' ')=false;
%--------------------------------------------------------------------------------

% - only allow one number per string
il=~any(t,2);
ib=strfind(reshape(t.',1,[]),[0 1]);
if	~isempty(ib)
	ixe=cell(3,1);
	n=reshape(char(t.*c).',1,[]);
	for	i=1:numel(ctmpl)
		id=strfind(n,ctmpl{i});
		if	~isempty(id)
			[dum,dum,ixu{i},ixe{i}]=dupinx(id,tc);
		end
	end
	in=false(tr,1);
	im=in;
	%   must check for anomalous cases like <'.d'>
	id=sort(...
		[find(n>='0' & n<='9'),...
		strfind(n,'inf'),...
		strfind(n,'nan')]);
	%		[ibu,ibd,ixbu,ixe{i+1}]=dupinx(id,tc);
	[ibu,ibd,ixbu,ixbd]=dupinx(id,tc);
	in(ixbu)=true;
	in(ixbd)=true;
	[ibu,ibd,ixbu,ixbd]=dupinx(ib,tc);
	im(ixbu)=true;
	in=in&im;
	in([ixe{:}])=false;
	il=~any(t,2);
	ia=~(in|il);

	% - read valid strings
	n=t(in,:).*c(in,:);
	n(n==0)=' ';
	n=char(n);
	dn=strread(n.','%n');
	if	numel(dn) ~= numel(find(in))
		%disp(sprintf('ASORT> unexpected fatal error reading input!'));
		if	nargout
			s.c=c;
			s.t=t;
			s.n=n;
			s.d=dn;
			varargout{1}=s;
		end
		return;
	end

	% - sort numbers
	[ds,dx]=sort(dn,1,smod);
	in=find(in);
	anr=ins(in(dx));
	snr=ins(ia);
end
str=ins(il);
to=clock;

% - prepare output
if	nargout < 3 || sflg
	s.magic='ASORT';
	s.ver='30-Mar-2005 11:57:07';
	s.time=datestr(clock);
	s.runtime=etime(to,ti);
	s.input_class=winp.class;
	s.input_msize=winp.size;
	s.input_bytes=winp.bytes;
	s.strng_class=wins.class;
	s.strng_msize=wins.size;
	s.strng_bytes=wins.bytes;
	s.anr=anr;
	s.snr=snr;
	s.str=str;
	if	dflg
		s.c=c;
		s.t=t;
		s.n=n;
		s.d=ds;
	end
	varargout{1}=s;
else
	s={anr,snr,str};
	for	i=1:nargout
		varargout{i}=s{i};
	end
end

if	vflg
	inp=cstr(inp);
	an=[{'--- NUMERICAL'};		anr];
	as=[{'--- ASCII NUMBERS'};	snr];
	at=[{'--- ASCII STRINGS'};	str];
	nn=[{'--- NUMBERS'};		num2cell(ds)];
	ag={' ';' ';' '};
	u=[{'INPUT'};			inp;ag];
	v=[{'ASCII SORT'};		ins;ag];
	w=[{'NUM SORT'};		an;as;at];
	x=[{'NUM READ'};		nn;as;at];
	w=[u,v,w,x];
	disp(w);
end

return;
%--------------------------------------------------------------------------------
function	c=cstr(s)
% - bottleneck waiting for a good <cellstr> replacement
%   it consumes ~75% of <asort>'s processing time!

c=s;
if	ischar(s)
	sr=size(s,1);
	c=cell(sr,1);
	for	i=1:sr
		c{i}=s(i,:);	% no deblanking!
	end
end
return;
%--------------------------------------------------------------------------------
function	[idu,idd,ixu,ixd]=dupinx(ix,nc)
% - check for more than one entry/row in a matrix of column size <nc>
%   unique    indices:	idu / ixu
%   duplicate indices:	idd / ixd

if	isempty(ix)
	idu=[];
	idd=[];
	ixu=[];
	ixd=[];
	return;
end
id=fix(ix/nc)+1;
idi=diff(id)~=0;
ide=[true idi];
idb=[idi true];
idu=idb & ide;
idd=idb==1 & ide==0;
ixu=id(idu);
ixd=id(idd);
return;
%--------------------------------------------------------------------------------
function	inp=setinp(inp,tmpl,flg)
% - remove space(s) and/or templates

if	isempty(inp) || ~any(flg)
	return;
end

for	i=sort(flg)
	switch	i
		case	flg(1)
			if	ischar(tmpl)
				tmpl={tmpl};
			end
			for	i=1:numel(tmpl)
				inp=strrep(inp,tmpl{i},' ');
			end
		case	flg(2)
			inp=strrep(inp,' ','');
	end
end
return;
%--------------------------------------------------------------------------------

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
clear;clc

%% Designing a decimation filter.
% We start by desinging a decimation filter in double precision. The
% decimation filter is a low pass filter with a decimation rate of 64. The
% filter specification is from the AD1877 Anolog Devices Sigma-Delta
% data sheet. You can find the data sheet at
% http://www.analog.com/UploadedFiles/Data_Sheets/AD1877.pdf

% The desing of the filter takes about 10 seconds. As such, for a faster
% start up the design is stored in a file. If you need to redesign the
% filter, simply uncomment the lines below and run the script.

Decimation_Factor=64;
Passband_Ripple=.006; %dB
Stopband_Attenuation=90; %dB
Fs=48e3;
Passband=21.6e3;  %dB
Stopband=26.4e3; %dB

Input_Sampling_Rate = Decimation_Factor*Fs;
f=fdesign.decimator(Decimation_Factor,'lowpass',Passband,Stopband,...
    Passband_Ripple,Stopband_Attenuation,Input_Sampling_Rate);
h=design(f);

save one_stage h
fvtool(h,'Fs',Input_Sampling_Rate);

%% To have a more efficnet implementation, we break down the filter into

hm=design(f,'multistage','Nstages',3);

save multi_stage hm
fvtool(h,'Fs',Input_Sampling_Rate);

%% We then convert the filter into fixed-point.
hf=hm;

hf.stage(1).Arithmetic = 'fixed';
hf.stage(1).CoeffWordLength = 24;
hf.stage(1).InputWordLength = 2;
hf.stage(1).InputFracLength = 0;
specifyall(hf.stage(1));
hf.stage(1).OutputWordLength = 12;
hf.stage(1).OutputFracLength = 10;
hf.stage(1).ProductWordLength = 8;
hf.stage(1).ProductFracLength = 10;
hf.stage(1).AccumWordLength = 13;
hf.stage(1).AccumFracLength = 10;
hf.stage(1).RoundMode = 'nearest';

hf.stage(2).Arithmetic = 'fixed';
hf.stage(2).CoeffWordLength = 24;
hf.stage(2).InputWordLength = 13;
hf.stage(2).InputFracLength = 10;
specifyall(hf.stage(2));
hf.stage(2).OutputWordLength = 13;
hf.stage(2).OutputFracLength = 10;
hf.stage(2).ProductWordLength = 10;
hf.stage(2).ProductFracLength = 10;
hf.stage(2).AccumWordLength = 13;
hf.stage(2).AccumFracLength = 10;
hf.stage(2).RoundMode = 'nearest';

hf.stage(3).Arithmetic = 'fixed';
hf.stage(3).CoeffWordLength = 24;
hf.stage(3).InputWordLength = 13;
hf.stage(3).InputFracLength = 10;
specifyall(hf.stage(3));
hf.stage(3).ProductWordLength = 11;
hf.stage(3).ProductFracLength = 12;
hf.stage(3).AccumWordLength = 14;
hf.stage(3).AccumFracLength = 12;
hf.stage(3).OutputWordLength = 16;
hf.stage(3).OutputFracLength = 15;
hf.stage(3).RoundMode = 'nearest';
hf.stage(3).OverflowMode = 'saturate';

save multi_stage_fixed hf

fvtool(hf,'Fs',Input_Sampling_Rate);


text:
function [ max_adif, max_adif_i, max_adif_j, max_rdif, max_rdif_i, max_rdif_j ] = ...
  p00_jac_check ( problem, option, nvar, x )

%*****************************************************************************80
%
%% P00_JAC_CHECK compares the jacobian with a finite difference estimate.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the jacobian.
%
%    Output, real MAX_ADIF, the maximum absolute difference.
%
%    Output, integer MAX_ADIF_I, MAX_ADIF_J, the indices where
%    the maximmum absolute difference was found.
%
%    Output, real MAX_RDIF, the maximum relative difference.
%
%    Output, integer MAX_RDIF_I, MAX_RDIF_J, the indices where
%    the maximmum relative difference was found.
%
  rel = 0.0001;
%
%  Compute the jacobian.
%
  jac = p00_jac ( problem, option, nvar, x );
  jac_norm = max ( max ( abs ( jac ) ) );
%
%  Estimate the jacobian via finite differences.
%
  jac_dif = p00_jac_dif ( problem, option, nvar, x );
%
%  Compare the jacobians.
%
  max_rdif = 0.0;
  max_rdif_i = 0;
  max_rdif_j = 0;

  max_adif = 0.0;
  max_adif_i = 0;
  max_adif_j = 0;

  for i = 1 : nvar - 1
    for j = 1 : nvar

      dif = abs ( jac(i,j) - jac_dif(i,j) );

      if ( max_adif < dif )
        max_adif = dif;
        max_adif_i = i;
        max_adif_j = j;
      end

      if ( rel < abs ( jac(i,j) ) )
        if ( max_rdif * abs ( jac(i,j) ) < dif )
          max_rdif = dif / abs ( jac(i,j) );
          max_rdif_i = i;
          max_rdif_j = j;
        end
      end

    end
  end

  return
end


text:
classdef prtKernelRbfFixed < prtKernel

    properties (SetAccess = private)
        name = 'RBF Kernel'; % RBF Kernel
        nameAbbreviation = 'RBF'; % RBF
    end

    properties
        sigma = 1; % The inverse kernel width
        x0 = 0;
    end

    methods (Hidden = true)

        function nDimensions = nDimensions(self)
            nDimensions = 1;
        end

    end
    methods (Access = protected, Hidden = true)
        function self = trainAction(self,ds)
            self.internalDataSet = ds;
            self.isTrained = true;
        end

        function dsOut = runAction(self,ds)
            if ~self.isTrained
                error('prtKernelRbf:run','Attempt to run an untrained kernel; use kernel.train(ds) to train');
            end
            if self.internalDataSet.nObservations == 0
                dsOut = prtDataSetClass;
            else
                gram = prtKernelRbf.kernelFn(ds.getObservations,self.x0,self.sigma);
                dsOut = ds.setObservations(gram);
            end
        end
    end

    methods
        function self = prtKernelRbfFixed(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end

        function self = set.sigma(self,value)
            if ~prtUtilIsPositiveScalar(value)
                assert(isnumeric(value) && all(value>0) && isvector(value),'sigma must be a positive numeric vector')
                if isempty(self.internalDataSet) || self.internalDataSet.nObservations==0
                    error('prtKernelRbf:set','Value of sigma must be a positive scalar');
                else
                    assert(self.internalDataSet.nObservations==numel(value),'When setting sigma to be an array of values the internalDataSet must be set and the number of observations and the length of sigma must match');
                end
            end
            self.sigma = value(:);
        end
    end

    methods(Hidden = true)
        function varargout = plot(obj)
            x = obj.internalDataSet.getObservations;

            if size(x,2) <= 3
                if size(x,2) == 1 && obj.internalDataSet.isLabeled
                    xy = cat(2,x,obj.internalDataSet.getTargets);
                    h = prtPlotUtilScatter(xy, {}, obj.plotOptions.symbol, obj.plotOptions.markerFaceColor, obj.plotOptions.color, obj.plotOptions.symbolLineWidth, obj.plotOptions.symbolSize);
                else
                    h = prtPlotUtilScatter(x, {}, obj.plotOptions.symbol, obj.plotOptions.markerFaceColor, obj.plotOptions.color, obj.plotOptions.symbolLineWidth, obj.plotOptions.symbolSize);
                end
            else
                h = nan;
            end

            varargout = {};
            if nargout
                varargout = {h};
            end
        end
    end

    methods (Static, Hidden = true)
        function gram = kernelFn(x,y,sigma)
            [n1, d] = size(x);
            [n2, nin] = size(y);
            %             if d ~= nin
            %                 error('size(x,2) must equal size(y,2)');
            %             end
            %             keyboard
            %dist2 = prtDistanceLNorm(x,y,2);
            dist2 = repmat(sum((x.^2), 2), [1 n2]) + repmat(sum((y.^2),2), [1 n1]).' - 2*x*(y.');

            if numel(sigma) == 1
                gram = exp(-dist2/(sigma.^2));
            else
                gram = exp(-bsxfun(@rdivide,dist2,(sigma.^2)'));
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [varargout] = applyTransM(varargin)
%"applyTransM"
%   Apply transformation matrix transM to points defined in xyz by xV, yV,
%   zV.  In order to run as fast as possible, first computes the rotation
%   component and then adds the translation portion of transM last.
%
%JRA 1/18/04
%
%Usage:
%   function [xT, yT, zT] = applyTransM(transM, xV, yV, zV);
%   function [pointsT]    = applyTransM(transM, pointsM);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Parse input arguments.
if nargin == 2
    usingPointMatrix = 1;
    transM  = varargin{1};
    pointsM = varargin{2}';
elseif nargin == 4
    usingPointMatrix = 0;
    transM  = varargin{1};
    xV      = varargin{2};
    yV      = varargin{3};
    zV      = varargin{4};

    if length(xV) ~= length(yV) | length(xV) ~= length(zV)
        error('xV, yV, and zV must be vectors of the same length.');
    end

    pointsM = [reshape(xV,[],1), reshape(yV,[],1) reshape(zV,[],1)]';

else
    error('Invalid number of input arguments to applyTransM.');
end

%If blank transformation matrix, return originals.
if isempty(transM)
    if usingPointMatrix
        varargout{1} = pointsM;
    else
        varargout{1} = xV;
        varargout{2} = yV;
        varargout{3} = zV;
    end
    return;
end

nPts = size(pointsM, 2);
%If no points passed in, return empty.
if nPts == 0
    if usingPointMatrix
        varargout{1} = [];
    else
        varargout{1} = [];
        varargout{2} = [];
        varargout{3} = [];
    end
    return;
end

%Split the rotation and translation portions of transM.  This is done for
%speed, to avoid allocating a 4th column of ones to pointsM.

%Rotation.
pointsM = transM(1:3,1:3) * pointsM;

%Translation.
pointsM(1,:) = pointsM(1,:) + transM(1,4);
pointsM(2,:) = pointsM(2,:) + transM(2,4);
pointsM(3,:) = pointsM(3,:) + transM(3,4);

if usingPointMatrix
    varargout{1} = pointsM';
else
    varargout{1} = pointsM(1,:);
    varargout{2} = pointsM(2,:);
    varargout{3} = pointsM(3,:);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function pts2annotations(directory)
% pts2annotations(directory)
%
% Converts annotations from pts format to our native .mat format.
%
%                         PARAMETERS
%
% directory Path to directory where the .pts files are stored,
%           along with the associated images.
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

landmark_files = dir(sprintf('%s/*.pts', directory));

app = imread(sprintf('%s/%s', directory, landmark_files(1).name(1:end-4)));

for i=1:numel(landmark_files)
	fid = fopen(sprintf('%s/%s', directory, landmark_files(i).name), 'r');

	% First line must be "version: 1"
	data = fgetl(fid);
	assert(strcmp(data, 'version: 1'));

	% Second line must be "n_points: <n>"
	num_points = fscanf(fid, 'n_points: %d\n');
	data = fgetl(fid);

	% Third line is "{"
	assert(strcmp(data, '{'));

	% Read annotation data
	for j=1:num_points
		P = fscanf(fid, '%f', 2);
		shape(j,1) = P(1);
		shape(j,2) = P(2);
	end

	% Match end of line
	data = fgetl(fid);
  assert(strcmp(data, ''));

  % Last line is "}"
  data = fgetl(fid);
	assert(strcmp(data, '}'));

	fclose(fid);

	% Convert to cartesian coordinates
	annotations(:,1) = shape(:,1);
	annotations(:,2) = repmat(size(app, 1), [num_points 1]) - shape(:,2);

	save(sprintf('%s/%s.mat', directory, landmark_files(i).name(1:end-4)), 'annotations');
end


text:
function [y,bw] = erbspace(fmin,fmax,n)
%ERBSPACE  Equidistantly spaced points on erbscale
%   Usage: y=erbspace(fmin,fmax,n);
%
%   This is a wrapper around |audspace| that selects the erb-scale. Please
%   see the help on |audspace| for more information.
%
%   See also: audspace, freqtoaud

%   AUTHOR : Peter L. Søndergaard

[y,bw] = audspace(fmin,fmax,n,'erb');

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [minVal, minPos] = min(f)
%MIN   Global minimum of a CLASSICFUN on [a,b].
%   MINVAL = MIN(F) returns the global minimum of the CLASSICFUN F on [a,b].  If F is
%   an array-valued CLASSICFUN, MINVAL is a row vector whose Kth entry is the global
%   minimum of the Kth column of F.
%
%   [MINVAL, MINPOS] = MIN(F) returns also a value such that MINVAL = F(MINPOS).
%
%   If F is complex-valued then absolute values are taken to determine maxima
%   but the output value corresponds to that of the original function. That
%   is, MINVAL = FEVAL(F, MINPOS) where [~, MINPOS] = MIN(ABS(F)).
%
% See also MAX, MINANDMAX.

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

% To avoid code duplication, we simply call MINANDMAX:
[minVal, minPos] = minandmax(f);
minVal = minVal(1,:);
minPos = minPos(1,:);

end

author: chebfun

repo: chebfun


text:
function value = r8_abs ( x )

%*****************************************************************************80
%
%% R8_ABS returns the absolute value of an R8.
%
%  Discussion:
%
%    An R8 is a double precision real value.
%
%    MATLAB supplies the ABS function, which should be used instead
%    of this function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose absolute value is desired.
%
%    Output, real R8_ABS, the absolute value of X.
%
  if ( 0.0 <= x )
    value = + x;
  else
    value = - x;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [lpp] = tapas_sem_prosa_lpp(theta, ptheta)
%% Log prior probability of the parameters of the prosa model.
%
% Input
%   pv -- Parameters in numerical form.
%   theta -- Parameters of the model
%   ptheta -- Priors of parameters
%
% Output
%   lpp -- Log prior probability
%

% aponteeduardo@gmail.com
% copyright (C) 2015
%

nconst = ptheta.pconst;

lpp = zeros(1, numel(theta));
for i = 1:numel(theta)
    if isfield(ptheta, 'dkjm')
        err = sum((ptheta.mu - theta{i}) .* ptheta.dkjm .* ...
        (ptheta.mu - theta{i}));
    else
        err = (ptheta.mu - theta{i})' * ptheta.kjm * (ptheta.mu - theta{i});
    end

    alpha = ptheta.mu(ptheta.bdist);
    beta = ptheta.pm(ptheta.bdist);
    lcosh = log(cosh(theta{i}(ptheta.bdist)/2));
    lpp(i) = nconst - 0.5 * err + ...
       ...sum(betaln(alpha, beta) + ...
        sum((alpha - 1) .* (theta{i}(ptheta.bdist)/2 - log(2) - ...
            lcosh) + ...
        (beta - 1) .* (-theta{i}(ptheta.bdist)/2 - log(2) - ...
            lcosh)) - ...
        2 * sum(log(2) + lcosh);
end

end

author: translationalneuromodeling

repo: tapas


text:
function c = tapas_align_priors_fields(c)
% Aligns parameter fields with the explicit prior definitions with the
% content of the vectors c.priormus and c.priorsas (vice-versa of function
% 'tapas_align_priors.m').
%
%   Example:
%   >> c_prc = tapas_ehgf_binary_config;
%   >> c_prc.priormus(13) = -2;
%   >> c_prc = tapas_align_priors_fields(c_prc)
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Get fieldnames. If a name ends on 'mu', that field defines a prior mean.
% If it ends on 'sa', it defines a prior variance.
names = fieldnames(c);
pm = 0;
ps = 0;

% Loop over fields and overwrite fiels whose name ends on 'mu or 'sa'
for i = 1:length(names)
    if regexp(names{i}, 'mu$')
        c.(names{i}) = c.priormus(pm+1:pm+length(c.(names{i})));
        pm = pm+length(c.(names{i}));
    elseif regexp(names{i}, 'sa$')
        c.(names{i}) = c.priorsas(ps+1:ps+length(c.(names{i})));
        ps = ps+length(c.(names{i}));
    end
end

end


text:
function value = p13_autonomous ( )

%*****************************************************************************80
%
%% P13_AUTONOMOUS reports whether problem P13 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end


text:
function gauss = get_gaussian_kernel(Nsigmas, sigma_in_pix)
% generate gaussian kernel with given sigmas
    ndims = numel(sigma_in_pix);
    width = 2 * floor(Nsigmas .* sigma_in_pix) + 1;
    hw = (width - 1)/2;

    sigma_in_pix(sigma_in_pix == 0) = 1;
    if ndims == 3
        -hw(1):hw(1)
        hw(2)
        hw(3)
        [n1, n2, n3] = ndgrid(-hw(1):hw(1), -hw(2):hw(2), -hw(3):hw(3));
        size(n1)
        weight = exp(-n1.^2/(2*sigma_in_pix(1)^2) - n2.^2/(2*sigma_in_pix(2)^2) - n3.^2/(2*sigma_in_pix(3)^2) );
    elseif ndims == 2
        [n1, n2] = ndgrid(-hw(1):hw(1), -hw(2):hw(2));
        weight = exp(-n1.^2/(2*sigma_in_pix(1)^2) - n2.^2/(2*sigma_in_pix(2)^2));
    end

    gauss = weight / sum(weight(:));
end

repo: pTVreg


text:
function [cropsResized, preparationTime, resizeTime] = cropImagePatches(curImage, boundingBoxes, cropPad, outputSize, maxGpuImages, cropMode, jitterStd)
%cropImagePatches crops patches from an image and resizes them to the standard size
%
% cropsResized = cropImagePatches(curImage, boundingBoxes, cropPad, outputSize);
% [cropsResized, cropTime, resizeTime] = cropImagePatches(curImage, boundingBoxes, cropPad, outputSize, maxGpuImages, cropMode, jitterStd);
%
% Input:
%   cropImage - image to crop from (has to be single)
%   boundingBoxes - bounding boxes in format [y1 x1 y2 x2]. x is col, y is row
%   cropPad - padding of the bounding boxes (padding is measure after crop and resize), format: left, top, right, bottom
%   outputSize - target size of the crops
%   maxGpuImages - maximal number of patches to crop on a GPU at the same time (default : 256)
%   cropMode - 'square' or 'warp' (mimicing the R-CNN cropping code) (default: 'warp')
%   jitterStd - ammount of jittering to do: std of a gaussian distribution of a patch-border shift (default: 0)
%
% Output:
%   cropsResized - the cropped patches
%   preparationTime - time on the computations of the crop parameters
%   resizeTime - time spent on cropping and resizing
%
% This function depends on cropRectanglesMex

tStart = tic;

if ~exist( 'maxGpuImages', 'var') || isempty(maxGpuImages)
    maxGpuImages = 256;
end
if ~exist( 'cropMode', 'var') || isempty(cropMode)
    cropMode = 'warp';
end
if ~exist('jitterStd', 'var') || isempty(jitterStd)
    jitterStd = 0;
end
useSquare = false;
if isequal( cropMode, 'square')
    useSquare = true;
end

numChannels = size(curImage, 3);
numBoxes = size( boundingBoxes, 1 );
if numBoxes <= maxGpuImages
    resultOnGpuGlag = true;
else
    resultOnGpuGlag = false;
end

cropBoxes = nan(numBoxes, 4);
for iBox = 1 : numBoxes
    %% get crop position: take padding into account
    leftBorder = boundingBoxes(iBox, 2);
    rightBorder = boundingBoxes(iBox, 4);

    topBorder = boundingBoxes(iBox, 1);
    bottomBorder = boundingBoxes(iBox, 3);

    % add the square mode
    if useSquare
        halfWidth = (rightBorder - leftBorder + 1 ) / 2;
        halfHeight = (bottomBorder - topBorder + 1 ) / 2;

        centerX = leftBorder + halfWidth;
        centerY = topBorder + halfHeight;

        if halfHeight > halfWidth
            halfWidth = halfHeight;
        else
            halfHeight = halfWidth;
        end

        topBorder = centerY - halfHeight;
        bottomBorder = centerY + halfHeight;
        leftBorder = centerX - halfWidth;
        rightBorder = centerX + halfWidth;
    end

    % compute the transformation
    posOld = [  topBorder, leftBorder, 1; ... % top left corner
        topBorder, rightBorder, 1; ... % top right corner
        bottomBorder, leftBorder, 1; ... % bottom left corner
        bottomBorder, rightBorder, 1; ... % bottom right corner
        ]';

    % crop corners after padding
    posNew = [  1 + cropPad(2), 1 + cropPad(1), 1; ... % top left corner
        1 + cropPad(2), outputSize(2) - cropPad(3), 1; ... % top right corner
        outputSize(1) - cropPad(4), 1 + cropPad(1), 1; ... % bottom left corner
        outputSize(1) - cropPad(4), outputSize(2)- cropPad(3), 1; ... % bottom right corner
        ]';

    % solve linear system
    M = posOld / posNew;

    % actual crop corners
    posNew = [  1, 1, 1; ... % top left corner
        1, outputSize(2), 1; ... % top right corner
        outputSize(1), 1, 1; ... % bottom left corner
        outputSize(1), outputSize(2), 1; ... % bottom right corner
        ]';

    cropPos = (M * posNew)';

    leftBorder = (cropPos(1, 2) + cropPos(3, 2)) / 2;
    rightBorder = (cropPos(2, 2) + cropPos(4, 2)) / 2;

    topBorder = (cropPos(1, 1) + cropPos(2, 1)) / 2;
    bottomBorder = (cropPos(3, 1) + cropPos(4, 1)) / 2;

    cropBoxes(iBox, 1) = topBorder;
    cropBoxes(iBox, 2) = leftBorder;
    cropBoxes(iBox, 3) = bottomBorder;
    cropBoxes(iBox, 4) = rightBorder;

end

%% apply jittering
borderNoise = randn( size( cropBoxes) ) * jitterStd;
cropBoxes = cropBoxes + borderNoise;
% check if the bouding boxes are still valid. If not (might happen if the size of the box is too small) than remove jittering
badMask = cropBoxes(:, 1) > cropBoxes(iBox, 3) | cropBoxes(:, 2) > cropBoxes(iBox, 4);
cropBoxes( badMask, : ) = cropBoxes( badMask, : ) - borderNoise( badMask, : );

preparationTime = toc(tStart);

%% crop the image
if ~isequal( class( curImage ), 'single')
    curImage = single(curImage);
end
if resultOnGpuGlag
    cropsResized = cropRectanglesMex( curImage, cropBoxes, outputSize );
    cropsResized = single(cropsResized);
else
    cropsResized = zeros(outputSize(1), outputSize(2), numChannels, numBoxes, 'single');
    for iBatchStart = 1 : maxGpuImages : numBoxes
        curIds = iBatchStart : 1 : min( numBoxes, iBatchStart + maxGpuImages - 1);
        curBoxes = cropBoxes( curIds, : );
        curCrops = cropRectanglesMex( curImage, curBoxes, outputSize );
        cropsResized(:,:,:,curIds) = gather(curCrops);
    end
    cropsResized = single(cropsResized);
end

resizeTime = toc(tStart) - preparationTime;

end

author: aosokin


text:
%                               applyOpts.m
%
% This function checks if an element is in a given list. If it is not,
% an error will be raised. It is used in any place where a value is
% expected to be chosen from a list of options
% (e.g. It is used to check opts.indexChoice in the function
% validateOptions in solveCoordinateDescent.m).
%
% Inputs:
%         fieldName(string)       :  The name of the field to be checked.
%         element(string)         :  The name of the element to be checked.
%         list(cell array)        :  The list that consists of all the
%                                    valid names.
%         customizedMsg(string)   :  A customized message put at the
%                                    beginning of the error message.
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

function checkIfInList(fieldName,element,list,customizedMsg)
    if ~exist('customizedMsg','var')
        customizedMsg = '';
    end
    if ~any(strcmp(list,lower(element)))
        % Use num2str in case element is given as a number
        error('%s %s is invalid for field %s\n',customizedMsg,num2str(element),fieldName);
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
function kernel = create_kernel(kernel_type, pars, nMax, lb, ub, bound_pars)
%% create convolution kernel
%% inputs:
%   kernel_type: string, convolution kernel type. now support {'exp',
%       'exp2', 'vector'}
%   pars: parameters for the selected kernel type
%   nMax: length of the kernel
%   lb:     lower bound for each parameter
%   ub:     upper bound for each parameter
%   bound_pars: logical variable, bound the parameters or not {1, 0}
%% outputs
%   kernel: struct variable

%% Author: Pengcheng Zhou, Carnegie Mellon University, 2016

%% kernel size
if ~exist('nMax', 'var') || isempty(nMax)
    nMax = 50;
end
kernel.nMax = nMax;

%% initialize kernel
if ~exist('kernel_type', 'var') || isempty(kernel_type)
    kernel_type = 'exp2';
end
if strcmpi(kernel_type, 'exp')
    % single exponential function: ~ exp(-t/tau)
    kernel.type = 'exp';
    % parameter
    if ~exist('pars', 'var') || isempty(pars)
        kernel.pars = [5, .1];
    else
        kernel = kernel.pars;
    end
    % function handle
    kernel.fhandle = @(pars, t) exp(-t/pars) * (1-exp(-1/pars)) ...
        / (1-exp(-nMax/pars));
elseif strcmpi(kernel_type, 'vector')
    % single vector
    kernel.type = 'vector';
    % parameter
    if ~exist('pars', 'var') || isempty(pars)
        kernel.pars = exp(-(1:nMax)/10);
    else
        kernel.pars = pars;
    end
    % function handle
    kernel.fhandle = @(pars, t) pars/sum(pars);
else
    % differencing of two exponential function:
    % ~  exp(-t/tau_d)-exp(-t/tau_r)
    kernel.type = 'exp2';

    % parameters
    if ~exist('pars', 'var') || isempty(pars)
        kernel.pars = [10, 1];
    else
        kernel.pars = pars;
    end
    % function handle
    kernel.fhandle = @(pars, t) (exp(-t/pars(1)) - exp(-t/pars(2)))  ...
        /( (1-exp(-nMax/pars(1)))/(1-exp(-1/pars(1))) ...
        - (1-exp(-nMax/pars(2)))/(1-exp(-1/pars(2))));
end

% lower and upper bounds for parameters
if ~exist('lb', 'var') || isempty(lb)
    kernel.lb = 0.5*kernel.pars;
else
    kernel.lb = lb;
end
if ~exist('ub', 'var') || isempty(ub)
    kernel.ub = 2*kernel.pars;
else
    kernel.ub = ub;
end

% bound the parameters of not
if ~exist('bound_pars', 'var')||isempty(bound_pars)
    kernel.bound_pars = false;
else
    kernel.bound_pars = bound_pars;
end

author: zhoupc


text:
function trans = localToGlobal3d(varargin)
%LOCALTOGLOBAL3D Transformation matrix from local to global coordinate system.
%
%   TRANS = localToGlobal3d(CENTER, THETA, PHI, PSI)
%   Compute the transformation matrix from a local (or modelling)
%   coordinate system to the global (or world) coordinate system.
%   This is a low-level function, used by several drawing functions.
%
%   The transform is defined by:
%   - CENTER: the position of the local origin into the world coordinate
%       system
%   - THETA: colatitude, defined as the angle with the Oz axis (between 0
%       and 180 degrees), positive in the direction of the of Oy axis.
%   - PHI: azimut, defined as the angle of the normal with the Ox axis,
%       between 0 and 360 degrees
%   - PSI: intrinsic rotation, corresponding to the rotation of the object
%       around the direction vector, between 0 and 360 degrees
%
%   The resulting transform is obtained by applying (in that order):
%   - Rotation by PSI   around the Z-axis
%   - Rotation by THETA around the Y-axis
%   - Rotation by PHI   around the Z-axis
%   - Translation by vector CENTER
%   This corresponds to Euler ZYZ rotation, using angles PHI, THETA and
%   PSI.
%
%   The 'eulerAnglesToRotation3d' function may better suit your needs as
%   it is more 'natural'.
%
%   Example
%   localToGlobal3d
%
%   See also
%   transforms3d, eulerAnglesToRotation3d

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2009-06-19, using Matlab 7.7.0.471 (R2008b)
% Copyright 2009-2022 INRA - Cepia Software Platform

% extract the components of the transform
if nargin == 1
    % all components are bundled in  the first argument
    var     = varargin{1};
    center  = var(1:3);
    theta   = var(4);
    phi     = var(5);
    psi     = 0;
    if length(var) > 5
        psi = var(6);
    end

elseif nargin == 4
    % arguments = center, then the 3 angles
    center  = varargin{1};
    theta   = varargin{2};
    phi     = varargin{3};
    psi     = varargin{4};

elseif nargin > 4
    % center is given in 3 arguments, then 3 angles
    center  = [varargin{1} varargin{2} varargin{3}];
    theta   = varargin{4};
    phi     = varargin{5};
    psi     = 0;
    if nargin > 5
        psi = varargin{6};
    end
end

% conversion from degrees to radians
k = pi / 180;

% rotation around normal vector axis
rot1    = createRotationOz(psi * k);

% colatitude
rot2    = createRotationOy(theta * k);

% longitude
rot3    = createRotationOz(phi * k);

% shift center
tr      = createTranslation3d(center);

% create final transform by concatenating transforms
trans   = tr * rot3 * rot2 * rot1;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function view = rmInitView(viewPointer,roiFileName)
% rmInitView - initiate view struct from key pointers
%
% view = rmInitView(view,roiFileName);
%
% These key pointers now are:
%         pointer(1): viewType
%         pointer(2): dataType
%         roiFileName: if defined roi will be loaded
%
% SOD: wrote it.

if ~exist('viewPointer','var') || isempty(viewPointer), error('Need viewPointers'); end;
if ~exist('roiFileName','var'), roiFileName = []; end;

% default vista stuff
mrGlobals;
loadSession;

% set viewType
if viewPointer(1)==1
    view=initHiddenGray;
else
    view=initHiddenInplane;
end

% set dataType
view = viewSet(view,'curdatatype',viewPointer(2));
fprintf(1,'[%s]:DataType: %s\n',mfilename, ...
    dataTYPES(view.curDataType).name);

% if roiFileName load roi
if ~isempty(roiFileName)
    if ~strcmp(roiFileName,'0')
        view = loadROI(view,roiFileName,[],[],[],1);
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function tcpWithBoost = seanWalshTCP(paramS,doseBinsC,volHistC)
% Sean Walsh TCP model for prostate
% http://dx.doi.org/10.1118/1.4939260
%
% APA, 07/20/2016
% AI, 09/12/16 Modified for use with outcomeModelsGUI

% Get parameters
alpha = 0.25; %paramS.alpha.val;
beta = 0.1008;% paramS.beta.val;
sigmaAlpha = alpha*11.3/100;
sigmaBeta = beta*12.9/100;

n = paramS.numFractions.val;                      % number of fractions

DILparS = paramS.structures.DIL;
dilVolume = DILparS.dilVolume.val;                % percent
pDIL = DILparS.pDIL.val;                          % DIL density
dDIL = DILparS.dDIL.val;                          % dose per fraction to the DIL

PTVparS = paramS.structures.PTV_1;
hypoxicFraction = PTVparS.hypoxicFraction.val;    % percent
OER = PTVparS.OER.val;                            % Oxygen Enhancement Ratio
pCTV = PTVparS.pCTV.val;                          % Prostate density
ctvVolume = PTVparS.ctvVolume.val;                % 36 cm^3 for intermediate risk
                                                  % 72 cm^3 for high risk pts
numSimulations = PTVparS.numSimulations.val;      % number of (alpha,beta) simulations

%Compute mean dose to prostate
prost = calc_meanDose(doseBinsC{2},volHistC{2},1);
%Dose per fraction to the prostate
dProst = prost/n;

% Compute TCP
%Seed the random generator
rng(now,'twister');

%Generate normally distributed, positive aplha and beta
alphaV = alpha + randn(numSimulations*2, 1) * sigmaAlpha;
betaV = beta + randn(numSimulations*2, 1) * sigmaBeta;
alphaV = alphaV(alphaV > 0);
betaV = betaV(betaV > 0);
alphaV = alphaV(1:numSimulations);
betaV = betaV(1:numSimulations);

%Hypoxic cells alpha,beta
alphaPO2 = alphaV / OER;
betaPO2 = betaV / OER^2;

% Total initial clonogen number
N0 = ((100-dilVolume)*pCTV + dilVolume*pDIL) / 100 * ctvVolume;

%Initial number of clonogens in prostate
Nprost = (100-dilVolume)*pCTV / 100 * ctvVolume;

%Initial number of clonogens in DIL
Ndil = dilVolume*pDIL / 100 * ctvVolume;

%Surviving fraction for Prostate
SProstateV = exp(-alphaV*n*dProst -betaV*n*dProst^2);

%Surviving fraction for hypoxic Prostate
ShypoxicProstateV = exp(-alphaPO2*n*dProst -betaPO2*n*dProst^2);

%Surviving fraction for DIL
SdilV = exp(-alphaV*n*dDIL -betaV*n*dDIL^2);

%Surviving fraction for hypoxic DIL
ShypoxicDilV = exp(-alphaPO2*n*dDIL -betaPO2*n*dDIL^2);

%Total surviving clonogens
NumSurvivingV = Nprost*(1-hypoxicFraction)*SProstateV + ...
    Nprost*hypoxicFraction*ShypoxicProstateV + ...
    Ndil*(1-hypoxicFraction)*SdilV +  ...
    Ndil*hypoxicFraction*ShypoxicDilV;

% TCP
TCPv = exp(-NumSurvivingV);

%Record the TCP for this DIL dose
tcpWithBoost = mean(TCPv);

end


text:
function value = p55_param ( action, name, value )

%*****************************************************************************80
%
%% P55_PARAM sets or gets real scalar parameters for problem 55.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'GET' to get a parameter.
%    'RANDOMIZE' to randomize a parameter,
%    'SET' to set a parameter.
%
%    Input, string NAME, the name of the variable.
%    'C' is the coefficient.
%    'X0' is the base point.
%
%    Input/output, real VALUE.
%    * If ACTION = 'SET', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'GET', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent c;
  persistent x0;

  if ( s_eqi ( action, 'get' ) )

    if ( s_eqi ( name, 'c' ) )
      if ( isempty ( c ) )
        c = 3.0;
      end
      value = c;
    elseif ( s_eqi ( name, 'x0' ) )
      if ( isempty ( x0 ) )
        x0 = 0.75;
      end
      value = x0;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P55_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name.\n' );
      error ( 'P55_PARAM - Fatal error!' );
    end

  elseif ( s_eqi ( action, 'random' ) )

    if ( s_eqi ( name, 'c' ) )
      c = rand ( );
    elseif ( s_eqi ( name, 'x0' ) )
      x0 = rand ( );
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P55_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name.\n' );
      error ( 'P554_PARAM - Fatal error!' );
    end

  elseif ( s_eqi ( action, 'set' ) )

    if ( s_eqi ( name, 'c' ) )
      c = value;
    elseif ( s_eqi ( name, 'x0' ) )
      x0 = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P55_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name.\n' );
      error ( 'P55_PARAM - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P55_PARAM - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action.\n' );
    error ( 'P55_PARAM - Fatal error!' );

  end

  return
end


text:
function [Fp, Fi] = triangle_triangle_adjacency(F)
  % TRIANGLE_TRIANGLE_ADJACENCY Build a face adjacency data structure for a
  % **manifold** triangle mesh. From each face we can find where on its
  % neighboring faces it's incident.
  %
  % [Fp, Fi] = triangle_triangle_adjacency(F)
  %
  % Input:
  %   F  list of face indices, #faces by 3
  % Output:
  %   Fp  #faces by 3, where Fp(i,j) tells the index of the neighboring triangle
  %     to the jth edge of the ith triangle in F. -1 if the jth edge of the ith
  %     triangle in F is a border edge. (jth edge refers to the edge opposite
  %     the jth vertex: so triangle in a triangle (a,b,c), the 1st edge is
  %     b-->c, the 2nd is c-->b and the 3rd is a-->b
  %   Fi  #faces by 3, where Fi(i,j) tells the position on the neighboring
  %     triangle to the jth edge of the ith triangle in F. -1 if the jth edge
  %     if the ith triangle in F is a border edge. Uses the same indexing of
  %     positions of edges on a triangle as above.
  %
  % For example:
  %
  % F = [ 1 3 2;
  %       2 3 4];
  % [Fp, Fi] = triangle_triangle_adjacency(F);
  % Fp =
  %     2    -1    -1
  %    -1    -1     1
  % Fi =
  %     3    -1    -1
  %    -1    -1     1
  %

  if size(F,2) == 3

    % get list of edges (1st edges then 2nd edges then 3rd edges)
    E = [F(:,2) F(:,3); F(:,3) F(:,1); F(:,1) F(:,2)];

    % sparse adjacency matrix where (i,j) = k, k is the index of i-->j in edge
    % list if i-->j AND j-->i exist in edge list. This way is slightly faster
    % than building a proper adjacency list first.
    Ei = sparse(E(:,1),E(:,2),1:size(E,1));
    % adjacency list of edges as if edges represent undirected graph
    unadj = Ei>0;
    % "non-adjacency" list, (i,j) > 0 iff i-->j exists but j-->i does not exist
    nonadj = unadj-(unadj');
    adj = ((unadj + nonadj)==1).*Ei;

    % need to get mapping from sparse ordering to original order
    [ii jj si] = find(adj);
    % this is slightly faster than sorting the above by ii, which is the same
    [ii jj v] = find(adj');
    % build map from edges to their corresponding faces
    E2F = repmat(1:size(F,1),1,3)';
    % initialize adjacency map to -1
    Fp = -ones(size(E,1),1);
    Fp(si) = E2F(v);
    Fp = reshape(Fp,size(F));
    % build map from edges to edge positions
    I = reshape(repmat((1:3),size(F,1),1),1,3*size(F,1))';
    % use corresponding edges to find positions of correponding faces
    Fi = -ones(size(E,1),1);
    Fi(si) = I(v);
    Fi = reshape(Fi,size(F));
  else
    error('Not supported yet...');
  end

end

author: alecjacobson

repo: gptoolbox


text:
function [X,E,obj,err,iter] = lrtcR_tnn(M,omega,lambda,opts)

% Solve the Noisy Low-Rank Tensor Completion (LRTC) problem by ADMM
%
% min_{X,E} ||X||_*+lambda*loss(E), s.t. P_Omega(X) + E = M.
% loss(E) = ||E||_1 or 0.5*||E||_F^2
%
% ---------------------------------------------
% Input:
%       M       -    d1*d2*d3 tensor
%       omega   -    index of the observed entries
%       lambda  -    >=0, parameter
%       opts    -    Structure value in Matlab. The fields are
%           opts.loss       -   'l1' (default): loss(E) = ||E||_1
%                               'l2': loss(E) = 0.5*||E||_F^2
%           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:
%       X       -    d1*d2*d3 tensor
%       E       -    d1*d2*d3 tensor
%       obj     -    objective function value
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 27/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, '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

dim = size(M);
X = zeros(dim);
Z = X;
E = X;
Y1 = X;
Y2 = X;
omegac = setdiff(1:prod(dim),omega);

iter = 0;
for iter = 1 : max_iter
    Xk = X;
    Zk = Z;
    Ek = E;
    % first super block {X,E}
    [X,tnnX] = prox_tnn(Z-Y2/mu,1/mu);
    temp = M-Y1/mu;
    temp(omega) = temp(omega)-Z(omega);
    if strcmp(loss,'l1')
        E = prox_l1(temp,lambda/mu);
    elseif strcmp(loss,'l2')
        E = temp*(mu/(lambda+mu));
    else
        error('not supported loss function');
    end

    % second super block {Z}
    Z(omega) = (-E(omega)+M(omega)-(Y1(omega)-Y2(omega))/mu+X(omega))/2;
    Z(omegac) = X(omegac)+Y2(omegac)/mu;

    dY1 = E-M;
    dY1(omega) = dY1(omega)+Z(omega);
    dY2 = X-Z;
    chgX = max(abs(Xk(:)-X(:)));
    chgE = max(abs(Ek(:)-E(:)));
    chgZ = max(abs(Zk(:)-Z(:)));
    chg = max([chgX chgE chgZ max(abs(dY1(:))) max(abs(dY2(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = tnnX+lambda*comp_loss(E,loss);
            err = sqrt(norm(dY1(:))^2+norm(dY2(:))^2);
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y1 = Y1 + mu*dY1;
    Y2 = Y2 + mu*dY2;
    mu = min(rho*mu,max_mu);
end
obj = tnnX+lambda*comp_loss(E,loss);
err = sqrt(norm(dY1(:))^2+norm(dY2(:))^2);

author: canyilu

repo: LibADMM-toolbox


text:
%DRAWMARKER  Draws a marker on a predefined position in an image
%
%     img = cv.drawMarker(img, pos)
%     img = cv.drawMarker(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __img__ Input image.
% * __pos__ The point where the crosshair is positioned `[x,y]`.
%
% ## Output
% * __img__ Output image.
%
% ## Options
% * __Color__ Line color. Default is a black color.
% * __MarkerType__ The specific type of marker you want to use. One of:
%   * __Cross__ (`+`) (default) A crosshair marker shape.
%   * __TiltedCross__ (`x`) A 45 degree tilted crosshair marker shape.
%   * __Star__ (`*`) A star marker shape, combination of cross and tilted
%     cross.
%   * __Diamond__ (`d`) A diamond marker shape.
%   * __Square__ (`s`) A square marker shape.
%   * __TriangleUp__ (`^`) An upwards pointing triangle marker shape.
%   * __TriangleDown__ (`v`) A downwards pointing triangle marker shape.
% * __MarkerSize__ The length of the marker axis. Default is 20 pixels.
% * __Thickness__ Line thickness. default 1
% * __LineType__ Type of the line. One of:
%   * __4__ 4-connected line
%   * __8__ 8-connected line (default)
%   * __AA__ anti-aliased line
%
% The function cv.drawMarker draws a marker on a given position in the image.
% For the moment several marker types are supported, see `MarkerType` for more
% information.
%
% See also: cv.line, cv.circle
%

author: kyamagu

repo: mexopencv


text:
function [KE, PE] = energy(q,dq,p)
%[KE, PE] = energy(q,dq,p)
%
% This function computes the mechanical energy for the five-link biped
%
% INPUTS:
%   q = [5,n] = configuration
%   dq = [5,n] = rates
%   p = parameter struct
%
% OUTPUTS:
%   Fx = [1,n] = horizontal contact force acting on robot
%   Fy = [1,n] = vertical contact force acting on robot
%

[KE, PE] = autoGen_energy(...
 q(1,:),q(2,:),q(3,:),q(4,:),q(5,:),...
    dq(1,:),dq(2,:),dq(3,:),dq(4,:),dq(5,:),...
    p.m1, p.m2, p.m3, p.m4, p.m5, p.I1, p.I2, p.I3, p.I4, p.I5, p.l1, p.l2, p.l3, p.l4, p.c1, p.c2, p.c3, p.c4, p.c5, p.g);

end

author: MatthewPeterKelly

repo: OptimTraj