text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [threshold] = get_evaluation_threshold(appliance, household)

  thresholdMatrix = [15 15 0 15 15 15;        %fridge
                     15 15 0 15 15 15;        %freezer
                     15 15 0 50 50 15;        %microwave
                     15 15 0 15 15 15;        %dishwasher
                     15 15 0 15 15 15;        %entertainment
                     15 15 0 15 15 15;        %water kettle
                     15 15 0 15 15 15;        %cooker
                     15 15 0 15 15 15;        %coffee machine
                     15 15 0 15 15 15;        %washing machine
                     15 15 0 15 15 15;        %dryer
                     15 15 0 15 15 15;        %lamp
                     15 15 0 15 15 15;        %pc
                     15 15 0 15 15 15;        %laptop
                     15 15 0 15 15 15;        %tv
                     15 15 0 15 15 15;        %stereo
                     15 15 0 15 15 15;        %tablet
                     15 15 0 15 15 15;        %router
                     15 15 0 15 15 15];       %illuminated fountain

    threshold = thresholdMatrix(appliance, household);

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

htheta = sigmoid(X * theta);
J = 1 / m * sum(-y .* log(htheta) - (1 - y) .* log(1 - htheta));
for i = 1:size(theta, 1)
    grad(i) = 1 / m * sum((htheta - y) .* X(:, i));
end

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

end

author: merwan

repo: ml-class

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


text:
function B = bernstein_eval(i,n,u)
  % BERNSTERIN_EVAL evaluate a Bernstein polynomial. B_i,n(u)
  %
  % B = bernstein_eval(i,n,u)
  %
  % Inputs:
  %  i  #i scalar index of polynomial
  %  n  scalar order of polynomial
  %  u  #u evaluation points
  % Outputs:
  %  B  #u by #i values
  %
  u = reshape(u,numel(u),1);
  i = reshape(i,1,numel(i));
  nck = factorial(n)./(factorial(i).*factorial(n-(i)));
  B = nck .*  u.^i .* (1-u).^(n-i);
end

author: alecjacobson

repo: gptoolbox


text:
function clipboardData = clipboardpaste(format)
    %clipboardpaste v0.1.3
    %   Usage:
    %       clipboardData = clipboardpaste(format)
    %
    %   Input Arguments:
    %       format (optional)
    %           Requested image format. If the clipboard contains an image,
    %           then the image will be automatically set to format. Valid
    %           formats are:
    %               'png'
    %               'jpg' (or 'jpeg')
    %               'tif' (or 'tiff')
    %               'bmp'
    %               'gif'
    %               'ico'
    %               'cur'
    %               'pnm'
    %               'ppm'
    %               'pbm'
    %               'ras'
    %               'xwd'
    %               'fts'/'fits'
    %
    %           By default, the format is 'png'.
    %
    %   Output Arguments:
    %       clipboardData
    %           Struct representing the data on the clipboard. It has three
    %           fields:
    %               (1) primaryType: String denoting which kind of data was
    %               retrieved from the clipboard. It can be one of several
    %               possibilities:
    %                   -- 'text'
    %                   -- 'binary'
    %
    %               (2) subType: String denoting specifically which kind of
    %               data was retrieved from the clipboard. It can be one of
    %               several possibilities:
    %                   -- 'url'
    %                   -- 'image'
    %                   -- 'plain-text'
    %                   -- 'email'
    %
    %               (3) data: Data from the clipboard. This can be of
    %               several types, depending on the primaryType/subType
    %               combination:
    %                   -- 'text'/'url': A string denoting a URL
    %                   -- 'text'/'email': A string denoting an email
    %                   address
    %                   -- 'text'/'plain-text': A string
    %                   -- 'binary'/'image': An image of the specified
    %                   format (uint8 array). It is viewable by using the
    %                   command imshow()
    %
    %   Description:
    %       Retrieves data from the clipboard. The clipboard data may be of
    %       several types and it is encapsulated in a self-describing
    %       structure. This function is based on clipboardimage, by Saurabh
    %       Kumar, saurabhkumar_@rediffmail.com, with heavy modifications to
    %       address various cross-platform issues.
    %
    %   Example:
    %       First, copy something on to the clipboard. Say that you copy
    %       the number 14584 onto the clipboard.
    %       >> clipData = clipboardpaste %now paste it into the workspace
    %       clipData =
    %
    %                  data: '14584'
    %           primaryType: 'text'
    %               subType: 'plain-text'
    %       >> data = clipData.data; %get the actual data from the clipboard
    %
    %   =======================
    %   Written by Bryan Raines on May 5, 2008
    %   Last updated on May 6, 2008.
    %   Email: rainesb@ece.osu.edu
    %
    %   See also clipboard.

    usesXClipboard = isunix && ~strcmpi(computer, 'mac');

    if usesXClipboard
        if isempty(get(0,'DefaultFigureXDisplay'))
            error('MATLAB:clipboard:NoXDisplay', 'There is no X display set.');
        end
    else
        err = javachk('awt', 'Clipboard access');
        if (~isempty(err))
            error('MATLAB:clipboard:UnsupportedPlatform', err.message);
        end
    end

    clipboardData     = [];
    error(nargchk(0,1,nargin));

    if nargin < 1
        format = 'png';
    end

    try
        tKit                =       java.awt.Toolkit.getDefaultToolkit();
        clipboardHandle     =       tKit.getSystemClipboard();
        reqObj              =       java.lang.Object;
        clipboardContents   =       clipboardHandle.getContents(reqObj); %get all objects on clipboard
        Dflavor             =       clipboardContents.getTransferDataFlavors();

        %Go through the various DataFlavors until you find one that is
        %recognized
        for ii = 1:length(Dflavor)
            dataType = char(Dflavor(ii).getMimeType);

            switch dataType
                case 'image/x-pict; class=java.io.InputStream'
                    imarray = clipboardContents.getTransferData(Dflavor(ii));
                    clipboardData = encapsulateImageData(getImageFromClipboardData(imarray,format));
                    break;

                case 'image/x-java-image; class=java.awt.Image'
                    imarray = clipboardContents.getTransferData(Dflavor(ii));
                    clipboardData = encapsulateImageData(getImageFromClipboardData(imarray,format));
                    break;

                case 'application/x-java-file-list; class=java.util.List'
                    fileList = clipboardContents.getTransferData(Dflavor(ii));
                    clipboardData = encapsulateFileListData(fileList);
                    break;

                case 'application/x-java-url; class=java.net.URL'
                    urlText = char(clipboardContents.getTransferData(Dflavor(end))');
                    urlText = strtrim(urlText);

                    if isURLType(urlText)
                        if isImageType(urlText)
                            clipboardData = encapsulateImageData(downloadImage(urlText,format));
                        else
                            clipboardData = encapsulateURLData(urlText);
                        end
                        break;
                    end

                case 'application/x-java-serialized-object; class=java.lang.String'
                    textData = char(clipboardContents.getTransferData(Dflavor(ii)));
                    textData = strtrim(textData);

                    if isURLType(textData)
                        clipboardData = encapsulateURLData(textData);
                    elseif isEmailType(textData)
                        [tf emailAddress] = isEmailType(textData);
                        clipboardData = encapsulateEmailData(emailAddress);
                    else
                        clipboardData = encapsulatePlainTextData(textData);
                    end

                    break;
            end
        end
    catch
        disp('Error occurred in retrieving clipboard data');
    end
end

%=========== Encapsulate Clipboard Data in a clipboardStruct ===========
function clipboardData = encapsulateImageData(data)
    clipboardData.data = data;
    clipboardData.primaryType = 'binary';
    clipboardData.subType = 'image';
end

function clipboardData = encapsulateURLData(data)
    clipboardData.data = data;
    clipboardData.primaryType = 'text';
    clipboardData.subType = 'url';
end

function clipboardData = encapsulatePlainTextData(data)
    clipboardData.data = data;
    clipboardData.primaryType = 'text';
    clipboardData.subType = 'plain-text';
end

function clipboardData = encapsulateEmailData(data)
    clipboardData.data = data;
    clipboardData.primaryType = 'text';
    clipboardData.subType = 'email';
end

function clipboardData = encapsulateFileListData(data)
    fileListArray = data.toArray;

    for ii = 1:length(fileListArray)
        fileList{ii,1} = char(fileListArray(ii).toString);
    end

    clipboardData.data = fileList;
    clipboardData.primaryType = 'text';
    clipboardData.subType = 'file-list';
end

%=========== Determine Clipboard Data Type ===========
function tf = isImageType(textData)
    [var1 var2 ext] = fileparts(textData); %determine the file type
    imageTypes = {'.png','.jpg','.tif','.tiff','.bmp','.gif'};
    tf = any(strcmpi(imageTypes,ext));
end

function tf = isURLType(textData)
    %if the text just starts with www, then attach http:// as a prefix
    if ~isempty(regexp(textData,'^www','once'))
        textData = ['http://' textData];
    end

    %If the text starts with a protocol, then assume that it is a URL
    protocols = {'http','ftp','https','vnc','smb','file','irc'};

    for ii = 1:length(protocols)
        tf(ii) = ~isempty(regexp(textData,['^' protocols{ii}], 'once' ));
    end

    tf = any(tf);
end

function [tf emailAddress] = isEmailType(textData)
    emailAddress = regexp(textData,'^(mailto:)*(?<email>[\w\-\_\.]+@[\w\.\-\_]+)','tokens','once');
    tf = ~isempty(emailAddress);
    if tf
        emailAddress = emailAddress{2};
    end
end

%=========== Get Image from Clipboard ===========
function clipboardData = getImageFromClipboardData(imarr,format)
    %Grabs an image from the clipboard
    filename = [tempname '.' format];

    if ismac
        theImage = getImageFromPictStream(imarr);
    else
        theImage = imarr;
    end
    filehandle = java.io.File(filename);
    javax.imageio.ImageIO.write(theImage,format,filehandle);
    clipboardData = imread(filename);
    delete(filename);
end

function clipboardData = downloadImage(urlText,format)
    %Grabs an image by downloading it using the URL from the clipboard
    imarray = imread(urlText);
    filename = [tempname '.' format];
    imwrite(imarray,filename);
    clipboardData = imread(filename);
    delete(filename);
end

function awtImage = getImageFromPictStream(inputstream)
    %Exists for Mac OS X compabilibility. Requires Quicktime for Java to be
    %properly installed.

    %Check if QT for Java is installed
    try
        quicktime.QTSession.isInitialized();
    catch
        warning('Quicktime for Java is not installed.'); %#ok<WNTAG>
    end

    baos = java.io.ByteArrayOutputStream();

    isc = com.mathworks.mlwidgets.io.InterruptibleStreamCopier.getInterruptibleStreamCopier;
    isc.copyStream(inputstream,baos);
    baos.close;

    imgBytes = typecast(baos.toByteArray,'uint8');
    imgBytes = [cast(zeros(512,1),'uint8'); imgBytes];

    if ~quicktime.QTSession.isInitialized()
        quicktime.QTSession.open();
    end

    handle = quicktime.util.QTHandle(imgBytes);
    gi = quicktime.std.image.GraphicsImporter(quicktime.util.QTUtils.toOSType('PICT'));
    gi.setDataHandle(handle);
    qdRect = gi.getNaturalBounds();
    gid = quicktime.app.view.GraphicsImporterDrawer(gi);
    qip = quicktime.app.view.QTImageProducer(gid,java.awt.Dimension(qdRect.getWidth(),qdRect.getHeight()));
    awtImage = java.awt.Toolkit.getDefaultToolkit.createImage(qip);
    awtImage = getBufImage(awtImage);
end

function bufImage = getBufImage(clipboardContents)
    %Creates a rendered (buffered) image from an OSX image data flavor
    ww = clipboardContents.getWidth;
    hh = clipboardContents.getHeight;

    bufImage = java.awt.image.BufferedImage(ww,hh,java.awt.image.BufferedImage.TYPE_INT_RGB);

    g2 = bufImage.createGraphics;
    g2.setRenderingHint(java.awt.RenderingHints.KEY_INTERPOLATION,...
        java.awt.RenderingHints.VALUE_INTERPOLATION_BILINEAR);
    g2.drawImage(clipboardContents,0,0,ww,hh,[]);
    g2.dispose();
end

function tf = ismac
    %Detect if the computer is a Mac for versions before R2007a
    tf = strcmpi(computer,'mac');
end

author: Sable

repo: mcbench-benchmarks


text:
function out = dfun(in)

global X

out = L2_distance(X,X(:,in));
out = out';

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% read_eeglabevent() - import EEGLAB dataset events
%
% Usage:
%     >> event = read_eeglabevent(filename, ...);
%
% Inputs:
%     filename - [string] file name
%
% Optional inputs:
%     'header' - FILEIO structure header
%
% Outputs:
%     event - FILEIO toolbox event structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2008-

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

function event = read_eeglabevent(filename, varargin)

if nargin < 1
  help read_eeglabheader
  return
end

hdr = ft_getopt(varargin, 'header');

if isempty(hdr)
  hdr = read_eeglabheader(filename);
end

event    = [];                % these will be the output in FieldTrip format
oldevent = hdr.orig.event;    % these are in EEGLAB format

if ~isempty(oldevent)
  nameList = fieldnames(oldevent);
else
  nameList = [];
end
nameList = setdiff(nameList, {'type','value','sample','offset','duration','latency'});

for index = 1:length(oldevent)

  if isfield(oldevent,'code')
    type = oldevent(index).code;
  elseif isfield(oldevent,'value')
    type = oldevent(index).value;
  else
    type = 'trigger';
  end

  % events can have a numeric or a string value
  if isfield(oldevent,'type')
    value  = oldevent(index).type;
  else
    value = 'default';
  end

  % this is the sample number of the concatenated data to which the event corresponds
  sample = oldevent(index).latency;

  % a non-zero offset only applies to trial-events, i.e. in case the data is
  % segmented and each data segment needs to be represented as event. In
  % that case the offset corresponds to the baseline duration (times -1).
  offset = 0;

  if isfield(oldevent, 'duration')
    duration = oldevent(index).duration;
  else
    duration = 0;
  end

  % add the current event in FieldTrip format
  event(index).type     = type;     % this is usually a string, e.g. 'trigger' or 'trial'
  event(index).value    = value;    % in case of a trigger, this is the value
  event(index).sample   = sample;   % this is the sample in the datafile at which the event happens
  event(index).offset   = offset;   % some events should be represented with a shifted time-axix, e.g. a trial with a baseline period
  event(index).duration = duration; % some events have a duration, such as a trial

  %add custom fields
  for iField=1:length(nameList)
    eval(['event(index).' nameList{iField} '=oldevent(index).' nameList{iField} ';']);
  end

end % for oldevent

if hdr.nTrials>1
  % add the trials to the event structure
  for i=1:hdr.nTrials
    event(end+1).type     = 'trial';
    event(end  ).sample   = (i-1)*hdr.nSamples + 1;
    if isfield(oldevent,'setname') && (length(oldevent) == hdr.nTrials)
      event(end  ).value  = oldevent(i).setname; % accommodate Widmann's pop_grandaverage function
    else
      event(end  ).value  = [];
    end
    event(end  ).offset   = -hdr.nSamplesPre;
    event(end  ).duration =  hdr.nSamples;
  end
end % if nTrials>1

author: spm


text:
function [which] = mne_find_channel(info,name)
%
% [which] = mne_find_channel(info,name)
%
% Find a channel by name employing the info structure
% output by mne_raw2mat or mne_epochs2mat
%
% epoch - The data structure containing the channel information
% name  - name of the channel to look for
%
% Returns index of the channel in the data
% If the channel is not found, returns -1
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%     $Id: mne_find_channel.m 8776 2013-11-14 09:04:48Z roboos $
%
%     Revision 1.4  2006/04/23 15:29:40  msh
%     Added MGH to the copyright
%
%     Revision 1.3  2006/04/14 15:49:49  msh
%     Improved the channel selection code and added ch_names to measurement info.
%
%     Revision 1.2  2006/04/10 23:26:54  msh
%     Added fiff reading routines
%
%     Revision 1.1  2006/02/20 15:45:05  msh
%     Added mne_find_channel.m and mne_read_epoch.m
%
%
me='MNE:mne_find_channel';
if(nargin ~= 2)
   error(me,'Usage : [which] = mne_find_channel(info,which)');
end

which = strmatch(name,info.ch_names,'exact');

return;


text:
% op_makePhaseDrift.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,phDrift]=op_makePhaseDrift(in,totalDrift,noise);
%
% DESCRIPTION:
% Add phase drift to a dataset containing multiple averages.  This is generally
% used to generate simulated datasets with phase drift.
%
% INPUTS:
% in         = input data in matlab structure format.
% totalDrift = total amount of phase drift (in degrees) to add over the whole scan.
%               If totalDrift is a scalar, then a constant slope of drift
%               will be added.  If totalDrift is a vector or matrix with dimensions
%               equal to the dimensions of the input data, then this
%               vector specifies the drift applied to each average.
% noise      = the standard deviation of noise to add to the phase drift
%             function.
%
% OUTPUTS:
% out        = Output dataset with phase drift added.
% phDrift     = Vector of phase drift values that were added (in degrees).

function [out,phDrift]=op_makePhaseDrift(in,totalDrift,noise);
%out=op_makedrift(in,totalDrift);

%First make the matrices needed for multiplication
if any(totalDrift)
    if isequal(size(totalDrift),[in.averages/in.subspecs,in.subspecs])
        ph=totalDrift;
    else
        ph=linspace(0,totalDrift,in.averages*in.subspecs);
        ph=reshape(ph,in.subspecs,in.averages);
        ph=ph';
    end
else
    ph=zeros(in.averages,in.subspecs);
end
phnoise=noise*randn(size(ph));
phDrift=ph+phnoise;
%PH=repmat(phDrift',in.sz(1),1);
phDrift_shft=shiftdim(phDrift,-1);
PH=repmat(phDrift_shft,in.sz(1),1,1);

%Now apply the drift to the fids;
fids=in.fids.*exp(-1i*PH*pi/180);

%Now re-calculate specs using ifft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%FILLING IN DATA STRUCTURES
out=in;
out.fids=fids;
out.specs=specs;

%FILLING IN THE FLAGS
out.flags=in.flags;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%MINMAXSUBD Calculate subdomain minima and maxima.
%
%   [ MIN_VAL, MAX_VAL, MIN_COORD, MAX_COORD ] = MINMAXSUBD( S_EXPR, PROB, IND_S, IND_C, I_CUB, SOLNUM )
%   Evaluates the mimimum and maximum value of expression S_EXPR over the
%   subdomains indicated in IND_S or alternatively the cells in IND_C.
%   PROB is a valid finite element problem struct. Returns the minima and
%   maxima in MIN_VAL and MAX_VAL, and the corresponding coordinates in
%   MIN_COORD and MAX_COORD.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       s_expr      string                 Expression to evaluate
%       prob        struct                 Finite element problem struct
%       ind_s       [1,n_subd]             Subdomain numbers (default all)
%       ind_c       [1,n_cells]            Cell indices (default all)
%       i_cub       scalar                 Evaluation point rule (default 2)
%       solnum      scalar {n_sols}        Solution number/time to evaluate
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       min_val     scalar                 Minimum value of expression
%       max_val     scalar                 Maximum value of expression
%
%   See also MINMAXBDR, INTSUBD, INTBDR

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
classdef MOEADM2M_F5 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing MOEA/D-M2M

%------------------------------- Reference --------------------------------
% H. Liu, F. Gu, and Q. Zhang, Decomposition of a multiobjective
% optimization problem into a number of simple multiobjective subproblems,
% IEEE Transactions on Evolutionary Computation, 2014, 18(3): 450-455.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = X(:,2:end) - repmat(sin(pi/2*X(:,1)),1,size(X,2)-1);
            g = 2*abs(cos(pi*X(:,1))).*sum(-0.9*t.^2+abs(t).^0.6,2);
            PopObj(:,1) = (1+g).*X(:,1);
            PopObj(:,2) = (1+g).*(1-sqrt(X(:,1)));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function u=glotros(d,t,p)
%GLOTROS  Rosenberg glottal model U=(D,T,P)
% d is derivative of flow waveform
% t is in fractions of a cycle
% p has parameters
%	p(1)=closure time
%	p(2)=+ve/-ve slope ratio

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: glotros.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 2
  tt=(0:99)'/100;
else
  tt=mod(t,1);
end
u=zeros(size(tt));
de=[0.6 0.5]';
if nargin < 3
  p=de;
elseif length(p)<2
  p=[p(:); de(length(p)+1:2)];
end
pp=p(1)/(1+p(2));
ta=tt<pp;
tb=tt<p(1) & ~ta;
wa=pi/pp;
wb=0.5*pi/(p(1)-pp);
fb=wb*pp;
if d==0
  u(ta)=0.5*(1-cos(wa*tt(ta)));
  u(tb)=cos(wb*tt(tb)-fb);
elseif d==1
  u(ta)=0.5*wa*sin(wa*tt(ta));
  u(tb)=-wb*sin(wb*tt(tb)-fb);
elseif d==2
  u(ta)=0.5*wa^2.*cos(wa*tt(ta));
  u(tb)=-wb^2*cos(wb*tt(tb)-fb);
else
  error('Derivative must be 0,1 or 2');
end

author: decouples


text:
%Frame Coordinate frame object
%
% F = Frame(P, OPTIONS) creates an object that graphically renders
% a coordinate frame for SE(2), SO(2) or SE(3) represented by the
% pose P which can be:
% - homogeneous transform (3x3) for SE(2)
% - Quaternion for SO(3)
% - orthonormal rotation matrix (3x3) for SO(3)
% - homogeneous transform (4x4) for SE(3)
%
% Methods::
% move       move the graphical coordinate frame to a new pose
% animate    move the graphical coordinate frame to a new pose
% char
% display
% delete
%
% Options::
% 'color',C          The color to draw the axes, MATLAB colorspec C
% 'noaxes'           Don't display axes on the plot
% 'axis',A           Set dimensions of the MATLAB axes to A=[xmin xmax ymin ymax zmin zmax]
% 'frame',F          The frame is named {F} and the subscript on the axis labels is F.
% 'text_opts', opt   A cell array of MATLAB text properties
% 'handle',H         Draw in the MATLAB axes specified by the axis handle H
% 'view',V           Set plot view parameters V=[az el] angles, or 'auto'
%                    for view toward origin of coordinate frame
% 'arrow'            Use arrows rather than line segments for the axes
% 'width', w         Width of arrow tips
%
% Examples::
%
%       f_a = Frame(TA, 'frame', 'A')
%       f_b = Frame(TB, 'frame', 'B', 'color', 'b')
%       f_c = Frame(TC, 'frame', 'C', 'text_opts', {'FontSize', 10, 'FontWeight', 'bold'})
%
%       f_a.move(T);
%
% Notes::
% - The arrow option requires the third party package arrow3.
%
% See also TRPLOT2, TRANIMATE.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% TODO
%  need to decide how to handle scaling
%  what does hold on mean?  don't touch scaling?

classdef Frame < handle

    properties
        T
        se2
        name
        hg
    end

    methods

        function f = Frame(T, varargin)

            if size(T,3) > 1
                error('trplot cannot operate on a sequence');
            end

            if ~ishomog(T) && ~isrot(T)
                error('trplot operates only on transform (4x4) or rotation matrix (3x3)');
            end

            opt.color = 'b';
            opt.axes = true;
            opt.axis = [];
            opt.frame = [];
            opt.text_opts = [];
            opt.view = [];
            opt.width = 1;
            opt.arrow = false;
            opt.handle = [];
            opt.se2 = false;

            opt = tb_optparse(opt, varargin);

            f.se2 = opt.se2;
            f.name = opt.frame;

            % axis labels
            if isempty(opt.frame)
                fmt = '%c';
            else
                fmt = sprintf('%%c_{%s}', opt.frame);
            end

            % text label options
            if isempty(opt.text_opts)
                opt.text_opts = {};
            end

            if isempty(opt.axis)
                % determine some default axis dimensions

                d = 1.2;
                if opt.se2
                    c = transl(T);
                    d = 1.2;
                    opt.axis = [c(1)-d c(1)+d c(2)-d c(2)+d];
                else
                    % get the origin of the frame
                    if isrot(T)
                        c = [0 0 0];  % at zero for a rotation matrix
                    else
                        c = transl(T);
                    end
                    opt.axis = [c(1)-d c(1)+d c(2)-d c(2)+d c(3)-d c(3)+d];
                end
            end

            % create the axes
            if ~isempty(opt.handle)
                hax = opt.handle;
                hold(hax);
            else
                ih = ishold;
                if ~ih
                    % if hold is not on, then clear the axes and set scaling
                    cla
                    if ~isempty(opt.axis)
                        axis(opt.axis);
                    end
                    daspect([1 1 1]);

                    if opt.axes
                        xlabel( 'X');
                        ylabel( 'Y');
                        zlabel( 'Z');
                        rotate3d on
                    end
                    new_plot = true;
                end
                hax = gca;
                hold on
            end

            opt.text_opts = {opt.text_opts{:}, 'Color', opt.color};

            % create the transfor for the frame, this allows the whole
            % graphical structure to be easily moved
            hg = hgtransform('Parent', hax);
            f.hg = hg;
            set(f.hg, 'Tag', 'Frame');
            set(f.hg, 'UserData', f);

            if opt.se2
                % create unit vectors
                o =  [0 0]';
                x1 = [1 0]';
                y1 = [0 1]';

                % draw the axes
                mstart = [o o]';
                mend = [x1 y1]';

                if opt.arrow
                    % draw the 2 arrows
                    S = [opt.color num2str(opt.width)];
                    ha = arrow3(mstart, mend, S);
                    for h=ha'
                        set(h, 'Parent', hg);
                    end
                else
                    for i=1:2
                        plot2([mstart(i,1:2); mend(i,1:2)], ...
                            'Color', opt.color, 'Parent', hg);
                    end
                end

                % add the labels to each axis
                h = text(x1(1), x1(2), sprintf(fmt, 'X'), 'Parent', hg);
                if ~isempty(opt.text_opts)
                    set(h, opt.text_opts{:});
                end

                h = text(y1(1), y1(2), sprintf(fmt, 'Y'), 'Parent', hg);
                if ~isempty(opt.text_opts)
                    set(h, opt.text_opts{:});
                end

                % label the frame
                if ~isempty(opt.frame)
                    h = text(o(1)-0.04*x1(1), o(2)-0.04*y1(2), ...
                        ['\{' opt.frame '\}'], 'Parent', hg);
                    set(h, 'VerticalAlignment', 'middle', ...
                        'HorizontalAlignment', 'center', opt.text_opts{:});
                end
            else
                % create unit vectors
                o =  [0 0 0]';
                x1 = [1 0 0]';
                y1 = [0 1 0]';
                z1 = [0 0 1]';

                % draw the axes
                mstart = [o o o]';
                mend = [x1 y1 z1]';

                if opt.arrow
                    % draw the 3 arrows
                    S = [opt.color num2str(opt.width)];
                    ha = arrow3(mstart, mend, S);
                    for h=ha'
                        set(h, 'Parent', hg);
                    end
                else
                    for i=1:3
                        h = plot2([mstart(i,1:3); mend(i,1:3)], ...
                            'Color', opt.color, 'Parent', hg);
                    end
                end

                % add the labels to each axis
                h = text(x1(1), x1(2), x1(3), sprintf(fmt, 'X'), 'Parent', hg);
                set(h, opt.text_opts{:});

                h = text(y1(1), y1(2), y1(3), sprintf(fmt, 'Y'), 'Parent', hg);
                set(h, opt.text_opts{:});

                h = text(z1(1), z1(2), z1(3), sprintf(fmt, 'Z'), 'Parent', hg);
                set(h, opt.text_opts{:});

                % label the frame
                if ~isempty(opt.frame)
                    h = text(o(1)-0.04*x1(1), o(2)-0.04*y1(2), o(3)-0.04*z1(3), ...
                        ['\{' opt.frame '\}'], 'Parent', hg);
                    set(h, 'VerticalAlignment', 'middle', ...
                        'HorizontalAlignment', 'center', opt.text_opts{:});
                end
            end

            if ~opt.axes
                set(gca, 'visible', 'off');
            end
            if isstr(opt.view) && strcmp(opt.view, 'auto')
                cam = x1+y1+z1;
                view(cam(1:3));
            elseif ~isempty(opt.view)
                view(opt.view);
            end

            if isempty(opt.handle) && ~ih
                grid on
                hold off
            end

            % now place the frame in the desired pose
            f.move(T);
        end

        function move(f, T)
            if f.se2
                if ~all(size(T) == [3 3])
                    error('expecting SE(2) matrix');
                end
                T = [T(1:2,1:2) zeros(2,1) T(1:2,3); 0 0 1 0; 0 0 0 1];
            elseif isrot(T)
                T = r2t(T);
            elseif ~ishomog(T)
                error('expecting SO(3) or SE(3) matrix');
            end

            % search for this named frame in all figs
            set(f.hg, 'Matrix', T);
            f.T = T;
        end

        function animate(f, P2, varargin)
        %ANIMATE Animate a coordinate frame
        %
        % ANIMATE(P1, P2, OPTIONS) animates a 3D coordinate frame moving from pose P1
        % to pose P2.  Poses P1 and P2 can be represented by:
        %   - homogeneous transformation matrices (4x4)
        %   - orthonormal rotation matrices (3x3)
        %   - Quaternion
        %
        % ANIMATE(P, OPTIONS) animates a coordinate frame moving from the identity pose
        % to the pose P represented by any of the types listed above.
        %
        % ANIMATE(PSEQ, OPTIONS) animates a trajectory, where PSEQ is any of
        %   - homogeneous transformation matrix sequence (4x4xN)
        %   - orthonormal rotation matrix sequence (3x3xN)
        %   - Quaternion vector (Nx1)
        %
        % Options::
        %  'fps', fps    Number of frames per second to display (default 10)
        %  'nsteps', n   The number of steps along the path (default 50)
        %  'axis',A      Axis bounds [xmin, xmax, ymin, ymax, zmin, zmax]
        %
        % See also TRPLOT.

            opt.fps = 10;
            opt.nsteps = 50;
            opt.axis = [];

            [opt, args] = tb_optparse(opt, varargin);

            P1 = [];

            % convert quaternion and rotation matrix to hom transform
            if isa(P2, 'Quaternion')
                T2 = P2.T;   % convert quaternion to transform
                if ~isempty(args) && isa(args{1},'Quaternion')
                    P1 = T2;
                    Q2 = args{1};
                    T2 = Q2.T;
                    args = args(2:end);
                else
                    T1 = eye(4,4);
                end
            elseif isrot(P2)
                T2 = r2t(P2);
                if ~isempty(args) && isrot(args{1})
                    P1 = T2;
                    T2 = r2t(args{1});
                    args = args(2:end);
                else
                    T1 = eye(4,4);
                end
            elseif ishomog(P2)
                T2 = P2;
                if ~isempty(args) && ishomog(args{1})
                    P1 = T2;
                    T2 = args{1};
                    args = args(2:end);
                else
                    T1 = eye(4,4);
                end
            end

            % at this point
            %   T1 is the initial pose
            %   T2 is the final pose
            %
            %  T2 may be a sequence

            if size(T2,3) > 1
                % tranimate(Ts)
                % we were passed a homog sequence
                if ~isempty(P1)
                    error('only 1 input argument if sequence specified');
                end
                Ttraj = T2;
            else
                % tranimate(P1, P2)
                % create a path between them
                Ttraj = ctraj(T1, T2, opt.nsteps);
            end

            if isempty(opt.axis)
                % create axis limits automatically based on motion of frame origin
                t = transl(Ttraj);
                mn = min(t) - 1.5;  % min value + length of axis + some
                mx = max(t) + 1.5;  % max value + length of axis + some
                axlim = [mn; mx];
                axlim = axlim(:)';
                args = [args 'axis' axlim];
            end

            % animate it for all poses in the sequence
            for i=1:size(Ttraj,3)
                T = Ttraj(:,:,i);
                f.move(T);
                pause(1/opt.fps);
            end
        end

        function delete(f)
        % DELETE Delete the coordinate frame
            children = get(f.hg, 'Children');
            for child=children
                delete(child);
            end
        end

        function s = char(f)
        %Link.char String representation of parameters
        %
        % s = L.char() is a string showing link parameters in compact single line format.
        % If L is a vector of Link objects return a string with one line per Link.
        %
        % See also Link.display.
            ts = trprint(f.T, 'rpy');
            s = sprintf(' {%s} %s', f.name, ts);
            if f.se2
                s = [s ' :: SE(2)'];
            else
                s = [s ' :: SE(3)'];
            end
        end

        function display(l)
        %Frame.display Display parameters
        %
        % F.display() display link parameters in compact single line format.  If L is a
        % vector of Link objects display one line per element.
        %
        % Notes::
        % - this method is invoked implicitly at the command line when the result
        %   of an expression is a Link object and the command has no trailing
        %   semicolon.
        %
        % See also Link.char, Link.dyn, SerialLink.showlink.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(l) );
        end % display()

        function rescale(f)
            mn = [Inf Inf Inf];
            mx = -[Inf Inf Inf];
            for frame=findobj('Tag', 'Frame')'
                %%T = frame.

            end
        end

    end % methods
end

author: petercorke

repo: robotics-toolbox-matlab


text:
clear all;
close all;
clc;

targetFolder = './images/';
candidateImages = {'skewed1.gif', 'trapezoid.jpg', 'QR2.jpg', 'QR3.jpg', 'book.jpg', 'square.jpg', 'skewed1.gif' };
candidateChoice = 5;

targetImageFile = strcat(targetFolder , char(candidateImages(candidateChoice )));
targetImageData = imread(targetImageFile);

% 5.1 Undoing Perspective Distortion of Planar Surface
% a)

%book=imread('trapezoid.jpg');
imshow(targetImageData);

% b)
fprintf('Corner selection must be clockwise or anti-clockwise.\n');
[X Y] = ginput(4);

%if ispolycw(X, Y) % if the coordinates are not clockwise, sort them clockwise
%    [X, Y] = poly2cw(X, Y)
%end

%X = uint8(X);
%Y = uint8(Y);
[X Y] = sortPolyFromClockwiseStartingFromTopLeft( X, Y );

x=[1;210;210;1];
y=[1;1;297;297];

% c)
A=zeros(8,8);
A(1,:)=[X(1),Y(1),1,0,0,0,-1*X(1)*x(1),-1*Y(1)*x(1)];
A(2,:)=[0,0,0,X(1),Y(1),1,-1*X(1)*y(1),-1*Y(1)*y(1)];

A(3,:)=[X(2),Y(2),1,0,0,0,-1*X(2)*x(2),-1*Y(2)*x(2)];
A(4,:)=[0,0,0,X(2),Y(2),1,-1*X(2)*y(2),-1*Y(2)*y(2)];

A(5,:)=[X(3),Y(3),1,0,0,0,-1*X(3)*x(3),-1*Y(3)*x(3)];
A(6,:)=[0,0,0,X(3),Y(3),1,-1*X(3)*y(3),-1*Y(3)*y(3)];

A(7,:)=[X(4),Y(4),1,0,0,0,-1*X(4)*x(4),-1*Y(4)*x(4)];
A(8,:)=[0,0,0,X(4),Y(4),1,-1*X(4)*y(4),-1*Y(4)*y(4)];

v=[x(1);y(1);x(2);y(2);x(3);y(3);x(4);y(4)];

u=A\v;
%which reshape

U=reshape([u;1],3,3)';

w=U*[X';Y';ones(1,4)];
w=w./(ones(3,1)*w(3,:));

% d)
%which maketform
T=maketform('projective',U');

%which imtransform
P2=imtransform(targetImageData,T,'XData',[1 210],'YData',[1 297]);

% e)
figure,
imshow(P2)


text:
classdef FileFolderState
    %FILEFOLDERSTATE Represent file/folder selections
    %   Enumerates a list of choices

    % Copyright 2020-2021 The MathWorks Inc.

    %% Enumerations
    enumeration
        file
        folder
    end %enumeration

end %classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%                               stopNow.m
%
% This function is used in the main loop of many solvers (i.e.solve*.m) to
% check if the stopping condition(time, residual and reconstruction error)
% has been met and thus loop should be breaked.
%
%
% Note:
% This function does not check for max iterations since the for-loop
% in the solver already gurantee it.
%
% Inputs:
%         opts(struct)                   :  consists of options. It is as
%                                           defined in solverPhaseRetrieval.
%                                           See its header or User Guide
%                                           for details.
%         currentResid(real number)      :  Definition depends on the
%                                           specific algorithm used see the
%                                           specific algorithm's file's
%                                           header for details.
%         currentReconError(real number) :  norm(xt-x)/norm(xt), where xt
%                                           is the m x 1 true signal,
%                                           x is the n x 1 estimated signal
%                                           at current iteration.
% Outputs:
%         ifStop(boolean)                :  If the stopping condition has
%                                           been met.
%
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

function ifStop = stopNow(opts, currentTime, currentResid, currentReconError)
    if currentTime >= opts.maxTime
        ifStop = true;
        return
    end
    if ~isempty(opts.xt)
            assert(~isempty(currentReconError),'If xt is provided, currentReconError must be provided.');
            ifStop = currentReconError < opts.tol;
    else
        assert(~isempty(currentResid),'If xt is not provided, currentResid must be provided.');
        ifStop = currentResid < opts.tol;
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
function [bbox_pred, bbox_pred_in_region, bbox_pred_in_region_quantized] = decode_loc_probs_to_bbox_targets(...
    bbox_in, class_indices, loc_prob_vectors, conf)
% decode_loc_probs_to_bbox_targets: given the input bounding boxes (bbox_in),
% the category ids of each bounding box (class_indices), and the predicted
% probability vectors (loc_prob_vectors) of each input bounding box it
% returns the location of the predicted bounding boxes.
%
% INPUTS:
% 1) bbox_in: a N x 4 array with the input bounding box coordinates in the
% form of [xi0,yi0,xi1,yi1] where (xi0,yi0) is the tot-left corner and
% (xi1,yi1) is the bottom-right corner. N is the number of bounding boxes.
% 2) class_indices: a N x 1 array with the category id of each input
% bounding box.
% 3) loc_prob_vectors: a M x K x C x N array with the predicted probability
% vectors of each input bounding box; N is the number of bounding boxes, C
% is the number of categoris, M is the resolution of the target
% probabilities specified by the input parameter conf.resolution and K is
% the number of target probability vectors per bounding box and per
% category. For the InOut model K = 2 (1 vector for each axis), for the
% Border model K = 4 (2 vectors for each axis; e.g. for the x axis we have
% 1 probability vector for the left border and 1 probability vector for the
% right bordr), and for the Combined model K = 6 (3 vectors for each axis;
% e.g. for the x axis we have 1 probability vector for the in-out elements,
% 1 probability vector for the left border, and 1 probability vector
% for the right border).
% 4) conf: struct with the configuration parameters of LocNet:
%   4.1) conf.scale_ratio: (scalar value) the scaling factor of the search region
%        w.r.t. the input bounding boxes. Specifically, given an input
%        bounding boxes, in order for LocNet to localize the target
%        bounding box it "look" in a search region that is obtained by scaling
%        the input bounding box by a scaling factor. The scaling factor is
%        given by the scale_ratio parameter.
%   4.2) conf.resolution: scalar value; In order for the LocNet to localize
%   the target bounding boxes inside the search region, it considers a division
%   of the search region in M horizontal stripes (rows) and M vertical stripes
%   (columns). The value of M is given by the parameter conf.resolution.
%   4.3) conf.num_classes: (scalar value) number of categories
%   4.4) conf.loc_type: string wiht the type of the localization model that
%   the LocNet implements. The possible options are: 'inout','borders', or 'combined'
%
% OUTPUTS:
% 1) bbox_pred: a N x 5 array with the predicted bounding box coordinates in the
% form of [xt0,yt0,xt1,yt1,c] where (xt0,yt0) is the tot-left corner, (xt1,yt1)
% is the bottom-right corner, and c is the category id of the bounding box.
% 3) bbox_pred_in_region: a N x 5 array with the predicted bounding box
% coordinates in the form of [xtr0,ytr0,xtr1,ytr1] where (xtr0,ytr0) is
% the tot-left corner, (xtr1,ytr1) is the bottom-right corner, and c is the
% category id of the bounding box. The coordinates are w.r.t. the
% corresponding search regions.
% 2) bbox_target_quantized: a N x 5 array with the predicted bounding box
% coordinates in the form of [xtq0,ytq0,xtq1,ytq1] where (xtq0,ytq0) is
% the tot-left corner, (xtq1,ytq1) is the bottom-right corner, and c is the
% category id of the bounding box. The coordinates are w.r.t. the
% corresponding search regions and quantized in discrite values between 1
% and resolution.
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

resolution  = conf.resolution;
scale_ratio = conf.scale_ratio;
loc_type    = conf.loc_type;

assert(size(loc_prob_vectors,1) == conf.resolution);

% get the coordinates of the search regions
bbox_search_region = scale_bboxes(bbox_in, scale_ratio);

% given the predicted probability vectors and the category ids, estimate
% the most likely bounding box location of the objects of interest; the
% predicted bounding boxes locations are w.r.t. the search region
% bbox_search_region and quantized in discrite values between 1 and
% resolution (in the paper the resolution parameter is referred as M)
switch loc_type
    case 'inout'
        bbox_pred_in_region_quantized = decode_in_out_probabilities(loc_prob_vectors, class_indices);
    case 'borders'
        bbox_pred_in_region_quantized = decode_borders_probabilities(loc_prob_vectors, class_indices);
    case 'combined'
        bbox_pred_in_region_quantized = decode_combined_probabilities(loc_prob_vectors, class_indices);
    otherwise
        error('Invalid localization type %s',loc_type)
end

% get the coordinates of the predicted bounding boxes w.r.t. the original
% image and in pixel units.
[bbox_pred, bbox_pred_in_region] = decode_quantized_location(...
    bbox_pred_in_region_quantized, resolution, bbox_search_region);
bbox_pred = [bbox_pred, class_indices];
end

function [real_coord, real_coord_region] = decode_quantized_location(quantized, resolution, bbox_search_region)
search_width  = bbox_search_region(:,3)-bbox_search_region(:,1)+1;
search_height = bbox_search_region(:,4)-bbox_search_region(:,2)+1;

% transform the coordinates of the predicted bounding boxes in pixel units
real_coord_region      = quantized;
real_coord_region(:,1) = (quantized(:,1)-1).*(search_width /resolution) + 0.5;
real_coord_region(:,2) = (quantized(:,2)-1).*(search_height/resolution) + 0.5;
real_coord_region(:,3) = (quantized(:,3)-1).*(search_width /resolution) + 0.5;
real_coord_region(:,4) = (quantized(:,4)-1).*(search_height/resolution) + 0.5;

% get the coordinates of the predicted bounding boxes w.r.t. the image
real_coord      = real_coord_region;
real_coord(:,1) = real_coord_region(:,1)+bbox_search_region(:,1);
real_coord(:,3) = real_coord_region(:,3)+bbox_search_region(:,1);
real_coord(:,2) = real_coord_region(:,2)+bbox_search_region(:,2);
real_coord(:,4) = real_coord_region(:,4)+bbox_search_region(:,2);
end

function best_location = decode_in_out_probabilities(loc_maps_per_cls, class_indices)

% keep from each input bounding box the probability vectors that correspond
% to its category id
loc_maps_size = size(loc_maps_per_cls);
assert(max(class_indices) <= loc_maps_size(3));
loc_maps_size(3) = 1;
loc_prob_vectors = zeros(loc_maps_size,'single');
unique_class_indices = unique(class_indices);
for i = 1:length(unique_class_indices)
    class_idx   = unique_class_indices(i);
    is_this_cls = class_idx == class_indices;
    loc_prob_vectors(:,:,1,is_this_cls) = loc_maps_per_cls(:,:,class_idx,is_this_cls);
end
loc_prob_vectors = squeeze(loc_prob_vectors);

% given the predicted probability vectors, for each input bounding box,
% estimate the most likely location (by minimizing the negative log likelihood)
% of bounding box of the actual object independently for the x and y axis.
[best_locationsX] = minimizeNegLogLikelihoodInOut(loc_prob_vectors(:,1,:)); % x axis
[best_locationsY] = minimizeNegLogLikelihoodInOut(loc_prob_vectors(:,2,:)); % y axis
best_location = single([best_locationsX(:,1),best_locationsY(:,1),best_locationsX(:,2),best_locationsY(:,2)]);
end

function best_location = decode_combined_probabilities(loc_maps_per_cls, class_indices)

% keep from each input bounding box the probability vectors that correspond
% to its category id
loc_maps_size = size(loc_maps_per_cls);
assert(max(class_indices) <= loc_maps_size(3));
loc_maps_size(3) = 1;
loc_prob_vectors = zeros(loc_maps_size,'single');
unique_class_indices = unique(class_indices);
for i = 1:length(unique_class_indices)
    class_idx   = unique_class_indices(i);
    is_this_cls = class_idx == class_indices;
    loc_prob_vectors(:,:,1,is_this_cls) = loc_maps_per_cls(:,:,class_idx,is_this_cls);
end
loc_prob_vectors = squeeze(loc_prob_vectors);

% given the predicted probability vectors, for each input bounding box,
% estimate the most likely location (by minimizing the negative log likelihood)
% of bounding box of the actual object independently for the x and y axis.
[best_locationsX] = minimizeNegLogLikelihoodCombined(loc_prob_vectors(:,1:3,:)); % x axis
[best_locationsY] = minimizeNegLogLikelihoodCombined(loc_prob_vectors(:,4:6,:)); % y axis
best_location = single([best_locationsX(:,1),best_locationsY(:,1),best_locationsX(:,2),best_locationsY(:,2)]);
end

function best_location = decode_borders_probabilities(loc_maps_per_cls, class_indices)

% keep from each input bounding box the probability vectors that correspond
% to its category id
loc_maps_size = size(loc_maps_per_cls);
assert(max(class_indices) <= loc_maps_size(3));
loc_maps_size(3) = 1;
loc_prob_vectors = zeros(loc_maps_size,'single');
unique_class_indices = unique(class_indices);
for i = 1:length(unique_class_indices)
    class_idx   = unique_class_indices(i);
    is_this_cls = class_idx == class_indices;
    loc_prob_vectors(:,:,1,is_this_cls) = loc_maps_per_cls(:,:,class_idx,is_this_cls);
end
loc_prob_vectors = squeeze(loc_prob_vectors);

% given the predicted probability vectors, for each input bounding box,
% estimate the most likely location (by minimizing the negative log likelihood)
% of bounding box of the actual object independently for the x and y axis.
[best_locationsX] = minimizeNegLogLikelihoodBorders(loc_prob_vectors(:,1:2,:)); % x axis
[best_locationsY] = minimizeNegLogLikelihoodBorders(loc_prob_vectors(:,3:4,:)); % y axis
best_location = single([best_locationsX(:,1),best_locationsY(:,1),best_locationsX(:,2),best_locationsY(:,2)]);
end

function [best_location, MinNegLogLikelihood, NegLogLikelihoodPerLocation, all_locations] = minimizeNegLogLikelihoodCombined(Probs)
% Given the predicted Combined probability vectors for a single dimension
% (either the x or y axis), estimate the most likely location of the actual
% bounding box in the corresponding axis by minimizing the negative
% log-likelihood. For the mimization we use exaustive search.

min_prob         = 0.0001;
PositiveProbs = max(  Probs,min_prob);
NegativeProbs = max(1-Probs,min_prob);

LogInsideProbs   = -log(squeeze(PositiveProbs(:,1,:)));
LogBordersProbs0 = -log(squeeze(PositiveProbs(:,2,:)))';
LogBordersProbs1 = -log(squeeze(PositiveProbs(:,3,:)))';

LogOutsideProbs     = -log(squeeze(NegativeProbs(:,1,:)));
LogNonBordersProbs0 = -log(squeeze(NegativeProbs(:,2,:)))';
LogNonBordersProbs1 = -log(squeeze(NegativeProbs(:,3,:)))';

resolution            = size(LogInsideProbs,1);
num_elems             = size(LogOutsideProbs,2);
LogInsideProbsCumSum  = cumsum([zeros(1,num_elems,'single');LogInsideProbs], 1)';
LogOutsideProbsCumSum = cumsum([zeros(1,num_elems,'single');LogOutsideProbs],1)';

NegLogLikelihoodPerLocation = zeros([num_elems, resolution*resolution],'single');
all_locations = zeros([resolution*resolution, 2],'single');
c = 0;
% compute the negative log-likelihood of each possible location of the
% bounding box
for a = 1:(resolution-1) % the double loop iterates over all possible (discrete) locations of the bounding box
    for b = (a+1):resolution
        c = c + 1;
        % NegLogLikelihoodPerLocation(:,c) is the negative log-likelihood of
        % the bounding box with coordines on [a, b] in the x (or y) axis.
        % For instance, in the x axis, a is the location of the left border
        % of the bounding box and b is the location of the right border of
        % the bounding box.
        NegLogLikelihoodPerLocation(:,c) = (LogInsideProbsCumSum(:,b+1)-LogInsideProbsCumSum(:,a))-...
            (LogOutsideProbsCumSum(:,b+1)-LogOutsideProbsCumSum(:,a)) + ...
            (LogBordersProbs0(:,a) + LogBordersProbs1(:,b)) - ...
            (LogNonBordersProbs0(:,a) + LogNonBordersProbs1(:,b));
        all_locations(c,1) = a;
        all_locations(c,2) = b;
    end
end
all_locations = all_locations(1:c,:);
NegLogLikelihoodPerLocation = NegLogLikelihoodPerLocation(:,1:c);

% pick the location with the minimum negative log-likelihood
[MinNegLogLikelihood, min_idx] = min(NegLogLikelihoodPerLocation,[],2);
best_location = all_locations(min_idx,:);
end

function [best_location, MinNegLogLikelihood, NegLogLikelihoodPerLocation, all_locations] = minimizeNegLogLikelihoodInOut(Probs)
% Given the predicted InOut probability vectors for a single dimension
% (either the x or y axis), estimate the most likely location of the actual
% bounding box in the corresponding axis by minimizing the negative
% log-likelihood. For the mimization we use exaustive search.

min_prob         = 0.0001;
InsideProbs      = max(  Probs,min_prob);
OutsideProbs     = max(1-Probs,min_prob);
LogInsideProbs   = -log(squeeze(InsideProbs(:,1,:)));
LogOutsideProbs  = -log(squeeze(OutsideProbs(:,1,:)));

resolution            = size(LogInsideProbs,1);
num_elems             = size(LogOutsideProbs,2);
LogInsideProbsCumSum  = cumsum([zeros(1,num_elems,'single');LogInsideProbs], 1)';
LogOutsideProbsCumSum = cumsum([zeros(1,num_elems,'single');LogOutsideProbs],1)';

NegLogLikelihoodPerLocation = zeros([num_elems, resolution*resolution],'single');
all_locations = zeros([resolution*resolution, 2],'single');
c = 0;
% compute the negative log-likelihood of each possible location of the
% bounding box
for a = 1:(resolution-1) % the double loop iterates over all possible (discrete) locations of the bounding box
    for b = (a+1):resolution
        c = c + 1;
        % NegLogLikelihoodPerLocation(:,c) is the negative log-likelihood of
        % the bounding box with coordines on [a, b] in the x (or y) axis.
        % For instance, in the x axis, a is the location of the left border
        % of the bounding box and b is the location of the right border of
        % the bounding box.
        NegLogLikelihoodPerLocation(:,c) = (LogInsideProbsCumSum(:,b+1)-LogInsideProbsCumSum(:,a))-...
            (LogOutsideProbsCumSum(:,b+1)-LogOutsideProbsCumSum(:,a));
        all_locations(c,1) = a;
        all_locations(c,2) = b;
    end
end
all_locations = all_locations(1:c,:);
NegLogLikelihoodPerLocation = NegLogLikelihoodPerLocation(:,1:c);

% pick the location with the minimum negative log-likelihood
[MinNegLogLikelihood, min_idx] = min(NegLogLikelihoodPerLocation,[],2);
best_location = all_locations(min_idx,:);
end

function [best_location, MinNegLogLikelihood, NegLogLikelihoodPerLocation, all_locations] = minimizeNegLogLikelihoodBorders(Probs)
% Given the predicted Borders probability vectors for a single dimension
% (either the x or y axis), estimate the most likely location of the actual
% bounding box in the corresponding axis by minimizing the negative
% log-likelihood. For the mimization we use exaustive search.

min_prob         = 0.0001;
PositiveProbs    = max(  Probs,min_prob);
LogBordersProbs0 = -log(squeeze(PositiveProbs(:,1,:)))';
LogBordersProbs1 = -log(squeeze(PositiveProbs(:,2,:)))';

resolution = size(LogBordersProbs0,2);
num_elems = size(LogBordersProbs0,1);

NegLogLikelihoodPerLocation = zeros([num_elems, resolution*resolution],'single');
all_locations = zeros([resolution*resolution, 2],'single');
c = 0;
% compute the negative log-likelihood of each possible location of the
% bounding box
for a = 1:(resolution-1) % the double loop iterates over all possible (discrete) locations of the bounding box
    for b = (a+1):resolution
        c = c + 1;
        % NegLogLikelihoodPerLocation(:,c) is the negative log-likelihood of
        % the bounding box with coordines on [a, b] in the x (or y) axis.
        % For instance, in the x axis, a is the location of the left border
        % of the bounding box and b is the location of the right border of
        % the bounding box.
        NegLogLikelihoodPerLocation(:,c) = (LogBordersProbs0(:,a) + LogBordersProbs1(:,b));
        all_locations(c,1) = a;
        all_locations(c,2) = b;
    end
end
all_locations = all_locations(1:c,:);
NegLogLikelihoodPerLocation = NegLogLikelihoodPerLocation(:,1:c);

% pick the location with the minimum negative log-likelihood
[MinNegLogLikelihood, min_idx] = min(NegLogLikelihoodPerLocation,[],2);
best_location = all_locations(min_idx,:);
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function selectInteractive(job,varargin)
% compute votes from grain boundaries
%
% Syntax
%
%   % compute votes from all p2c and c2c boundaries
%   job.calcGBVotes('threshold', 2*degree)
%
%   % compute votes only from p2c boundaries -> growth algorithm
%   job.calcGBVotes('p2c', 'threshold', 3*degree, 'tol', 1.5*degree)
%
% Input
%  job - @parentGrainReconstructor
%
% Output
%  job.votes - table of votes
%
% Options
%  p2c  - consider only parent / child grain boundaries
%  c2c  - consider only child / child grain boundaries
%  threshold - threshold fitting angle between job.p2c and the boundary OR
%  tolerance - range over which the probability increases from 0 to 1 (default 1.5)
%  numFit - number of fits to be computed
%

% datacursormode does not work with grains due to a Matlab bug
datacursormode off

% define a hand written selector
set(gcf,'WindowButtonDownFcn',{@doSelection,job});

end

function doSelection(src,eventdata,job)

% remove old selection
ax = gca;
handleSelected = getappdata(ax,'handleSelected');
try delete(handleSelected); end %#ok<TRYNC>

% new grainid
pos = get(ax,'CurrentPoint');
localId = findByLocation(job.grains,[pos(1,1) pos(1,2)]);

if isempty(localId), return; end
grain = job.grains(localId);

% mark grains
hold on
handleSelected = plot(grain.boundary,'lineColor','w','linewidth',4);
hold off
setappdata(ax,'handleSelected',handleSelected);

votesFit = job.calcGBVotes(grain.id,'bestFit','reconsiderAll');
votesProb = job.calcGBVotes(grain.id,'reconsiderAll','numFit',24,'tolerance',5*degree,'curvatureFactor',1);

fig = figure(100);
clf(fig)
set(fig,'name',['grain: ' xnum2str(grain.id)])
numV = size(votesFit.parentId,2);
cKey = ipfHSVKey(job.csParent);

oriPV = variants(job.p2c,job.grainsPrior(localId).meanOrientation,votesFit.parentId);

bgColor = cKey.orientation2color(oriPV);
fgColor = bgColor .* sum(bgColor,2) < 1.5;

for n=1:numV

  handles.b{n} = uicontrol('Style','PushButton','Units','normalized',...
    'Position',[0.1 (n-1)/numV 0.8 0.9*1/numV],...
    'backGroundColor',bgColor(n,:),'ForegroundColor',fgColor(n,:),...
    'String',[xnum2str(votesFit.fit(n)./degree) ' - ' ...
    xnum2str(votesProb.prob(votesProb.parentId==votesFit.parentId(n))) ],...
    'Callback',{@setOri,job,localId,oriPV(n),ax});
end

end

function setOri(a,b,job,id,pOri,ax)

job.grains(id).meanOrientation = pOri;
job.grains.update;

hold(ax,'on')
plot(job.grains(id),pOri,'parent',ax);
hold(ax,'off')

end

author: mtex-toolbox

repo: mtex

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


text:
function evalPCK_multi(expidxs,bUseHeadSize,bHeadSizeFromRect,bTrain)

fprintf('evalPCK()\n');

if (nargin < 2)
    bUseHeadSize = true;
end

if (nargin < 3)
    bHeadSizeFromRect = false;
end

if (nargin < 4)
    bTrain = false;
end

if (bUseHeadSize)
    prefix = 'pck-head-';
    range = 0:0.01:0.5;
else
    prefix = 'pck-torso-';
    range = 0:0.01:0.2;
end

if (bTrain)
    prefix = [prefix 'train-'];
end

if bTrain
    image_set = 'train';
else
    image_set = 'test';
end

legendName = cell(0);

[~,parts] = util_get_parts24();
jidxsUpperBody = 7:12;
nJoints = 14;
% nJoints = length(p.cidxs);
for expidx = expidxs
    fprintf('**********************************************************************************\n');
    % load experiment parameters
    p = exp_params(expidx);

    % load ground truth
    if (~bTrain)
        load(p.testGT);
    else
        load(p.trainGT);
        annolist = annolist(1:17000);
    end

    if (~exist('annolist','var'))
        annolist = single_person_annolist;
    end
    annolist_gt = annolist;%(1:100)

    if (~exist(p.outputDir,'dir')), mkdir(p.outputDir); end
    fnameEvalRes = [fileparts(p.outputDir) '/evalAll'];

    clear boxes boxes_nonms image_ids keypointsAll;

    if (~exist(p.latexDir,'dir')), mkdir(p.latexDir); end
    if (~exist(p.plotsDir,'dir')), mkdir(p.plotsDir); end

%     assert(length(p.pidxs) == nJoints);
%     assert(length(p.cidxs) == nJoints);

    if (~isfield(p,'bidxs'))
        bidxs = 1:length(p.pidxs)+1;
    else
        bidxs = p.bidxs;
    end

    % Load keypoints
    fn = fullfile(p.exp_dir, image_set, 'keypointsAll.mat');
    if exist(fn, 'file') ~= 2
        scoremaps2keypoints(expidx, image_set);
    end
    fprintf('loading %s\n', fn);
    load(fn, 'keypointsAll');

%     visBoxes(annolist_gt,keypointsAll,p.pidxs,parts);
    nrects = 0;
    for imgidx = 1:length(annolist_gt)
        nrects = nrects + length(annolist_gt(imgidx).annorect);
    end

    distAll = nan(nrects,nJoints);
    accAll = zeros(length(range),nJoints+2);

    for i = 1:length(p.pidxs)
      pidx = p.pidxs(i);
      % part is a joint
      assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
      jidx = parts(pidx+1).pos(1);

      % compute distance to the ground truth jidx
      distAll(:,i) = getNormGTJointDistMulti(keypointsAll,annolist_gt,jidx,bUseHeadSize,bHeadSizeFromRect,jidx,nrects);
    end
    matchPCK = double(distAll <= range(end));

    totalPCK = zeros(size(matchPCK,1),1);
    nCorrect = zeros(size(matchPCK,1),1);
    for imgidx = 1:length(totalPCK)
      idxs = ~isnan(distAll(imgidx,:));
%       idxs = ~isnan(matchPCK(imgidx,:));
      totalPCK(imgidx) = sum(matchPCK(imgidx,idxs))/sum(idxs);
%       totalPCK(imgidx) = sum(matchPCK(imgidx,idxs))/length(idxs);
      nCorrect(imgidx) = sum(matchPCK(imgidx,idxs));
    end

    for i = 1:length(p.pidxs)
      dist = distAll(:,i);
      % remove the cases without the ground truth
      dist(isnan(dist)) = [];
      % compute accuracy for each threshold
      for k = 1:numel(range)
        accAll(k,i) = 100*mean(dist<=range(k));
      end
    end

    % compute avg PCKh upper body
    dist = reshape(distAll(:,jidxsUpperBody),size(distAll,1)*length(jidxsUpperBody),1);
    dist(isnan(dist)) = [];
    for k = 1:numel(range)
      accAll(k,end-1) = 100*mean(dist<=range(k));
    end

    % compute avg PCKh full body
    dist = reshape(distAll,size(distAll,1)*size(distAll,2),1);
    dist(isnan(dist)) = [];
    for k = 1:numel(range)
      accAll(k,end) = 100*mean(dist<=range(k));
    end

    pidxs = p.pidxs;
%     save(fnameEvalRes,'distAll','accAll','pidxs','keypointsAll','matchPCK');

    auc = area_under_curve(scale01(range),accAll(:,end));
    legendName{end+1} = sprintf('%s, AUC: %1.1f%%', p.name, auc);

    tableFilename = [p.latexDir '/' prefix 'expidx' num2str(expidx) '.tex'];
    [row, header] = genTableNew(accAll(end,:),auc,p.name);
    row = strrep(row,'NaN','-');
    fid = fopen(tableFilename,'wt');assert(fid ~= -1);
    fprintf(fid,'%s\n',row{1});fclose(fid);
    fid = fopen([p.latexDir '/' prefix 'header.tex'],'wt');assert(fid ~= -1);
    fprintf(fid,'%s\n',header);fclose(fid);

%     fprintf('%s\n',legendName{end});
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% Helper class for GridLayout
classdef CellProps < handle

    properties
        HAlign = 'Stretch'
        VAlign = 'Stretch'
        LMargin = 0
        RMargin = 0
        TMargin = 0
        BMargin = 0
    end
    properties (Hidden, Dependent)
        Margin
    end

    methods (Static)
        function Names = GetCellParamNames()
            Names = cellfun(@(x)(['Cell' x]),CellProps.GetParamNames(), ...
                'UniformOutput',false);
        end
        function Names = GetParamNames()
            Names = { ...
                'HAlign' ...
                'VAlign' ...
                'Margin' ...
                'LMargin' ...
                'RMargin' ...
                'TMargin' ...
                'BMargin'};
        end
    end

    methods
        function Obj = CellProps(varargin)
            if length(varargin) == 1 && isstruct(varargin{1})
                P = varargin{1};
                FieldNames = fieldnames(P);
                for i = 1:length(FieldNames)
                    FieldName = FieldNames{i};
                    if strcmp(FieldName(1:4),'Cell')
                        ObjFieldName = FieldName(5:end);
                    else
                        ObjFieldName = FieldName;
                    end
                    Field = P.(FieldName);
                    if ~isempty(Field)
                        Obj.(ObjFieldName) = Field;
                    end
                end
            else
                for i = 1:2:length(varargin)
                    Obj.(varargin{i}) = varargin{i+1};
                end
            end

        end

        function set.HAlign(Obj, Value)
            Set = {'Left','Center','Right','Stretch'};
            assert(any(strcmp(Value,Set)), ...
                'HAlign must be one of: ''Left'',''Center'',''Right'',''Stretch''.');
            Obj.HAlign = Value;
        end
        function set.VAlign(Obj, Value)
            Set = {'Top','Center','Bottom','Stretch'};
            assert(any(strcmp(Value,Set)), ...
                'VAlign must be one of: ''Top'',''Center'',''Bottom'',''Stretch''.');
            Obj.VAlign = Value;
        end

        function set.Margin(Obj, Value)
            Obj.LMargin = Value;
            Obj.RMargin = Value;
            Obj.TMargin = Value;
            Obj.BMargin = Value;
        end

        function set.LMargin(Obj, Value)
            assert(isscalar(Value) && Value >= 0, ...
                'LMargin must a non-negative scalar.');
            Obj.LMargin = Value;
        end
        function set.RMargin(Obj, Value)
            assert(isscalar(Value) && Value >= 0, ...
                'RMargin must a non-negative scalar.');
            Obj.RMargin = Value;
        end
        function set.TMargin(Obj, Value)
            assert(isscalar(Value) && Value >= 0, ...
                'TMargin must a non-negative scalar.');
            Obj.TMargin = Value;
        end
        function set.BMargin(Obj, Value)
            assert(isscalar(Value) && Value >= 0, ...
                'BMargin must a non-negative scalar.');
            Obj.BMargin = Value;
        end
    end

end

% EOF

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
% forcelocs() - rotate location in 3-D so specified electrodes
%               match specified locations.
%               CAUTION: Only for use on electrodes in
%               and remaining in the upper spherical hemisphere,
%               otherwise it will work improperly. Written primarily for
%               adjusting all electrodes homogenously with Cz.
%
% Usage:
%   >> chanlocs = forcelocs( chanlocs ); % pop-up window mode
%   >> chanlocs = forcelocs( chanlocs, loc1, loc2, ... );
% Example:
%   >> chanlocs = forcelocs( chanlocs, { 0.78, 'x', 'A1' }, { 0.023, 'x', ...
%   'B1','B2','Cz' } );
%
% Inputs:
%   chanlocs  - EEGLAB channel structure. See help readlocs()
%
% Optional inputs:
%   loc1      - cell array: { location, axis, channame1, channame2, .. }
%               'location' is new cartesian coordinate of channame1 along 'axis'
%               'axis' is either
%                 'X'   New x-coordinate of mean of channame1, channame2,
%                       etc. Used to calculate the X-Z plane angle by
%                       which to rotate all channels.
%                       Note that all rotations are to the corresponding positive
%                       Z-value, since theta=atan(z/x).
%                 'Y'   New x-coordinate of mean of channame1, channame2,
%                       etc.
%
%               'channame#'  Name of channel(s) to be rotated, as they appear in
%                       chanlocs.label
%   loc2      - same as loc1
%
% Outputs:
%   chanlocs  - updated EEGLAB channel structure.
%
%
% Author: Arnaud Delorme, CNL / Salk Institute, 15 April 2003
%
% See also: readlocs()

% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [chanlocs,options] = forcelocs( chanlocs, varargin)

    NENTRY = 1; % number of lines in GUI
    FIELDS = { 'X' 'Y' };

    options = [];
    if nargin < 1
        help forcelocs;
        return;
    end;
    if nargin < 2
        geom = { [0.4 1 1 0.3] };
        uilist = { { 'style' 'text' 'string' 'X/Y value' 'tag' 'valstr' } ...
                   { 'style' 'text' 'string' 'Coordinate' } ...
                   { 'style' 'text' 'string' 'Electrode list' } ...
                   { } };
        for index = 1:NENTRY
            tag = [ 'c' int2str(index) ];
            geom = { geom{:}  [0.3 1 1 0.3] };
            uilist = { uilist{:} { 'style' 'edit' 'string' fastif(index==1, '0','') } ...
                       { 'style' 'listbox' 'string' 'X (rotate X-Z plane)|Y (rotate Y-Z plane)' ...
                         'callback' [ 'if get(gco, ''value'') == 1,' ...
                                      '     set(findobj(gcbf, ''tag'', ''valstr''), ''string'', ''Y value'');' ...
                                      'else set(findobj(gcbf, ''tag'', ''valstr''), ''string'', ''X value'');' ...
                                      'end;' ] } ...
                       { 'style' 'edit' 'string'  fastif(index==1, 'Cz','') 'tag' tag } ...
                       { 'style' 'pushbutton' 'string' 'Pick' ...
                         'callback', [ 'tmp3 = get(gcbf, ''userdata'');' ...
                                       '[tmp1 tmp2] = pop_chansel({tmp3.labels}, ''selectionmode'', ''single'');' ...
                                       'if ~isempty(tmp1) set(findobj(gcbf, ''tag'', ''' tag '''), ''string'', tmp2); end;' ...
                                       'clear tmp1 tmp2;' ] } };
        end;

        results = inputgui( geom, uilist, 'pophelp(''forcelocs'');', 'Force electrode location -- forcelocs()', chanlocs );
        if length(results) == 0, return; end;

        options = {};
        for index = 1:NENTRY
            tmpi = 3*(index-1)+1;
            if ~isempty(results{tmpi})
                tmpchans = parsetxt(results{tmpi+2});
                options = { options{:} { str2num(results{tmpi}) FIELDS{results{tmpi+1}} tmpchans{:} }};
            end;
        end;
    else
        options = varargin;
    end;

    % scan all locations
    % ------------------
    channelnames = lower(strvcat({chanlocs.labels}));
    for index = 1:length(options)

        val   = options{index}{1};
        type  = options{index}{2};
        chans = getchans(options{index}(3:end), channelnames);

        % rotate X-Z plane
        % ----------------
        if strcmpi(type, 'x')
            curx   = mean([ chanlocs(chans).X ]);
            curz   = mean([ chanlocs(chans).Z ]);
            newx = val;
            rotangle = solvesystem(curx, curz, newx);

            for chanind = 1:length(chanlocs)
                [chanlocs(chanind).X chanlocs(chanind).Z]= rotation(chanlocs(chanind).X, chanlocs(chanind).Z, rotangle);
            end;
            chanlocs = convertlocs(chanlocs, 'cart2all');
        end;

        % rotate Y-Z plane
        % ----------------
        if strcmpi(type, 'y')
            cury   = mean([ chanlocs(chans).Y ]);
            curz   = mean([ chanlocs(chans).Z ]);
            newy = val;
            rotangle = solvesystem(cury, curz, newy);

            for chanind = 1:length(chanlocs)
                [chanlocs(chanind).Y chanlocs(chanind).Z]= rotation(chanlocs(chanind).Y, chanlocs(chanind).Z, rotangle);
            end;
            chanlocs = convertlocs(chanlocs, 'cart2all');
        end;

    end;

% get channel indices
% -------------------
function chanlist = getchans(chanliststr, channelnames);
    chanlist = [];
    for index = 1:length(chanliststr)
        i = strmatch (lower(chanliststr{index}), channelnames, 'exact');
        chanlist  = [chanlist i];
    end;

% function rotate coordinates
% ---------------------------
function [X,Y] = rotation(x,y,rotangle)
    X = real((x+j*y)*exp(j*rotangle));
    Y = imag((x+j*y)*exp(j*rotangle));

% function solvesyst
% ------------------
function theta = solvesystem(x,y,nx)
    % Original Solution
    %eq(1,:) = [x -y]; res(1) = nx;
    %eq(2,:) = [y x];  res(2) = sqrt(x^2+y^2-nx^2);
    %sol = eq\res';
    %theta = atan2(sol(2), sol(1));

    % simplified solution
    ny = sqrt(x^2+y^2-nx^2);
    ang1 = angle(x+j*y);
    ang2 = angle(nx+j*ny);
    theta = ang2-ang1;

    % Even simpler solution     Toby 03/05/2007
    % theta = atan(y/x);

author: PatternRecognition

repo: OpenBMI


text:
function path = vl_root
% VL_ROOT  Obtain VLFeat root path
%   PATH = VL_ROOT() returns the path to the VLFeat installation.
%
%   See also: VL_SETUP(), VL_HELP().

% Authors: Andrea Vedaldi and Brian Fulkerson

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

[a,b,c] = fileparts(mfilename('fullpath')) ;
[a,b,c] = fileparts(a) ;
path = a ;

author: yihui-he

repo: panorama

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


text:
% The COBRAToolbox: testAnalyzeGCdesign.m
%
% Purpose:
%     - test the analyzeGCdesign function
%
% Authors:
%     - Jacek Wachowiak

% save the current path
currentDir = pwd;

%This tests requires gurobi
requireOneSolverOf = {'tomlab_cplex', 'ibm_cplex', 'gurobi'};
solverPkgs = prepareTest('needsLP',true,'needsQP',true,'requireOneSolverOf', requireOneSolverOf);

% initialize the test
fileDir = fileparts(which('testAnalyzeGCdesign'));
cd(fileDir);

% test variables
model = getDistributedModel('ecoli_core_model.mat');
changeCobraSolver(solverPkgs.LP{1}, 'LP');
fprintf('Reducing model with %s ...\n',solverPkgs.LP{1})

modelRed = reduceModel(model);
selectedRxns = modelRed.rxns(22:25);
target = modelRed.rxns(20); % 'EX_ac(e)'
deletions = modelRed.rxns(21);

% function outputs
% solver change due to instability of qpng
% to be changed with gurobi

for k = 1:numel(solverPkgs.QP)
    changeCobraSolver(solverPkgs.QP{k}, 'QP');
    if any(ismember(solverPkgs.QP{k}, solverPkgs.LP))
        changeCobraSolver(solverPkgs.QP{k}, 'LP');
    end
    fprintf('Testing AnalyzeGCdesign with %s ...\n',solverPkgs.LP{k})
    improvedRxnsM = {};
    intermediateSlnsM = {};
    [improvedRxns, intermediateSlns] = analyzeGCdesign(modelRed, selectedRxns, target, deletions);

    for i=2:8
        [improvedRxns2, intermediateSlns2] = analyzeGCdesign(modelRed, selectedRxns, target, deletions, i, i);
        improvedRxnsM{end+1} = improvedRxns2;
        intermediateSlnsM{end+1} = intermediateSlns2;
    end

    % tests
    assert(isequal({'EX_akg(e)', 'EX_co2(e)'}, improvedRxns));
    assert(isequal({'EX_acald(e)'}, intermediateSlnsM{7}{1}));
end

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function halham_write ( dim_num, n, step, seed, leap, base, r, file_out_name )

%*****************************************************************************80
%
%% HALHAM_WRITE writes a Halton or Hammersley dataset to a file.
%
%  Discussion:
%
%    The initial lines of the file are comments, which begin with a
%    '#' character.
%
%    Thereafter, each line of the file contains the DIM_NUM-dimensional
%    components of the next entry in the dataset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer N, the number of elements in the subsequence.
%
%    Input, integer STEP, the index of the subsequence element.
%    0 <= STEP is required.
%
%    Input, integer SEED(DIM_NUM), the sequence index corresponding
%    to STEP = 0.
%
%    Input, integer LEAP(DIM_NUM), the successive jumps in the sequence.
%
%    Input, integer BASE(DIM_NUM), the bases.
%
%    Input, real R(DIM_NUM,N), the points.
%
%    Input, string FILE_OUT_NAME, the name of
%    the output file.
%
  file_out_unit = fopen ( file_out_name, 'w' );

  if ( file_out_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HALHAM_WRITE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file:\n' );
    fprintf ( 1, '  "%s".\n', file_out_name );
    error ( 'HALHAM_WRITE - Fatal error!' );
  end

  fprintf ( file_out_unit, '#  %s\n', file_out_name );
  fprintf ( file_out_unit, '#  created by HALHAM_WRITE.M\n' );
  fprintf ( file_out_unit, '#\n' );
  fprintf ( file_out_unit, '#  DIM_NUM = %d\n', dim_num );
  fprintf ( file_out_unit, '#  N =    %d\n', n );
  fprintf ( file_out_unit, '#  STEP = %d\n', step );

  for mlo = 1 : 5 : dim_num
    mhi = min ( mlo + 5 - 1, dim_num );
    if ( mlo == 1 )
      fprintf ( file_out_unit, '#  SEED = ' );
    else
      fprintf ( file_out_unit, '#         ' );
    end
    for i = mlo : mhi
      fprintf ( file_out_unit, '  %12d', seed(i) );
    end
    fprintf ( file_out_unit, '\n' );
  end
  for mlo = 1 : 5 : dim_num
    mhi = min ( mlo + 5 - 1, dim_num );
    if ( mlo == 1 )
      fprintf ( file_out_unit, '#  LEAP = ' );
    else
      fprintf ( file_out_unit, '#         ' );
    end
    for i = mlo : mhi
      fprintf ( file_out_unit, '  %12d', leap(i) );
    end
    fprintf ( file_out_unit, '\n' );
  end
  for mlo = 1 : 5 : dim_num
    mhi = min ( mlo + 5 - 1, dim_num );
    if ( mlo == 1 )
      fprintf ( file_out_unit, '#  BASE = ' );
    else
      fprintf ( file_out_unit, '#         ' );
    end
    for i = mlo : mhi
      fprintf ( file_out_unit, '  %12d', base(i) );
    end
    fprintf ( file_out_unit, '\n' );
  end
  fprintf ( file_out_unit, '#\n' );

  for j = 1 : n
    for i = 1 : dim_num
      fprintf ( file_out_unit, '  %10f', r(i,j) );
    end
    fprintf ( file_out_unit, '\n' );
  end

  fclose ( file_out_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x_best, psi_best, out] = GLM_FY(mapp, lin_sym_solver, x0, options)
% GLM_FY is a Levenberg-Marquardt algorithm for solving systems of
% nonlinear equations :math:`h(x) = 0`, `x` in :math:`R^m`
% using the nonlinear unconstrained minimization :math:`\textrm{min} \psi(x) = 1/2 ||h(x)||^2`
% s.t. `x` in :math:`R^m`.
%
% USAGE:
%
%    [x_best, psi_best, out] = GLM_FY(mapp, lin_sym_solver, x0, options)
%
% INPUTS:
%    mapp:              function handle provides `h(x)` and gradient `h(x)`
%    lin_sym_solver:    function handle for solving the linear system
%    x0:                initial point
%    options:           structure including the parameteres of scheme
%
%                         * .eta - parameter of the scheme
%                         * .MaxNumIter - maximum number of iterations
%                         * .MaxNumMapEval - maximum number of function evaluations
%                         * .MaxNumGmapEval - maximum number of subgradient evaluations
%                         * .TimeLimit - maximum running time
%                         * .epsilon - accuracy parameter
%                         * .x_opt - optimizer
%                         * .psi_opt - optimum
%                         * .adaptive - update lambda adaptively
%                         * .flag_x_error - 1: saves :math:`x_{error}`, 0: do not saves :math:`x_{error}` (default)
%                         * .flag_psi_error - 1: saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                         * .flag_time - 1: saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                         * .Stopping_Crit - stopping criterion
%
%                           1. stop if :math:`||grad|| \leq \epsilon`
%                           2. stop if :math:`||nhxk|| \leq \epsilon`
%                           3. stop if `MaxNumIter` is reached
%                           4. stop if `MaxNumMapEval` is reached
%                           5. stop if `MaxNumGmapEval` is reached
%                           6. stop if `TimeLimit` is reached
%                           7. stop if :math:`||grad|| \leq \textrm{max}(\epsilon, \epsilon^2 * ngradx0)`
%                           8. stop if :math:`||nhxk|| \leq \textrm{max}(\epsilon, \epsilon^2 * nhx0)`
%                           9. stop if (default) :math:`||hxk|| \leq \epsilon` or `MaxNumIter` is reached
%
% OUTPUTS:
%    x_best:            the best approximation of the optimizer
%    psi_best:          the best approximation of the optimum
%    out:               structure including more output information
%
%                         * .T - running time
%                         * .Niter - total number of iterations
%                         * .Nmap - total number of mapping evaluations
%                         * .Ngmap - total number of mapping gradient evaluations
%                         * .merit_func - array including all merit function values
%                         * .x_error - relative error :math:`\textrm{norm}(x_k(:)-x_{opt}(:))/\textrm{norm}(x_{opt})`
%                         * .psi_error - relative error :math:`(\psi_k-\psi_{opt})/(\psi_0-\psi_{opt}))`
%                         * .Status - reason of termination
%
% .. REFERENCE:
% .. J. Fan, Y. Yuan, On  the  quadratic  convergence  of  the Levenberg-Marquardt method without nonsingularity assumption, Computing, 74(1), 23-39  (2005)
% .. Author: - Masoud Ahookhosh, System Biochemistry Group, Luxembourg Center for System Biomedicine, University of Luxembourg, Luxembourg
%            - Update: July 2017, M. Ahookhosh

format longG ;

% ================ Error messages for input and output =================
if nargin > 4
    error('The number of input arguments is more than what is needed');
elseif nargin < 4
    error('The number of input arguments is not enough');
end;

if isempty(mapp)
    error('the function handle mapp has to be defined');
elseif ~isa(mapp,'function_handle')
    error('mapp should be a function handle');
end

if isempty(lin_sym_solver)
    error('the function handle lin_sym_solver has to be defined');
elseif ~isa(lin_sym_solver,'function_handle')
    error('lin_sym_solver should be a function handle');
end

if isempty(x0)
    error('The starting point x0 has to be defined');
elseif ~isa(x0,'numeric')
    error('x0 should be a numeric vector');
end

% =================== initializing the parameters ======================
% ===== user has requested viewing the default values of "options" =====
[eta,epsilon,MaxNumIter,MaxNumMapEval,MaxNumGmapEval, adaptive, ...
    TimeLimit,flag_x_error,flag_psi_error,flag_time,Stopping_Crit] ...
    = Initialization(options);

if ~isa(eta,'numeric') || (eta <= 0)
    error('eta should be numeric and eta in (0,4*delta)');
end

if isfield(options,'x_opt')
    x_opt=options.x_opt;
elseif flag_x_error==1
    error('x_error requires to x_opt be specified');
end

if flag_x_error == 1
    Nxopt      = sqrt(sum(x_opt(:).^2));
    x_error(1) = sqrt(sum((x0(:)-x_opt(:)).^2))/Nxopt;
end

if flag_psi_error == 1
    psi_error(1) = 1;
end

if flag_time == 1
    Time(1) = 0;
end

xk         = x0;
Niter      = 1;
[hxk,ghxk] = mapp(x0);
Nmap       = 1;
Ngmap      = 1;
grad       = ghxk*hxk;
nhxk       = sqrt(sum(hxk.^2));
nhx0       = nhxk;
ngradx0    = sqrt(sum(grad.^2));
I          = eye(length(xk));
psik       = 0.5*nhxk^2;
merit_func = psik;
sigma      = 1e-2;
rho        = 0.2;
gamma      = 0.5; % ???
max_inner  = 5;
StopFlag   = 0;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Main body of GLM_FY.m %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T0 = tic;

% ======================= start of the main loop =======================
while ~StopFlag

    muk   = nhxk;
    Hk    = ghxk*ghxk'+muk*I;
    dk    = lin_sym_solver(Hk,grad);
    xk1   = xk+dk;
    hxk1  = mapp(xk1);
    Nmap  = Nmap+1;
    nhxk1 = norm(hxk1);
    if nhxk1 <= gamma*nhxk
        xk    = xk1;
        Niter = Niter+1;
        hxk   = mapp(xk);
        Nmap  = Nmap+1;
        nhxk  = norm(hxk);
        psik  = 0.5*nhxk^2;
    else
        psik1       = 0.5*nhxk1^2;
        alphak      = 1;
        inner_count = 0;
        while ((psik1-psik)>=sigma*alphak*(grad'*dk)&& ...
                inner_count<= max_inner)
            alphak      = rho*alphak;
            xk1         = xk+alphak*dk;
            hxk1        = mapp(xk1);
            Nmap        = Nmap+1;
            nhxk1       = norm(hxk1);
            psik1       = 0.5*nhxk1^2;
            inner_count = inner_count+1;
        end
        xk    = xk1;
        Niter = Niter+1
        nhxk  = nhxk1;
        psik  = psik1;
    end
    [hxk,ghxk] = mapp(xk);
    Nmap       = Nmap+1;
    Ngmap      = Ngmap+1;
    grad       = ghxk*hxk;

    % ================= Gathering output information ===================
    merit_func(Niter) = psik;
    if flag_time == 1
        Time(Niter+1) = toc(T0);
    end

    if flag_x_error == 1
        Nx_opt = norm(x_opt);
        x_error(Niter+1) = sqrt(sum((xk(:)-x_opt(:)).^2))/Nx_opt;
    end

    if flag_psi_error == 1
        psi_error(Niter+1) = (psik-psi_opt)/(psi0-psi_opt);
    end

    % ================== checking stopping criteria ====================
    T = toc(T0);

    [StopFlag,Status] = StopCriterion(grad,nhxk,Niter,Nmap, ...
    Ngmap,MaxNumIter,MaxNumMapEval,MaxNumGmapEval,T,TimeLimit, ...
    epsilon,nhx0,ngradx0,Stopping_Crit);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Status
x_best         = xk;
psi_best       = psik;
out.T          = T;
out.nhx        = nhxk;
out.merit_func = merit_func';
out.Niter      = Niter;
out.Nmap       = Nmap;
out.Ngmap      = Ngmap;
out.Status     = Status;

if flag_x_error == 1
    out.x_error = x_error;
end
if flag_psi_error == 1
    out.psi_error = psi_error;
end
if flag_time == 1
    out.Time = Time;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% End of GLM_FY.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function cp_ptp = cat_cp_ptp(cp, ptp)
% Concatenate CP and PTP trajectory
%
% [Inputs]
%	cp			: CP trajectory
%	ptp			: PTP trajectory
% [Outputs]
%	cp_ptp		: CP + PTP trajectory

num_cp = length(cp);
num_ptp = length(ptp);
cp_ptp = cell(1, num_cp + num_ptp);
for n = 1:num_cp
	cp_ptp{2 * n} = cp{n};
	cp_ptp{2 * n - 1} = ptp{n};
end
cp_ptp{end} = ptp{end};

author: Sable


text:
function [Recognized_index filename] = Recognition(Test_image, Mean, A, eigenfaces)

ProjectedImages = [];
Train_Number = size(eigenfaces,2);
for i = 1 : Train_Number
    temp = eigenfaces'*A(:,i);
    ProjectedImages = [ProjectedImages temp];
end

% Extract the PCA features from test image

temp = Test_image;

[irow icol] = size(temp);
InImage = reshape(temp',irow*icol,1);
Difference = double(InImage)-Mean;
ProjectedTestImage = eigenfaces'*Difference;

% Calculate Euclidean distances
% Test image is supposed to have minimum distance with its corresponding
% image in the training database.

dist = [];
for i = 1 : Train_Number
    q = ProjectedImages(:,i);
    temp = ( norm( ProjectedTestImage - q ) )^2;
    dist = [dist temp];
end

[dist_min , Recognized_index] = min(dist);

 filename=sprintf('s%d.1.tif', Recognized_index);
   disp(['matched image is ',filename]);

subplot(1,2,1); imshow(Test_image);
title('Test Image');
subplot(1,2,2); imshow(filename);
title('Matched Image')

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function LFs = LRFforPntCldFunc(mesh, keypntIdx, neighborSize)

%  Author: Yulan Guo {yulan.guo@nudt.edu.cn}
%  NUDT, China & CSSE, UWA, Australia
%
% This function takes a point cloud as an input, and generate local reference frames (LRFs) for a set of
% keypoints as an output.
%
% Arguments : mesh - with vertices and faces
%                       keypntIdx - the indices of keypoints on a mesh
%                       neighborSize - the size of neighborhood to define a
%                                                  local surface for a selected keypoint
% Return :         LRFs - local reference frames (LRFs) corresponding to all
%                                     keypoints
% Copyright : This code is written by Yulan Guo {yulan.guo@nudt.edu.cn}, NUDT.
%               The code may be used, modified and distributed for research purposes with
%              acknowledgement of the author and inclusion this copyright information.
%References:
%           [1] Yulan Guo, Ferdous Sohel, Mohammed Bennamoun, Min Lu, Jianwei Wan.
%           Rotational Projection Statistics for 3D Local Surface Description and Object Recognition.
%           Internation Journal of Computer Vision. 2013, 105 (1), 63-86
%           [2] Yulan Guo, Mohammed Bennamoun, Ferdous A Sohel, Min Lu, Jianwei Wan.
%           3D Object Recognition in Cluttered Scenes with Local Surface Features: A Survey.
%           IEEE Transactions on Pattern Analysis and Machine Intelligence,2014
%
% Disclaimer : This code is provided as is without any warranty.

%%%%%%%%%%%% parallel computing to accelerate.
p = gcp('nocreate'); % If no pool, do not create new one.
if isempty(p)   % if exists no parallel computing handle, create one.
    poolsize = parpool;
    % disp('Create a Parallel Computing Pool to accelarate process' );
else
    poolsize = p.NumWorkers;
    % disp(sprintf( 'Parpool has %d workers', poolsize) );
end

%%%%%%
% BucketSize = floor(length(mesh.vertices)/100);
% kdtreeVertices = KDTreeSearcher(mesh.vertices,'Distance','euclidean','BucketSize',BucketSize);
kdtreeVertices = KDTreeSearcher(mesh.vertices,'Distance','euclidean');
%%%%%%
parfor i = 1:length(keypntIdx)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    keypnt = mesh.vertices(keypntIdx(i),:);
    [neighborIdx,neighborDis]= rangesearch(kdtreeVertices,keypnt,neighborSize);
    neighborIdx = cell2mat(neighborIdx);
    neighborIdx = neighborIdx(2:end);
    neighborDis = cell2mat(neighborDis);
    neighborDis = neighborDis(2:end);
    %%%%%
    M = zeros(3,3);
    dis = 0;
    for j = 1:length(neighborIdx)
        M = M+(mesh.vertices(neighborIdx(j),:)-mesh.vertices(keypntIdx(i),:))'*(mesh.vertices(neighborIdx(j),:)-mesh.vertices(keypntIdx(i),:))*(neighborSize-neighborDis(j));
        dis = dis+(neighborSize-neighborDis(j));
    end
    M = M/dis;
    if isnan(M(1,1)) ==1
        LFs{i,1} = eye(3,3);
        continue;
    end
    [V,D] = eig(M);
    lamda = [D(1,1),D(2,2),D(3,3)];
    [temp, idxMinLam] = min(lamda);
    [temp, idxMaxLam] = max(lamda);
    xtemp = V(:,idxMaxLam);
    ztemp = V(:,idxMinLam);
    xPlus = 0;
    xMinus = 0;
    zPlus = 0;
    zMinus = 0;
    for j = 1:length(neighborIdx)
        if (mesh.vertices(neighborIdx(j),:)-mesh.vertices(keypntIdx(i),:))*xtemp>0
            xPlus = xPlus+1;
        else
            xMinus = xMinus+1;
        end
        if (mesh.vertices(neighborIdx(j),:)-mesh.vertices(keypntIdx(i),:))*ztemp>0
            zPlus = zPlus+1;
        else
            zMinus = zMinus+1;
        end
    end
    if xPlus>xMinus
        xAxis = xtemp;
    else
        xAxis = -xtemp;
    end
   if zPlus>zMinus
        zAxis = ztemp;
    else
        zAxis = -ztemp;
    end
    yAxis = cross(zAxis,xAxis);
    rotation = [xAxis';yAxis';zAxis'];
    LFs{i,1} = rotation;
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR