text:
%SerialLink.IKINE Inverse manipulator kinematics
%
% Q = R.ikine(T) is the joint coordinates corresponding to the robot
% end-effector pose T which is a homogenenous transform.
%
% Q = R.ikine(T, Q0, OPTIONS) specifies the initial estimate of the joint
% coordinates.
%
% Q = R.ikine(T, Q0, M, OPTIONS) specifies the initial estimate of the joint
% coordinates and a mask matrix.  For the case where the manipulator
% has fewer than 6 DOF the solution space has more dimensions than can
% be spanned by the manipulator joint coordinates.  In this case
% the mask matrix M specifies the Cartesian DOF (in the wrist coordinate
% frame) that will be ignored in reaching a solution.  The mask matrix
% has six elements that correspond to translation in X, Y and Z, and rotation
% about X, Y and Z respectively.  The value should be 0 (for ignore) or 1.
% The number of non-zero elements should equal the number of manipulator DOF.
%
% For example when using a 5 DOF manipulator rotation about the wrist z-axis
% might be unimportant in which case  M = [1 1 1 1 1 0].
%
% In all cases if T is 4x4xM it is taken as a homogeneous transform sequence
% and R.ikine() returns the joint coordinates corresponding to each of the
% transforms in the sequence.  Q is MxN where N is the number of robot joints.
% The initial estimate of Q for each time step is taken as the solution
% from the previous time step.
%
% Options::
%
% Notes::
% - Solution is computed iteratively using the pseudo-inverse of the
%   manipulator Jacobian.
% - The inverse kinematic solution is generally not unique, and
%   depends on the initial guess Q0 (defaults to 0).
% - Such a solution is completely general, though much less efficient
%   than specific inverse kinematic solutions derived symbolically.
% - This approach allows a solution to obtained at a singularity, but
%   the joint angles within the null space are arbitrarily assigned.
%
% See also SerialLink.fkine, tr2delta, SerialLink.jacob0, SerialLink.ikine6s.

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, 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

function qt = ikine(robot, tr, varargin)
    %  set default parameters for solution
    opt.ilimit = 100;
    opt.tol = 1e-6;
    opt.alpha = 1;
    opt.plot = false;
    opt.pinv = false;

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

    n = robot.n;

    % robot.ikine(tr, q)
    if length(args) > 0
        q = args{1};
        if isempty(q)
            q = zeros(1, n);
        else
            q = q(:)';
        end
    else
        q = zeros(1, n);
    end

    % robot.ikine(tr, q, m)
    if length(args) > 1
        m = args{2};
        m = m(:);
        if numel(m) ~= 6
            error('Mask matrix should have 6 elements');
        end
        if numel(find(m)) ~= robot.n
            error('Mask matrix must have same number of 1s as robot DOF')
        end
    else
        if n < 6
            error('For a manipulator with fewer than 6DOF a mask matrix argument must be specified');
        end
        m = ones(6, 1);
    end
    % make this a logical array so we can index with it
    m = logical(m);

    npoints = size(tr,3);    % number of points
    qt = zeros(npoints, n);  % preallocate space for results
    tcount = 0;              % total iteration count
    eprev = Inf;

    save.e = Inf;
    save.q = [];

    history = [];
    for i=1:npoints
        T = tr(:,:,i);
        nm = Inf;
        count = 0;

        optim = optimset('Display', 'iter', 'TolX', 0, 'TolFun', opt.tol, 'MaxIter', opt.ilimit);
        optim
        q = fminsearch(@(x) ikinefunc(x, T, robot, m), q, optim);

    end
    qt = q;
end

function E = ikinefunc(q, T, robot, m)
        e = tr2delta(fkine(robot, q'), T);
        E = norm(e(m));
end

author: petercorke

repo: robotics-toolbox-matlab


text:
classdef RectangularSensorState < AbstractSensorState
    %RectangularSensorState Summary of this class goes here
    %   Detailed explanation goes here

    properties
        sensor RectangularSensor

        activeTf(1,1) logical = true;
        steeringModel(1,1) AbstractSensorSteeringModel
        azAngle(1,1) double {mustBeGreaterThan(azAngle,0)} = deg2rad(10);
        decAngle(1,1) double {mustBeGreaterThan(decAngle,0)} = deg2rad(10);
        range(1,1) double {mustBeGreaterThan(range,0)} = 1;
    end

    methods
        function obj = RectangularSensorState(sensor, activeTf, steeringModel, azAngle, decAngle, range)
            arguments
                sensor(1,1) RectangularSensor
                activeTf(1,1) logical
                steeringModel(1,1) AbstractSensorSteeringModel
                azAngle(1,1) double {mustBeGreaterThanOrEqual(azAngle,0)}
                decAngle(1,1) double {mustBeGreaterThanOrEqual(decAngle,0)}
                range(1,1) double {mustBeGreaterThan(range, 0)}
            end
            obj.sensor = sensor;
            obj.activeTf = activeTf;
            obj.steeringModel = steeringModel;
            obj.azAngle = azAngle;
            obj.decAngle = decAngle;
            obj.range = range;
        end

        function sensor = getSensor(obj)
            sensor = obj.sensor;
        end

        function activeTf = getSensorActiveState(obj)
            activeTf = obj.activeTf;
        end

        function setSensorActiveState(obj, activeTf)
            obj.activeTf = activeTf;
        end

        function steeringModel = getSensorSteeringMode(obj)
            steeringModel = obj.steeringModel;
        end

        function setSteeringModel(obj, steeringModel)
            obj.steeringModel = steeringModel;
        end

        function range = getSensorMaxRange(obj)
            range = obj.range;
        end

        function setSensorMaxRange(obj, range)
            obj.range = range;
        end

        function azAngle = getSensorAzAngle(obj)
            azAngle = obj.azAngle;
        end

        function setSensorAzAngle(obj, azAngle)
            obj.azAngle = azAngle;
        end

        function decAngle = getSensorDecAngle(obj)
            decAngle = obj.decAngle;
        end

        function setSensorDecAngle(obj, decAngle)
            obj.decAngle = decAngle;
        end

        function angle = getMaxAngle(obj)
            angle = max([obj.getSensorAzAngle, obj.getSensorDecAngle()]);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% TTEST_CELL - compute paired t-test. Allow fast computation of
%                multiple t-test using matrix manipulation.
%
% Usage:
%    >> [F df] = ttest_cell( { a b } );
%    >> [F df] = ttest_cell(a, b);
%
% Inputs:
%   a,b       = data consisting of PAIRED arrays to be compared. The last
%               dimension of the data array is used to compute the t-test.
% Outputs:
%   T   - T-value
%   df  - degree of freedom (array)
%
% Example:
%   a = { rand(1,10) rand(1,10)+0.5 }
%   [T df] = ttest_cell(a)
%   signif = 1-tcdf(T, df(1))
%
%   % for comparison, the same using the Matlab t-test function
%   [h p ci stats] = ttest(a{1}', b{1}');
%   [ stats.tstat' p]
%
%   % fast computation (fMRI scanner volume 100x100x100 and 10 subjects in
%   % two conditions). The computation itself takes 0.5 seconds instead of
%   % half an hour using the standard approach (1000000 loops and Matlab
%   % t-test function)
%   a = rand(100,100,100,10); b = rand(100,100,100,10);
%   [F df] = ttest_cell({ a b });
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005
%
% Reference:
%   Schaum's outlines in statistics (3rd edition). 1999. Mc Graw-Hill.

% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [tval, df] = ttest_cell(a,b)

    if nargin < 1
        help ttest_cell;
        return;
    end

    if iscell(a), b = a{2}; a = a{1}; end
    tmpdiff = a-b;
    diff = mymean(tmpdiff,    myndims(a));
    sd   = mystd( tmpdiff,[], myndims(a));
    tval = diff./sd*sqrt(size(a, myndims(a)));
    df   = size(a, myndims(a))-1;

    % check values againg Matlab statistics toolbox
    %[h p ci stats] = ttest(a', b');
    % [ tval stats.tstat' ]

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end
    end;

function res = mymean( data, varargin) % deal with complex numbers
    res = mean( data, varargin{:});
    if ~isreal(data)
        res = abs( res );
    end

function res = mystd( data, varargin) % deal with complex numbers
    if ~isreal(data)
        res = std( abs(data), varargin{:});
    else
        res = sqrt(sum( bsxfun(@minus, data, mean( data, varargin{2})).^2, varargin{2})/(size(data,varargin{2})-1)); % 8 percent speedup
        %res = std( data, varargin{:});
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function removed = plotlycleanup

% cleans up any old Plotly API MATLAB library files and folders

% initialize output
removed = {};

%----REMOVE WRAPPER FILES----%
REMOVEFILES = {'plotly.m'};

%----REMOVE WRAPPER FOLDERS----%
REMOVEFOLDERS = {'fig2plotly_aux'};

%----check for local Plotly instances----%
try
    plotlyScriptDirs = which('plotly.m','-all');

    if isempty(plotlyScriptDirs);
        error('plotly:missingScript',...
            ['\n\nWe were unable to locate plotly.m. Please Add this\n',...
            'script to your MATLAB search path and try again.\n\n']);
    end

catch exception %locating plotly error catch...
    fprintf(['\n\n' exception.identifier exception.message '\n\n']);
    return
end

% plotly toolbox directory
plotlyToolboxDir = fullfile(matlabroot,'toolbox','plotly');

% find the location of all plotly/ directories
dircount = 1;
for d = 1:length(plotlyScriptDirs)
    %parse filepath string at the Plotly directory
    plotlyLoc = strfind(fileparts(plotlyScriptDirs{d}),fullfile('MATLAB-api-master','plotly'));
    plotlyToolboxLoc = strfind(fileparts(plotlyScriptDirs{d}),plotlyToolboxDir);
    if ~isempty(plotlyLoc)
        plotlyDirs{dircount} = fullfile(plotlyScriptDirs{d}(1:plotlyLoc-1),'MATLAB-api-master','plotly');
        dircount = dircount + 1;
    elseif ~isempty(plotlyToolboxLoc)
        plotlyDirs{dircount} = plotlyToolboxDir;
        dircount = dircount + 1;
    end
end

for d = 1:length(plotlyDirs)

    % add plotlydirs to searchpath (will be removed in future once handled by plotlyupdate)
    addpath(genpath(plotlyDirs{d}));

    % delete files from plotly directory
    removefiles = fullfile(plotlyDirs{d}, REMOVEFILES);

    for f = 1:length(removefiles)

        % remove removefiles filepath from searchpath
        rmpath(fileparts(removefiles{f}));

        if exist(removefiles{f},'file')
            delete(removefiles{f});
            % update removed list
            removed = [removed removefiles{f}];
        end

        % add removefiles filepath back to searchpath
        addpath(fileparts(removefiles{f}));

    end

    % remove folders from plotly directory
    removefolders = fullfile(plotlyDirs{d},REMOVEFOLDERS);

    for f = 1:length(removefolders)
        if exist(removefolders{f},'dir')

            %remove folder from path
            rmpath(genpath(removefolders{f}));

            %delete folder/subfolders
            try
                status = rmdir(removefolders{f},'s');

                if (status == 0)
                    error('plotly:deletePlotlyAPI',...
                        ['\n\nShoot! It looks like something went wrong removing the Plotly API ' ...
                        'from the MATLAB toolbox directory \n' ...
                        'Please contact your system admin or post a topic on https://community.plotly.com/c/api/matlab/ for more information. \n\n']);
                end

                % update removed list
                removed = [removed removefolders{f}];

            end
        end
    end
end
end

author: plotly


text:
function ELF(option, up)
%ELF eliminates very low frequencies from resampled respiratory
% signals.
%	            ELF(option, up)
%
%	Inputs:
%		option          the option which has led to this function being used
%       up              universal parameters structure
%

fprintf('\n--- Eliminating VLFs ');
wave_type = option(1:3);

for subj = up.paramSet.subj_list
    %% Skip if this processing has been done previously
    log_int_respSig = 0;             % Has value 1 unless this is a final respiratory signal
    iden_resp_sig_file_ending
    savepath = [up.paths.data_save_folder, num2str(subj), ending];
    log_int_respSig = 1;             % Has value 1 unless this is a final respiratory signal
    iden_resp_sig_file_ending
    initloadpath = [up.paths.data_save_folder, num2str(subj), ending];
    filecontents = whos('-file', initloadpath);
    var_names = extractfield(filecontents, 'name');
    rel_log = zeros(size(var_names));
    for s = 1 : length(var_names)
        if strfind(var_names{s}, [option(1:3), up.paths.filenames.resampler])
            rel_log(s) = 1;
        end
    end
    rel_var_names = var_names(logical(rel_log));
    for rel_var_name_no = 1 : length(rel_var_names)
        eval(['save_name = ''' option(1:3), up.paths.filenames.elim_vlf2, rel_var_names{rel_var_name_no}(4:end) ''';']);
        exist_log = check_exists(savepath, save_name);
        if exist_log
            continue
        end

        %% Load relevant data
        rel_name = rel_var_names{rel_var_name_no};
        loadpath = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.int_respSigs];
        load(loadpath, rel_name);
        eval(['old_data = ' rel_name ';']);

        %% Eliminate VLFs
        data.hpfilt.t = old_data.t;
        try
            data.hpfilt.v = elim_vlfs(old_data, up);
        catch
            % if there aren't enough points to use the filter, simply carry forward the previous data
            data.hpfilt.v = old_data.v;
        end
        data.hpfilt.fs = old_data.fs;
        eval([save_name ' = data.hpfilt;']);

        %% Save processed data
        save_or_append_data
    end
end

end

author: peterhcharlton

repo: RRest


text:
classdef PTKDivideLungsIntoSagittalBins < PTKPlugin
    % PTKDivideLungsIntoSagittalBins. Plugin for dividing the lungs into bins along
    % the left-right axis
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2014.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Sagittal<br>bins'
        ToolTip = 'Divides the lungs into bins along the left-right axis'
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Context = PTKContextSet.Lungs
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            whole_lung_mask = dataset.GetTemplateImage(PTKContext.Lungs);
            orientation = PTKImageOrientation.Sagittal;
            results = PTKDivideVolumeIntoSlices(whole_lung_mask, orientation, reporting);
            results.SetPreferredOrientation(PTKImageOrientation.Coronal);
            results.SetPreferredSliceNumber(round(results.ImageSize/2));
        end

        function results = GenerateImageFromResults(bins_results, image_templates, reporting)
            results = bins_results.BinImage;
            results.SetPreferredOrientation(PTKImageOrientation.Coronal);
            results.SetPreferredSliceNumber(round(results.ImageSize/2));
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [imgOut, L_wa] = YPWardGlobalTMO(img, Ld_max, maxLayers)
%
%       [imgOut, L_wa] = YPWardGlobalTMO(img, Ld_max, maxLayers)
%
%
%       Input:
%           -img: input HDR image
%           -Ld_max: maximum monitor LDR luminance in cd/m^2
%           -L_wa: world adpatation luminance in cd/m^2
%
%       Output
%           -imgOut: tone mapped image
%           -L_wa: luminance adaptation
%
%     Copyright (C) 2020 Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     This is the TumblinTMO with Yee and Pattanaik Luminance Adaptation;
%     see the file tmo/util/YeePattanaikLuminanceAdaptation.m
%

if(~exist('Ld_max', 'var'))
    Ld_max = 100;
end

if(~exist('maxLayers', 'var'))
    maxLayers = 32;
end

if(Ld_max <= 0.0)
    Ld_max = 100;
end

L_wa = YeePattanaikLuminanceAdaptation(img, maxLayers);

imgOut = WardGlobalTMO(img, Ld_max, L_wa);

end

author: banterle


text:
function [keyPresses , elapsedTime] = RobotHardKeyBoard(CommPort)
%function [keyPresses , elapsedTime] = RobotHardKeyBoard(CommPort);
%permits control of robot via numerical key pad on a physical key board
% Numerical keypad presets:
%8-moves Roomba forward 10cm
%5-moves Roomba backward 5cm
%4-turns Roomba left 15 degrees
%7-slightly turns Roomba left 5 degrees
%6-turns Roomba right 15 degrees
%9-slightly turns Roomba right 5 degrees

%Alphabetical key presets:
%y-moves Roomba forward 10cm
%h-moves Roomba backward 5cm
%g-turns Roomba left 15 degrees
%t-slightly turns Roomba left 5 degrees
%j-turns Roomba right 15 degrees
%u-slightly turns Roomba right 5 degrees
%
% q quits
%
% input is a number for the comm port ex.   use 1, for COM1
% if omitted, set to default to 9.
% outputs: keypresses = total number of key pressed  (q doesn't count, but
% all other keys, including inactive ones ex 'a' do
% elapsedTime time in seconds since first keypress. can be 0 if q pressed
% immediately
disp('Wait for robot beep...')
a = instrfind();
if ~isempty(a)
    fclose(a);
    pause(1);
    delete(a);
    pause(1);
end

if nargin==0
    CommPort = 9;
end

ser = RoombaInit(CommPort);

im = imread('KeyBoardControl.png');

%logging variables
keyPresses = [0 0];
startTime = [];
elapsedTime = [];
%-------------------------------------------
%        SWITCH SPEEDS AND ANGLES HERE
%
% bD=distance forward sD=distance backward
% sA=small turn angle bA=big turn angle
% sV=backwardspeed bV=forwardspeed
% sW=small angle speed bW=big angle speed
%
%-------------------------------------------
%
% speeds and angles
sV = .1; sD = 0.05; % vel in m/s angles in deg W is angluar speed but in m/s?
bV = .2; bD = 0.1;
sW = .05; sA = 5;
bW = .1; bA = 15;

% set up hidden key logger window
callstr = ['set(gcbf,''Userdata'',double(get(gcbf,''Currentcharacter''))) ; uiresume '];
fh = figure('keypressfcn',callstr, ...
    'windowstyle','modal',...
    'units', 'normalized',...
    'position',[.1 .1 .8 .8],...
    'Name','GETKEY', ...
    'userdata','timeout') ;
imshow(im);

disp('Starting Keyboard Mode.   Press q to quit.   ')

while(1)
  uiwait ;
  try
      %disp('here');  why does a single keypress appear to trigger two
      %events.   First with empty, then actual character????   JME 5/26
      k = char( get(fh,'Userdata') );

  catch
    k = nan ;
  end

  if ~isempty(k)&& k~='q'
    keyPresses(1) = keyPresses(1) +1;
    %disp(k)
  end
if isempty(startTime)
    startTime = tic;
end
%-------------------------------------
%       SWITCH KEYBOARD KEYS HERE
%-------------------------------------
if ~isempty(k)
switch k
    %2
    case {'5','h'}
        travelDist(ser, sV, -sD);
    case {'4','g'}
        turnAngle(ser, bW, bA);
    case {'6','j'}
        turnAngle(ser, bW, -bA);
    case {'7','t'}
        turnAngle(ser, sW, sA);
    case {'8','y'}
        travelDist(ser, bV, bD);
    case {'9','u'}
        turnAngle(ser, sW, -sA);
    case '0'
        elapsedTime = toc(startTime);
        fprintf('Sec:  %.2f,  Keystrokes: %d, NegKey: %d \n ', [elapsedTime keyPresses(1) keyPresses(2)]);
    case 'q'
        if ~isempty(startTime)
        elapsedTime = toc(startTime);
        else
            elapsedTime =0;
        end
        fprintf('Sec:  %.2f,  Keystrokes: %d, NegKey: %d \n ', [elapsedTime keyPresses(1) keyPresses(2)]);
        delete(fh);
        break
    otherwise
        disp('Otherwise')
        keyPresses(2) = keyPresses(2) +1;
end % end switch k
end % if ~isempty(k)
end % end while(1)


text:
function S=PartLabelCompSize(S,varargin)

%% PartLabel Assigns particle labels to particles
%% based on features in the carbon edge spectrum. Note that DiffMaps.m must
%% be run on S prior to execution of this function.
%% OUTPUT NOTES:
%% S.CompSize is a matrix of areas for each component with particles in the
%% rows and components in the columns. the components are ordered as
%% follows: ['OCArea','InArea','Sp2Area','KArea','CO3Area','Total Area']

%% 081201 Ryan C. Moffet

if isempty(S.DiffMap)
    return
end
figsav=0;
particle=S.particle;
if ~isempty(varargin)
    rootdir=varargin{1};
    sample=varargin{2};
    figsav=1;
end
ColorVec=[0,170,0;0,255,255;255,0,0;255,170,0;255,255,255]; %% rgb colors of the different components

MatSiz=size(S.DiffMap(:,:,1));
RgbMat=zeros([MatSiz,3]);
RedMat=zeros(MatSiz);
GreMat=zeros(MatSiz);
BluMat=zeros(MatSiz);

cnt=1;

for i=[5,1,3,2,4]; %% cooh, In, sp2, K, co3
        GrayImage=mat2gray(S.DiffMap(:,:,i)); %% Turn into a greyscale with vals [0 1]
        GrayImage=imadjust(GrayImage,[0 1],[0 1],0.1); %% increase contrast
%         GrayImage=medfilt2(GrayImage);  %% something weird was happening
                                          % %to the pixels on the corner
        Thresh=graythresh(GrayImage); %% Otsu thresholding
        Mask=im2bw(GrayImage,Thresh); %% Give binary image
        BinCompMap{cnt}=bwlabel(Mask,8);
        cnt=cnt+1;
end

% [j,k]=find(BinCompMap{3}>0); %% find index of >0 components
% [l,m]=find(S.LabelMat>0);
% S.BinCompMap=BinCompMap;
% if ~isempty(j) || ~isempty(k)
%     [orgspec]=ComponentSpec(S,'COOH');
%     [inorgspec]=ComponentSpec(S,'Inorg');
%     [ecspec]=ComponentSpec(S,'sp2');
%     ecspec(:,2)=ecspec(:,2)-mean(ecspec(1:8,2));
%     svdM = stackSVD(S,ecspec,orgspec,inorgspec);
%     svdM(svdM<0)=0;
%     figure,imagesc(svdM(:,:,1)),colormap gray, colorbar
%     BinCompMap{3}=threshimage(svdM(:,:,1));
% end
%
%% This first loop creates masks for the individual components over the
%% entire field of view. Each component is then defined as a colored
%% component for visualization.
cnt=1;
for i=[5,1,3,2,4]%i=[5,1,2,3] %% loop over chemical components
%         GrayImage=mat2gray(S.DiffMap(:,:,i)); %% Turn into a greyscale with vals [0 1]
%         GrayImage=imadjust(GrayImage,[0 1],[0 1],0.1); %% increase contrast
% %         GrayImage=medfilt2(GrayImage);  %% something weird was happening
%                                           % %to the pixels on the corner
%         Thresh=graythresh(GrayImage); %% Otsu thresholding
%         Mask=im2bw(GrayImage,Thresh); %% Give binary image
%         BinCompMap{cnt}=bwlabel(Mask,8);
%     end
    [j,k]=find(BinCompMap{cnt}>0); %% find index of >0 components
    [l,m]=find(S.LabelMat>0);
    if ~isempty(j) || ~isempty(k)  %% this conditional defines a color for the component areas (if it exists)
        linidx=sub2ind(size(BinCompMap{cnt}),j,k); %% change to linear in dex
        labidx=sub2ind(size(BinCompMap{cnt}),l,m);
        rejidx=setdiff(linidx,labidx);
        BinCompMap{cnt}(rejidx)=0;
        linidx=sub2ind(size(BinCompMap{cnt}),find(BinCompMap{cnt}>0));
        RedMat(linidx)=ColorVec(cnt,1); %% define the color for the ith component
        GreMat(linidx)=ColorVec(cnt,2);
        BluMat(linidx)=ColorVec(cnt,3);
        trmat=zeros(size(RedMat));
        tgmat=zeros(size(RedMat));
        tbmat=zeros(size(RedMat));
        trmat(linidx)=ColorVec(cnt,1);
        tgmat(linidx)=ColorVec(cnt,2);
        tbmat(linidx)=ColorVec(cnt,3);
        ccmap{cnt}(:,:,1)=trmat;
        ccmap{cnt}(:,:,2)=tgmat;
        ccmap{cnt}(:,:,3)=tbmat;
        clear trmat tgmat tbmat
    else
        ccmap{cnt}=zeros(MatSiz(1),MatSiz(2),3);
    end
    cnt=cnt+1;  %% this counter keeps track of the indivdual components.
    clear j k l m linidx GrayImage Thresh Mask rejidx;
end

%% This second loop assigns labels over individual particles defined
%% previously in Diffmaps.m
NumPart=max(max(S.LabelMat));
% LabelStr={'OC','In','K','EC'};
LabelStr={'OC','In','EC','K','CO3'};

CompSize=zeros(NumPart,length(LabelStr)+1);
PartLabel={};
for i=1:NumPart  %% Loop over particles defined in Diffmaps.m
    PartLabel{i}='';
    for j=1:5  %% Loop over chemical components
        [a1,b1]=find(S.LabelMat==i);  %% get particle i
        [a2,b2]=find(BinCompMap{j}>0); %% get component j
        if ~isempty([a1,b1]) && ~isempty([a2,b2])
            linidx1=sub2ind(size(S.LabelMat),a1,b1); %% Linear index for particle
            linidx2=sub2ind(size(BinCompMap{j}),a2,b2); %% Linear index for component
            IdxCom=intersect(linidx1,linidx2); %% find common indices
            if length(IdxCom)>3%0.05*length(linidx1) %% if component makes up greater than 2% of the pixels of the particle...
                    PartLabel{i}=strcat(PartLabel{i},LabelStr{j}); %% give label of component.
                    CompSize(i,j)=length(IdxCom); %% number of pixels in the component.
            end
        end
    end
    if isempty(PartLabel{i})
        PartLabel{i}='NoID'; %% Particles identified by Otsu's mehod here but not in Particle map script
    end
    CompSize(i,j+1)=length(linidx1);  %% number of pixels in the particle
    clear linidx1 linidx2 IdxCom a1 b1 a2 b2;
end
if isempty(PartLabel)
    S.PartLabel='NoID';
else
    S.PartLabel=PartLabel;
end
S=ParticleSize(S);

MatSiz=size(S.LabelMat);
XSiz=S.Xvalue/MatSiz(1);
YSiz=S.Yvalue/MatSiz(2);
CompSize=CompSize.*(XSiz*YSiz); %% Area of components in um^2
S.CompSize=CompSize;

RgbMat(:,:,1)=RedMat;
RgbMat(:,:,2)=GreMat;
RgbMat(:,:,3)=BluMat;
S.RGBCompMap=RgbMat;
for i=1:length(BinCompMap)
    temp{i}=BinCompMap{i};
    temp{i}(temp{i}>1)=1;
end
xdat=[0:XSiz:S.Xvalue];
ydat=[0:YSiz:S.Yvalue];

% Show DiffMaps
figure('Name',particle,'NumberTitle','off','Position',[1,1,634,872])
subplot(3,2,1)
% imagesc(temp{1})
image(xdat,ydat,uint8(ccmap{1}))
title('Carbox')
axis image
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,2)
% imagesc(temp{2})
image(xdat,ydat,uint8(ccmap{2}))
title('Inorg.')
axis image
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,3)
% imagesc(temp{3})
image(xdat,ydat,uint8(ccmap{3}))
title('%sp2>35')
axis image
xlabel('X (\mum)');
ylabel('Y (\mum)');

subplot(3,2,4)
% imagesc(temp{4})
image(xdat,ydat,uint8(ccmap{4}))
title('Potassium')
axis image
xlabel('X (\mum)');
ylabel('Y (\mum)');

% LabelMat=raw2mask(S);
subplot(3,2,5)
image(xdat,ydat,uint8(ccmap{5}))
title('CO_{3}')
axis image
xlabel('X (\mum)');
ylabel('Y (\mum)');

% title('Particle Map')
% axis image

subplot(3,2,6)
image(xdat,ydat,uint8(RgbMat))
title('Mixed')
axis image
xlabel('X (\mum)');
ylabel('Y (\mum)');

if figsav==1
filename=sprintf('%s%s%s%s',rootdir,sample,particle,'_f2_thresh');
saveas(gcf,filename,'tiff');
end
S.BinCompMap=temp;
S=MultPartAvSpec(S);
if figsav==1
filename=sprintf('%s%s%s%s',rootdir,sample,particle,'_f3_spec');
saveas(gcf,filename,'tiff');
close all
end

author: Sable

repo: mcbench-benchmarks


text:
function lambda = householder_eigenvalues ( n, x )

%*****************************************************************************80
%
%% HOUSEHOLDER_EIGENVALUES returns the eigenvalues of a HOUSEHOLDER matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(N,1), the vector that defines the
%    Householder matrix.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  x = x(:);

  lambda(1,1)   = -1.0;
  lambda(2:n,1) = +1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
x = load('ex4x.dat');
y = load('ex4y.dat');
[m, n] = size(x);
x = [ones(m, 1), x];
[m, n] = size(x);

figure;
hold on;
plot(x(y==1, 2), x(y==1, 3), '+')
plot(x(y!=1, 2), x(y!=1, 3), 'o')
xlabel('Exam 1 score')
ylabel('Exam 2 score')

g = inline('1.0 ./ (1.0 + exp(-z))')

theta = zeros(n, 1);
diff = ones(n, 1);
Js = []

while norm(diff) > 1e-6
  h = g(x * theta);
  Js(length(Js)+1) = (1/m) * sum(-y.*log(h) - (1-y).*log(1 - h));
  grad = (1/m) * x' * (h - y);
  H = (1/m) * (repmat(h .* (1 - h), 1, n) .* x)' * x;
  diff = H \ grad;
  theta = theta - diff;
end

iterations = length(Js)
theta
prob = 1 - g([1, 20, 80] * theta)

plot_x = [min(x(:,2))-2, max(x(:,2))+2];
plot_y = (-1./theta(3)).*(theta(2).*plot_x + theta(1));
plot(plot_x, plot_y);
legend('Admitted', 'Not admitted', 'Decision boundary')
hold off;

figure;
plot(0:iterations-1, Js, 'o--', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
ylabel('J(\theta)');
xlabel('Iteration')

author: zellyn


text:
function [d,a] =  training(a,d)

% [results,algorithm] =  training(algorithm,data,loss)

  disp(['normalizing with ' get_name(a) '.... '])

  x=get_x(d); l=get_dim(d);

  if a.scale_type==4
     y=get_y(d);
     corr = zeros(size(x,2),1);
     corr = (mean(x(y==1,:))-mean(x(y==-1,:))).^2;
     st   = std(x(y==1,:)).^2+std(x(y==-1,:)).^2;
     a.corr = corr ./ st;
  end

 if a.scale_type==1 | a.scale_type==3  | a.scale_type==4
  a.mean_vec  = mean(x);
  a.scale_vec = std(x);
  if a.scale_type==3 %% try to do both scalings
    x = x - ones(l,1)*a.mean_vec; x = x * diag(1./a.scale_vec);
  end
 end
 if a.scale_type==2 | a.scale_type==3
   a.mean_vec2  = mean(x')';
   a.scale_vec2 = std(x')';
 end

 d=test(a,d);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function terms = strsplit(s, delimiter)
%STRSPLIT Splits a string into multiple terms
%
%   terms = strsplit(s)
%       splits the string s into multiple terms that are separated by
%       white spaces (white spaces also include tab and newline).
%
%       The extracted terms are returned in form of a cell array of
%       strings.
%
%   terms = strsplit(s, delimiter)
%       splits the string s into multiple terms that are separated by
%       the specified delimiter.
%
%   Remarks
%   -------
%       - Note that the spaces surrounding the delimiter are considered
%         part of the delimiter, and thus removed from the extracted
%         terms.
%
%       - If there are two consecutive non-whitespace delimiters, it is
%         regarded that there is an empty-string term between them.
%
%   Examples
%   --------
%       % extract the words delimited by white spaces
%       ts = strsplit('I am using MATLAB');
%       ts <- {'I', 'am', 'using', 'MATLAB'}
%
%       % split operands delimited by '+'
%       ts = strsplit('1+2+3+4', '+');
%       ts <- {'1', '2', '3', '4'}
%
%       % It still works if there are spaces surrounding the delimiter
%       ts = strsplit('1 + 2 + 3 + 4', '+');
%       ts <- {'1', '2', '3', '4'}
%
%       % Consecutive delimiters results in empty terms
%       ts = strsplit('C,Java, C++ ,, Python, MATLAB', ',');
%       ts <- {'C', 'Java', 'C++', '', 'Python', 'MATLAB'}
%
%       % When no delimiter is presented, the entire string is considered
%       % as a single term
%       ts = strsplit('YouAndMe');
%       ts <- {'YouAndMe'}
%

%   History
%   -------
%       - Created by Dahua Lin, on Oct 9, 2008
%

%% parse and verify input arguments

assert(ischar(s) && ndims(s) == 2 && size(s,1) <= 1, ...
    'strsplit:invalidarg', ...
    'The first input argument should be a char string.');

if nargin < 2
    by_space = true;
else
    d = delimiter;
    assert(ischar(d) && ndims(d) == 2 && size(d,1) == 1 && ~isempty(d), ...
        'strsplit:invalidarg', ...
        'The delimiter should be a non-empty char string.');

    d = strtrim(d);
    by_space = isempty(d);
end

%% main

s = strtrim(s);

if by_space
    w = isspace(s);
    if any(w)
        % decide the positions of terms
        dw = diff(w);
        sp = [1, find(dw == -1) + 1];     % start positions of terms
        ep = [find(dw == 1), length(s)];  % end positions of terms

        % extract the terms
        nt = numel(sp);
        terms = cell(1, nt);
        for i = 1 : nt
            terms{i} = s(sp(i):ep(i));
        end
    else
        terms = {s};
    end

else
    p = strfind(s, d);
    if ~isempty(p)
        % extract the terms
        nt = numel(p) + 1;
        terms = cell(1, nt);
        sp = 1;
        dl = length(delimiter);
        for i = 1 : nt-1
            terms{i} = strtrim(s(sp:p(i)-1));
            sp = p(i) + dl;
        end
        terms{nt} = strtrim(s(sp:end));
    else
        terms = {s};
    end
end

author: PatternRecognition

repo: OpenBMI


text:
function plotTestSuiteResults(testResultsFolder,reconVersion)

% This function prints and summarizes the results of testAllReconstructionFunctions.

%% plot the computed reconstruction features
currentDir = pwd;
cd(testResultsFolder)

%% load all test results
fields = {
    'Carbon_sources_TruePositives'
    'Carbon_sources_FalseNegatives'
    'Fermentation_products_TruePositives'
    'Fermentation_products_FalseNegatives'
    'growsOnDefinedMedium'
    'Metabolite_uptake_TruePositives'
    'Metabolite_uptake_FalseNegatives'
    'Secretion_products_TruePositives'
    'Secretion_products_FalseNegatives'
    'Bile_acid_biosynthesis_TruePositives'
    'Bile_acid_biosynthesis_FalseNegatives'
    'Drug_metabolism_TruePositives'
    'Drug_metabolism_FalseNegatives'
    'PutrefactionPathways_TruePositives'
    'PutrefactionPathways_FalseNegatives'
    'AromaticAminoAcidDegradation_TruePositives'
    'AromaticAminoAcidDegradation_FalseNegatives'
    'Mass_imbalanced'
    'Charge_imbalanced'
    'Leaking_metabolites'
    };

Results=struct;

for i=1:length(fields)
    if isfile([testResultsFolder filesep fields{i} '_' reconVersion '.txt'])
        savedResults = readtable([testResultsFolder filesep fields{i} '_' reconVersion '.txt'], 'ReadVariableNames', false, 'Delimiter','tab', 'format', 'auto');
        Results.(fields{i}) = table2cell(savedResults);
        numberRecons=size(Results.(fields{i}),1);
    else
        Results.(fields{i}) = {};
    end
end

for j=1:length(fields)
    data=Results.(fields{j});
    if size(data,2)>1
        plotdata=[];
        % if there are no entries for data
        if size(data,2)==1
            for k=1:length(data)
                plotdata(k,1)=0;
            end
            label='Number of data points';
        else
            if strcmp(fields{j},'growsOnDefinedMedium')
                plotdata=data(:,2);
                plotdata(find(strcmp(plotdata(:,1),'NA')),:)=[];
                plotdata=str2double(plotdata);
            else
                if isnumeric(data{1,2}) && ~isempty(data{1,2})
                    % if the data is fluxes-plot the values
                    plotdata=cell2mat(data(:,2));
                    if ~any(strcmp(fields{j},{'Number_genes', 'Number_reactions', 'Number_metabolites'}))
                        label='Flux (mmol*gDW-1*hr-1)';
                    else
                        label='Number of data points';
                    end
                else
                    % count the non-empty data entries
                    for k=1:size(data,1)
                        plotdata(k,1)=length(find(~cellfun(@isempty,data(k,2:end))));
                    end
                    label='Number of data points';
                end
            end
        end
        figure;
        hold on
        % need a workaround if all data is zero
        if sum(plotdata)==0
            plotdata(1,1)=0.00001;
        end
        % does not work if all values are equal or there are too few values
        try
            violinplot(plotdata, {reconVersion});
            ylabel(label)
            if contains(pwd,'_refined')
                h=title([strrep(fields{j},'_',' ') ', refined reconstructions']);
            elseif contains(pwd,'_draft')
                h=title([strrep(fields{j},'_',' ') ', draft reconstructions']);
            end
            ylim([0 max(plotdata)+1])
            set(h,'interpreter','none')
            set(gca,'TickLabelInterpreter','none')
            set(gca, 'FontSize', 12)
            print(fields{j},'-dpng','-r300')
        end
    end
end

% calculate specificity and sensitivity
features=fields;
features(cellfun(@isempty,regexp(features,strjoin({'True','False'},'|'))),:)=[];
features=strrep(features,'_FalseNegatives','');
features=strrep(features,'_FalsePositives','');
features=strrep(features,'_TruePositives','');
features=strrep(features,'_TrueNegatives','');
features=unique(features);
Table{1,1}='Feature';
Table{1,2}=strcat('Sensitivity ',reconVersion);
cnt=2;
plotdata=[];
for i=1:length(features)
    TPs=zeros(numberRecons,1);
    FNs=zeros(numberRecons,1);
    Table{cnt,1}=features{i};
    if size(Results.(strcat(features{i},'_TruePositives')),2)==1
        for j=1:size(Results.(strcat(features{i},'_TruePositives')),1)
            TPs(j,1)=0;
        end
    else
        for j=1:size(Results.(strcat(features{i},'_TruePositives')),1)
            TPs(j,1)=length(find(~cellfun(@isempty,Results.(strcat(features{i},'_TruePositives'))(j,2:end))));
        end
    end
    if size(Results.(strcat(features{i},'_FalseNegatives')),2)==1
        for j=1:size(Results.(strcat(features{i},'_FalseNegatives')),1)
            FNs(j,1)=0;
        end
    else
        for j=1:size(Results.(strcat(features{i},'_FalseNegatives')),1)
            FNs(j,1)=length(find(~cellfun(@isempty,Results.(strcat(features{i},'_FalseNegatives'))(j,2:end))));
        end
    end

    Sensitivity=sum(TPs)/(sum(TPs)+sum(FNs));
    if isnan(Sensitivity)
        Table{cnt,2}='N/A';
    else
        Table{cnt,2}=Sensitivity;
    end
    % to plot all results
    plotdata(i,1)=sum(FNs);
    plotdata(i,2)=sum(TPs);
    labels{i,1}=strrep(features{i},'_',' ');
    labels{i,1}=strrep(labels{i,1},'AromaticAminoAcidDegradation','Amino acid degradation');
    labels{i,1}=strrep(labels{i,1},'PutrefactionPathways','Putrefaction pathways');
    cnt=cnt+1;
end
% remove drug metabolism-better plot separately
[C,IA]=intersect(labels,{'Drug metabolism'});
labels(IA)=[];
plotdata(IA,:)=[];

% remove entries with no experimental data
delInd=find(sum(plotdata,2)==0);
labels(delInd)=[];
plotdata(delInd,:)=[];

% worksround if there is only one row in plotdata
if size(plotdata,1) == 1
    plotdata(2,:) = NaN;
end

figure;
hold on
h=bar(plotdata);
barvalues(h)
h(1).FaceColor = [1 0 0];
h(2).FaceColor = [0 0 1];
set(findall(gcf,'-property','FontSize'),'FontSize',12)
xticklabels(labels);
set(gca,'XTick',1:numel(plotdata))
xtickangle(45)
set(gca,'YTickLabel',[])
ylabel('Total number of model predictions')
legend('Number of false negatives','Number of true positives')
set(gca,'TickLabelInterpreter','none')
set(gca,'FontSize',10)
if contains(pwd,'_refined')
    title('Features succesfully and unsuccessfully captured by refined reconstructions')
elseif contains(pwd,'_draft')
    title('Features succesfully and unsuccessfully captured by draft reconstructions')
end
print('FN_vs_TPs','-dpng','-r300')

close all

Table=cell2table(Table);
writetable(Table,[reconVersion '_Properties'],'FileType','spreadsheet','WriteVariableNames',false);

%% summarize more features
% report all unbalanced reactions
if size(Results.Mass_imbalanced,2)>1 || size(Results.Charge_imbalanced,2)>1
    Mass_imbalanced=Results.Mass_imbalanced(:,2:end);
%     Mass_imbalanced(cellfun(@isempty,Mass_imbalanced)==1)=[];
    Charge_imbalanced=Results.Charge_imbalanced(:,2:end);
%     Charge_imbalanced(cellfun(@isempty,Charge_imbalanced)==1)=[];
    Unbalanced_reactions=unique([Mass_imbalanced,Charge_imbalanced]);
    Unbalanced_reactions=cell2table(Unbalanced_reactions);
    writetable(Unbalanced_reactions,'Unbalanced_reactions','FileType','text','WriteVariableNames',false);
end
% report all leaking metabolites
if size(Results.Leaking_metabolites,2)>1
    Leaking_metabolites=Results.Leaking_metabolites(:,2:end);
    Leaking_metabolites(cellfun(@isempty,Leaking_metabolites)==1)=[];
    Leaking_metabolites=cell2table(Leaking_metabolites);
    writetable(Leaking_metabolites,'Leaking_metabolites','FileType','text','WriteVariableNames',false);
end

%% print percentage of reconstructions agreeing with data
Percentages={'Feature','StrainsWithData','PercentAgreeing'};
features={
    'Carbon_sources'
    'Fermentation_products'
    'PutrefactionPathways'
    'Secretion_products'
    'Metabolite_uptake'
    'AromaticAminoAcidDegradation'
    'Bile_acid_biosynthesis'
    'Drug_metabolism'
    'growsOnDefinedMedium'
    };
for i=1:length(features)
    Percentages{i+1,1} = features{i};
    cntData=0;
    cntAgreeing=0;
    if i < length(features)
        for j=1:size(Results.(strcat(features{i},'_TruePositives')),1)
            if size(Results.(strcat(features{i},'_TruePositives')),2) > 1
                if ~isempty(Results.(strcat(features{i},'_TruePositives')){j,2})
                    cntData=cntData+1;
                else
                    if size(Results.(strcat(features{i},'_FalseNegatives')),2) > 1
                        if ~isempty(Results.(strcat(features{i},'_FalseNegatives')){j,2})
                            cntData=cntData+1;
                        end
                    end
                end
            else
                if size(Results.(strcat(features{i},'_FalseNegatives')),2) > 1
                    if ~isempty(Results.(strcat(features{i},'_FalseNegatives')){j,2})
                        cntData=cntData+1;
                    end
                end
            end
            TP=length(find(~cellfun(@isempty,Results.(strcat(features{i},'_TruePositives'))(j,2:end))));
            if size(Results.(strcat(features{i},'_FalseNegatives')),2) > 1
                FN=length(find(~cellfun(@isempty,Results.(strcat(features{i},'_FalseNegatives'))(j,2:end))));
            else
                FN = 0;
            end
            if TP > 0 && FN == 0
                cntAgreeing=cntAgreeing+1;
            end
        end
    else
        growth=length(find(str2double(Results.(features{i})(:,2))==1));
        nogrowth=length(find(str2double(Results.(features{i})(:,2))==0));
        cntData=growth+nogrowth;
        cntAgreeing=growth;
    end
    Percentages{i+1,2} = cntData;
    Percentages{i+1,3} = cntAgreeing/cntData;
end
Percentages=cell2table(Percentages);
writetable(Percentages,[reconVersion '_PercentagesAgreement'],'FileType','spreadsheet','WriteVariableNames',false);

cd(currentDir)
end

author: opencobra

repo: cobratoolbox


text:
% function F = det_F_algebraic(x1,x2,L_COST,NORMALIZE)
% Determines the F by iteratively minimizing the Algebraic error
% Algorithm 11.2 in Hartley & Zisserman, Multiple View Geometry in Computer Vision
% Inputs:
%           x           3xN coordinates of matched points in image 1(homogeneous)
%           xp          3xN coordinates of matched points in image 2(homogeneous)
%           L_COST      1x1 (optional) controls penalization scheme
%                           for the cost function: L_COST 1 leads to L_1 minimization of
%                           the norm of the algebraic error (the one mentioned in the book)
%           NORMALIZE   1x1 (optional) determines if the algorithm should use the normalized points
% Outputs:
%           F           3x3 the fundamental matrix
%
% Author: Omid Aghazadeh, KTH(Royal Institute of Technology), 2010/05/09
function F = det_F_algebraic(x1,x2,L_COST,NORMALIZE)
if sum(size(x1)~=size(x2)), error('size of correspondences do not match!'), end
if size(x1,1) ~= 3, error('invalid points'), end
if nargin<3; L_COST = 1; end;
if nargin<4; NORMALIZE = 1; end
global A F_I EPSILON_I L_e TOL_X TOL_FUN MAX_FUN_EVAL MAX_ITER;
if NORMALIZE
    nmat1 = get_normalization_matrix(x1); % isotropic normalization (translation/scaling)
    nmat2 = get_normalization_matrix(x2);
    x1n = nmat1*x1;
    x2n = nmat2*x2;
    x1n = x1n./repmat(x1n(3,:),3,1); x2n = x2n./repmat(x2n(3,:),3,1); % normalizing points so their last coordinate is 1
else
    nmat1 = eye(3); nmat2 = eye(3); x1n = x1; x2n = x2;
end
F_0 = det_F_normalized_8point(x1n,x2n);

F_0 = F_0/norm(F_0(:));
F_0T = reshape(F_0',9,1); % f corresponds to F_0 in row order
F_I = [F_0T];
A = [(repmat(x2n(1,:),3,1).*x1n)',(repmat(x2n(2,:),3,1).*x1n)',x1n']; % the data matrix

EPSILON_I = norm(A*F_0T)^L_COST; %initial cost

%% (i) finding an initial estimate of the epipole from an initial estimate of F
[e_0,eprime_0] = get_epipole(F_0);
it=0;
L_e = e_0;

e_i = lsqnonlin(@(e)costAlgebraic(A,L_COST,e),e_0,[],[],optimset('Display','off','TolX',TOL_X,'TolFun',TOL_FUN,'MaxFunEval',MAX_FUN_EVAL,'MaxIter',MAX_ITER,'Algorithm',{'levenberg-marquardt' 0.01}));

F = reshape(F_I,3,3)';
F = nmat2'*F*nmat1;
end

%% function x = const_min_subject_span_space(A,G,r)
% finds the vector x that maximizes ||A x|| subject to the constraint ||x||=1 and x=G x_hat
% where G has rank r
function x = const_min_subject_span_space(A,G,r)
[u_g,s_g,v_g] = svd(G);
u_prime = u_g(:,1:r);

[u_a,s_a,v_a] = svd(A*u_prime);
x_p = v_a(:,end);

x = u_prime * x_p;
end

%% epsilon_cost_i = costAlgebraic(A,L_COST,e_i)
% the cost function for LM algorithm. The goal is to minimize || epsilon_i || by varying e_i
function epsilon_cost_i = costAlgebraic(A,L_COST,e_i)
global F_I EPSILON_I L_e;
%% (ii) find f_i that minimizes ||A f_i|| subject to ||f_i|| = 1
e_cross = get_x_cross(e_i);
E = [e_cross zeros(3) zeros(3); zeros(3) e_cross zeros(3); zeros(3) zeros(3) e_cross];
f_i = const_min_subject_span_space(A,E,6);

%% (iii) compute epsilon_i = A f_i and correct its sign
pos_sign = double(sign(e_i'*L_e));
f_i = f_i * pos_sign;

epsilon_i = A*f_i;
epsilon_cost_i = (norm(epsilon_i))^L_COST;

if epsilon_cost_i < EPSILON_I
    EPSILON_I = epsilon_cost_i;
    F_I = f_i;
    L_e = e_i;
end
end


text:
function gpcf = gpcf_additive(varargin)
%GPCF_ADDITIVE  Create a mixture over products of kernels for each dimension
%
%  Description
%    GPCF = GPCF_ADDITIVE('PARAM1',VALUE1, 'PARAM2,VALUE2, ...)
%    creates a mixture over all possible product combinations of given
%    covariance functions for each input dimension separately.
%
%    Parameters [default]:
%      cf           - cell array {CF_1, CF_2, ... , CF_N} of covariance
%                     functions for each dimension [no default]
%      max_deg      - maximum order of interaction (must be less or equal
%                     to the number of covariance functions N) [2]
%      sigma2       - 1D array of length max_deg defining the variance for
%                     each order of interaction [0.1*ones]
%      sigma2_prior - prior for sigma2 [prior_logunif]
%
%    For example N = 3, max_deg = 2:
%      k1 = CF_1 + CF_2 + CF_3
%      k2 = CF_1*CF_2 + CF_1*CF_3 + CF_2*CF_3
%      k3 = CF_1*CF_2*CF_3
%      GPCF = sigam2(1)*k1 + sigam2(2)*k2,
%
%  See also
%    GP_SET, GPCF_*
%
%  References:
%    Duvenaud, D. K., Nickisch, H., & Rasmussen, C. E. (2011). Additive
%    Gaussian processes. In Advances in neural information processing
%    systems (pp. 226-234).

% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2014 Tuomas Sivula

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

  ip=inputParser;
  ip.FunctionName = 'GPCF_ADDITIVE';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('cf',[], @iscell);
  ip.addParamValue('max_deg', 2, @(x) isscalar(x) && isnumeric(x) ...
                                      && x>0 && mod(x,1) == 0);
  ip.addParamValue('sigma2', [], @(x) isnumeric(x) && isvector(x) && all(x>0));
  ip.addParamValue('sigma2_prior', prior_logunif(), ...
                   @(x) isstruct(x) || isempty(x))
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

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

  % Initialize parameters

  % Kernels
  if init || ~ismember('cf',ip.UsingDefaults)
    if length(ip.Results.cf) < 2
      error('At least two covariance functions has to be given in cf');
    end
    gpcf.cf = ip.Results.cf;
  end
  ncf = length(gpcf.cf);

  % Max degree
  if init || ~ismember('max_deg',ip.UsingDefaults)
    if ip.Results.max_deg <= ncf
      gpcf.max_deg = ip.Results.max_deg;
    else
      warning('max_deg in additive kernel can not be greater than number of dimensions, max_deg truncated.');
      gpcf.max_deg = ncf;
    end
  end

  % Degree variances
  if init || ~ismember('sigma2',ip.UsingDefaults)
    if isempty(ip.Results.sigma2)
      gpcf.sigma2 = 0.1*ones(1,gpcf.max_deg);
    elseif length(ip.Results.sigma2) == gpcf.max_deg
      gpcf.sigma2 = ip.Results.sigma2;
    else
      error('Wrong number of elements in degree variance parameter vector sigma2')
    end
    % Ensure the right direction
    if size(gpcf.sigma2,1) ~= 1
      gpcf.sigma2 = gpcf.sigma2';
    end
  end

  % Degree variance priors
  if init || ~ismember('sigma2_prior',ip.UsingDefaults)
    gpcf.p.sigma2 = ip.Results.sigma2_prior;
  end

  % Ensure sigma2 matches max_deg
  if length(gpcf.sigma2) ~= gpcf.max_deg
    error('Parameters sigma2 and max_deg sizes does not match')
  end

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_additive_pak;
    gpcf.fh.unpak = @gpcf_additive_unpak;
    gpcf.fh.lp = @gpcf_additive_lp;
    gpcf.fh.lpg = @gpcf_additive_lpg;
    gpcf.fh.cfg = @gpcf_additive_cfg;
    gpcf.fh.ginput = @gpcf_additive_ginput;
    gpcf.fh.cov = @gpcf_additive_cov;
    gpcf.fh.trcov  = @gpcf_additive_trcov;
    gpcf.fh.trvar  = @gpcf_additive_trvar;
    gpcf.fh.recappend = @gpcf_additive_recappend;
  end

end

function [w, s, h] = gpcf_additive_pak(gpcf)
%GPCF_ADDITIVE_PAK  Combine GP covariance function parameters into one vector
%
%  Description
%    W = GPCF_ADDITIVE_PAK(GPCF, W) loops through all the covariance
%    functions and packs their parameters into one vector as
%    described in the respective functions. This is a mandatory
%    subfunction used for example in energy and gradient computations.
%
%  See also
%    GPCF_ADDITIVE_UNPAK

  ncf = length(gpcf.cf);
  w = []; s = {}; h=[];

  if ~isempty(gpcf.p.sigma2)
    w = [w log(gpcf.sigma2)];
    s = [s; sprintf('log(additive.sigma2 x %d)',numel(gpcf.sigma2))];
    h = [h ones(1,numel(gpcf.sigma2))];
    % Hyperparameters of lengthScale
    [wh, sh, hh] = gpcf.p.sigma2.fh.pak(gpcf.p.sigma2);
    sh=strcat(repmat('prior-', size(sh,1),1),sh);
    w = [w wh];
    s = [s; sh];
    h = [h 1+hh];
  end

  for i=1:ncf
    cf = gpcf.cf{i};
    [wi, si, hi] = cf.fh.pak(cf);
    w = [w wi];
    s = [s; si];
    h = [h hi];
  end
end

function [gpcf, w] = gpcf_additive_unpak(gpcf, w)
%GPCF_ADDITIVE_UNPAK  Sets the covariance function parameters into
%                 the structures
%
%  Description
%    [GPCF, W] = GPCF_ADDITIVE_UNPAK(GPCF, W) loops through all the
%    covariance functions and unpacks their parameters from W to
%    each covariance function structure. This is a mandatory
%    subfunction used for example in energy and gradient computations.
%
%  See also
%    GPCF_ADDITIVE_PAK
%
  ncf = length(gpcf.cf);

  if isfield(gpcf.p,'sigma2') && ~isempty(gpcf.p.sigma2)
    i1=1;
    i2=length(gpcf.sigma2);
    gpcf.sigma2 = exp(w(i1:i2));
    w = w(i2+1:end);
    % Hyperparameters of lengthScale
    [p, w] = gpcf.p.sigma2.fh.unpak(gpcf.p.sigma2, w);
    gpcf.p.sigma2 = p;
  end

  for i=1:ncf
    cf = gpcf.cf{i};
    [cf, w] = cf.fh.unpak(cf, w);
    gpcf.cf{i} = cf;
  end

end

function lp = gpcf_additive_lp(gpcf)
%GPCF_ADDITIVE_LP  Evaluate the log prior of covariance function parameters
%
%  Description
%    LP = GPCF_ADDITIVE_LP(GPCF, X, T) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%    GPCF_ADDITIVE_PAK, GPCF_ADDITIVE_UNPAK, GPCF_ADDITIVE_LPG, GP_E

  lp = 0;
  gpp = gpcf.p;

  if isfield(gpp,'sigma2') && ~isempty(gpp.sigma2)
    lp = lp +gpp.sigma2.fh.lp(gpcf.sigma2, gpp.sigma2) +sum(log(gpcf.sigma2));
  end

  ncf = length(gpcf.cf);
  for i=1:ncf
    cf = gpcf.cf{i};
    lp = lp + cf.fh.lp(cf);
  end

end

function lpg = gpcf_additive_lpg(gpcf)
%GPCF_ADDITIVE_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_ADDITIVE_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_ADDITIVE_PAK, GPCF_ADDITIVE_UNPAK, GPCF_ADDITIVE_LP, GP_G

  % Evaluate the gradients
  lpg = [];
  gpp=gpcf.p;

  if isfield(gpp,'sigma2') && ~isempty(gpp.sigma2)
    lll = length(gpcf.sigma2);
    lpgs = gpp.sigma2.fh.lpg(gpcf.sigma2, gpp.sigma2);
    lpg = [lpg lpgs(1:lll).*gpcf.sigma2+1 lpgs(lll+1:end)];
  end

  ncf = length(gpcf.cf);
  for i=1:ncf
    cf = gpcf.cf{i};
    lpg=[lpg cf.fh.lpg(cf)];
  end

end

function es = degrees(r, zs, inds)
% DEGREEs - Internal function used to calculate the covariance degree
% components e_k(z_inds(1), ..., z_inds(n))
%
%   Parameters:
%     r    - max degree
%     zs   - all the covariances of the underlying kernels
%     inds - indexes of the used kernels (empty means all)
%
%   N.B. Here all Cs and inds are given instead using slicing because the
%   latter is memory inefficient. Compare the memory usage in the
%   following:
%
%     inds = 1:100;
%     zs = rand(1000,1000,100);
%
%     % Memory inefficient
%     t = sum(zs(:,:,inds(inds~=13)),3);
%
%     % Memory efficient
%     t = zeros(1000,1000);
%     for i = inds(inds~=13)
%       t = t + zs(:,:,i);
%     end
%

  if isempty(inds)

    ncf = size(zs,3);

    if r == 1
      % Only one degree (implemented here for completeness)
      es = sum(zs,3);

    elseif r == 2
      % Only two degrees (implemented here for completeness)
      es = zeros(size(zs,1),size(zs,2),r);
      es(:,:,1) = sum(zs,3);
      for i1 = 1:ncf-1
        for i2 = i1+1:ncf
          es(:,:,2) = es(:,:,2) + zs(:,:,i1).*zs(:,:,i2);
        end
      end

    elseif r > 2
      % Over two degrees ... use Newton-Girard formulae
      sk = zeros(size(zs,1),size(zs,2),r);
      sk(:,:,1) = sum(zs,3);
      for i = 2:r
        sk(:,:,i) = sum(zs.^i, 3);
      end
      es = zeros(size(zs,1),size(zs,2),r);
      es(:,:,1) = sk(:,:,1);
      for i = 2:r
        presign = true;
        for k = 1:i-1
          if presign
            es(:,:,i) = es(:,:,i) + es(:,:,i-k).*sk(:,:,k);
            presign = false;
          else
            es(:,:,i) = es(:,:,i) - es(:,:,i-k).*sk(:,:,k);
            presign = true;
          end
        end
        if presign
          es(:,:,i) = es(:,:,i) + sk(:,:,i);
        else
          es(:,:,i) = es(:,:,i) - sk(:,:,i);
        end
        es(:,:,i) = es(:,:,i)./i;
      end

    else
      error('Invalid max_deg parameter')
    end

  else
    % Use only selected kernels

    % Check direction
    if r > length(inds)
      error('Invalid max degree relative to inds')
    end
    if size(inds,1) ~= 1
      inds = inds';
      if size(inds,1) ~= 1
        error('Invalid parameter inds')
      end
    end

    if r == 1
      % Only one degree
      es = 0;
      for i = inds
        es = es + zs(:,:,i);
      end

    elseif r == 2
      % Only two degrees
      es = zeros(size(zs,1),size(zs,2),r);
      for i = inds
        es(:,:,1) = es(:,:,1) + zs(:,:,i);
      end
      for i1 = 1:length(inds)-1
        for i2 = i1+1:length(inds)
          es(:,:,2) = es(:,:,2) + zs(:,:,inds(i1)).*zs(:,:,inds(i2));
        end
      end

    elseif r > 2
      % Over two degrees ... use Newton-Girard formulae
      sk = zeros(size(zs,1),size(zs,2),r);
      for i = inds
        sk(:,:,1) = sk(:,:,1) + zs(:,:,i);
      end
      for i = 2:r
        for j = inds
          sk(:,:,i) = sk(:,:,i) + zs(:,:,j).^i;
        end
      end
      es = zeros(size(zs,1),size(zs,2),r);
      es(:,:,1) = sk(:,:,1);
      for i = 2:r
        presign = true;
        for k = 1:i-1
          if presign
            es(:,:,i) = es(:,:,i) + es(:,:,i-k).*sk(:,:,k);
            presign = false;
          else
            es(:,:,i) = es(:,:,i) - es(:,:,i-k).*sk(:,:,k);
            presign = true;
          end
        end
        if presign
          es(:,:,i) = es(:,:,i) + sk(:,:,i);
        else
          es(:,:,i) = es(:,:,i) - sk(:,:,i);
        end
        es(:,:,i) = es(:,:,i)./i;
      end

    else
      error('Invalid max_deg parameter')
    end

  end

end

function DKff = gpcf_additive_cfg(gpcf, x, x2, mask, i1)
%GPCF_ADDITIVE_CFG  Evaluate gradient of covariance function
%               with respect to the parameters.
%
%  Description
%    DKff = GPCF_ADDITIVE_CFG(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_ADDITIVE_CFG(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_ADDITIVE_CFG(GPCF, X, [], MASK) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the diagonal of gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_ADDITIVE_CFG(GPCF, X, X2, [], i) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    hyperparameter. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%    GPCF_ADDITIVE_PAK, GPCF_ADDITIVE_UNPAK, GPCF_ADDITIVE_LP, GP_G

  [n, m] =size(x);
  ncf = length(gpcf.cf);
  r = gpcf.max_deg;

  DKff = {};

  if nargin==5
    % Use memory save option
    savememory=1;
    i3 = zeros(1,ncf);
    for k=1:ncf
      % Number of hyperparameters for each covariance function
      cf = gpcf.cf{k};
      i3(k) = cf.fh.cfg(cf,[],[],[],0);
    end
    if i1==0
      % Return number of hyperparameters
      DKff = sum(i3) + gpcf.max_deg;
      return
    end
    if i1 > gpcf.max_deg
      i1 = i1 - gpcf.max_deg;
      % The parameter belongs to one of the underlying kernels
      % Now i1 is [kernel_index, param_index_in_that_kernel]
      i3 = cumsum(i3);
      ind = find(i3 >= i1, 1);
      if ind > 1
        i1 = [ind, i1-i3(ind-1)];
      else
        i1 = [ind, i1];
      end
    end
  else
    savememory=0;
  end

  % Evaluate: DKff{1} = d Kff / d magnSigma2
  %           DKff{2} = d Kff / d lengthScale
  % NOTE! Here we have already taken into account that the parameters are transformed
  % through log() and thus dK/dlog(p) = p * dK/dp

  % evaluate the gradient for training covariance
  if nargin == 2 || (isempty(x2) && isempty(mask))

    % evaluate the individual covariance functions
    zs = zeros(n,n,ncf);
    for i=1:ncf
      cf = gpcf.cf{i};
      zs(:,:,i) = cf.fh.trcov(cf, x(:,i));
    end

    % Evaluate the gradients
    ind = 1:ncf;

    if ~savememory

      DKff = {};

      % Order variances sigma2
      es = degrees(r, zs, []);
      for i = 1:r
        % dlog(p) ... See NOTE above
        DKff{end+1} = gpcf.sigma2(i).*es(:,:,i);
      end

      % Subkernel hyperparameters
      for i=1:ncf
        cf = gpcf.cf{i};
        DK = cf.fh.cfg(cf, x(:,i));

        CC = gpcf.sigma2(1).*ones(n,n);
        if r > 1
          es = degrees(r-1, zs, ind(ind~=i));
          for j = 2:r
            CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
          end
        end
        for j = 1:length(DK)
          DKff{end+1} = DK{j}.*CC;
        end
      end

    else
      if length(i1) == 1
        % Order variance sigma2
        es = degrees(i1, zs, []);
        % dlog(p) ... See NOTE above
        DKff = gpcf.sigma2(i1).*es(:,:,end);
      else

        cf = gpcf.cf{i1(1)};
        DK = cf.fh.cfg(cf,x(:,i1(1)),[],[],i1(2));

        CC = gpcf.sigma2(1).*ones(n,n);
        if r > 1
          es = degrees(r-1, zs, ind(ind~=i1(1)));
          for j = 2:r
            CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
          end
        end
        DKff = DK.*CC;

      end

    end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 || isempty(mask)
    if size(x,2) ~= size(x2,2)
      error('gpcf_prod -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end

    % evaluate the individual covariance functions
    zs = zeros(size(x,1),size(x2,1),ncf);
    for i=1:ncf
      cf = gpcf.cf{i};
      zs(:,:,i) = cf.fh.cov(cf, x(:,i), x2(:,i));
    end

    % Evaluate the gradients
    ind = 1:ncf;

    if ~savememory

      DKff = {};

      % Order variances sigma2
      es = degrees(r, zs, []);
      for i = 1:r
        % dlog(p) ... See NOTE above
        DKff{end+1} = gpcf.sigma2(i).*es(:,:,i);
      end

      % Subkernel hyperparameters
      for i=1:ncf
        cf = gpcf.cf{i};
        DK = cf.fh.cfg(cf, x(:,i), x2(:,i));

        CC = gpcf.sigma2(1).*ones(n,n);
        if r > 1
          es = degrees(r-1, zs, ind(ind~=i));
          for j = 2:r
            CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
          end
        end
        for j = 1:length(DK)
          DKff{end+1} = DK{j}.*CC;
        end
      end

    else
      if length(i1) == 1
        % Order variance sigma2
        es = degrees(i1, zs, []);
        % dlog(p) ... See NOTE above
        DKff = gpcf.sigma2(i1).*es(:,:,end);
      else

        cf = gpcf.cf{i1(1)};
        DK = cf.fh.cfg(cf, x(:,i), x2(:,i), [], i1(2));

        CC = gpcf.sigma2(1).*ones(n,n);
        if r > 1
          es = degrees(r-1, zs, ind(ind~=i1(1)));
          for j = 2:r
            CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
          end
        end
        DKff = DK.*CC;

      end

    end

    % Evaluate: DKff{1}    = d mask(Kff,I) / d magnSigma2
    %           DKff{2...} = d mask(Kff,I) / d lengthScale
  elseif nargin == 4 || nargin == 5

    % evaluate the individual covariance functions
    zs = zeros(size(x,1),ncf);
    for i=1:ncf
      cf = gpcf.cf{i};
      zs(:,i) = cf.fh.trvar(cf, x(:,i));
    end

    % Evaluate the gradients
    ind = 1:ncf;

    if ~savememory

      DKff = {};

      % Order variances sigma2
      es = degrees(r, zs, []);
      for i = 1:r
        % dlog(p) ... See NOTE above
        DKff{end+1} = gpcf.sigma2(i).*es(:,:,i);
      end

      % Subkernel hyperparameters
      for i=1:ncf
        cf = gpcf.cf{i};
        DK = cf.fh.cfg(cf, x(:,i), [], 1);

        CC = gpcf.sigma2(1).*ones(n,n);
        if r > 1
          es = degrees(r-1, zs, ind(ind~=i));
          for j = 2:r
            CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
          end
        end
        for j = 1:length(DK)
          DKff{end+1} = DK{j}.*CC;
        end
      end

    else
      if length(i1) == 1
        % Order variance sigma2
        es = degrees(i1, zs, []);
        % dlog(p) ... See NOTE above
        DKff = gpcf.sigma2(i1).*es(:,:,end);
      else

        cf = gpcf.cf{i1(1)};
        DK = cf.fh.cfg(cf, x(:,i1(1)), [], 1, i1(2));

        CC = gpcf.sigma2(1).*ones(n,n);
        if r > 1
          es = degrees(r-1, zs, ind(ind~=i1(1)));
          for j = 2:r
            CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
          end
        end
        DKff = DK.*CC;

      end

    end

  end

end

function DKff = gpcf_additive_ginput(gpcf, x, x2, i1)
%GPCF_ADDITIVE_GINPUT  Evaluate gradient of covariance function with
%                  respect to x
%
%  Description
%    DKff = GPCF_ADDITIVE_GINPUT(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (cell array with matrix elements). This subfunction
%    is needed when computing gradients with respect to inducing
%    inputs in sparse approximations.
%
%    DKff = GPCF_ADDITIVE_GINPUT(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_ADDITIVE_GINPUT(GPCF, X, X2, i) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    covariate in X (cell array with matrix elements). This
%    subfunction is needed when using memory save option in
%    gp_set.
%
%  See also
%    GPCF_ADDITIVE_PAK, GPCF_ADDITIVE_UNPAK, GPCF_ADDITIVE_LP, GP_G

  [n, m] =size(x);
  r = gpcf.max_deg;

  if nargin==4
    % Use memory save option
    savememory=1;
    if i1==0
      % Return number of covariates
      if isfield(gpcf,'selectedVariables')
        DKff=length(gpcf.selectedVariables);
      else
        DKff=m;
      end
      return
    end
  else
    savememory=0;
  end
  % evaluate the gradient for training covariance
  if nargin == 2 || isempty(x2)

    ncf = length(gpcf.cf);

    % evaluate the individual covariance functions
    zs = zeros(n,n,ncf);
    for i=1:ncf
      cf = gpcf.cf{i};
      zs(:,:,i) = cf.fh.trcov(cf, x(:,i));
    end

    % Evaluate the gradients
    ind = 1:ncf;
    if ~savememory
      DKff=cellfun(@(a) zeros(n,n), cell(1,numel(x)), 'UniformOutput', 0);
    else
      DKff=cellfun(@(a) zeros(n,n), cell(1,n), 'UniformOutput', 0);
    end
    for i=1:ncf
      cf = gpcf.cf{i};
      if ~savememory
        DK = cf.fh.ginput(cf, x(:,i));
      else
        DK = cf.fh.ginput(cf, x(:,i), [], i1);
      end

      CC = gpcf.sigma2(1).*ones(n,n);
      if r > 1
        es = degrees(r-1, zs, ind(ind~=i));
        for j = 2:r
          CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
        end
      end
      for j = 1:length(DK)
        DKff{j} = DKff{j} + DK{j}.*CC;
      end
    end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 ||  nargin == 4
    if size(x,2) ~= size(x2,2)
      error('gpcf_prod -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end

    ncf = length(gpcf.cf);

    % evaluate the individual covariance functions
    zs = zeros(n,n,ncf);
    for i=1:ncf
      cf = gpcf.cf{i};
      zs(:,:,i) = cf.fh.cov(cf, x(:,i), x2(:,i));
    end

    % Evaluate the gradients
    ind = 1:ncf;
    if ~savememory
      DKff=cellfun(@(a) zeros(n,n), cell(1,numel(x)), 'UniformOutput', 0);
    else
      DKff=cellfun(@(a) zeros(n,n), cell(1,n), 'UniformOutput', 0);
    end
    for i=1:ncf
      cf = gpcf.cf{i};
      if ~savememory
        DK = cf.fh.ginput(cf, x(:,i), x2(:,i));
      else
        DK = cf.fh.ginput(cf, x(:,i), x2(:,i), i1);
      end

      CC = gpcf.sigma2(1).*ones(n,n);
      if r > 1
        es = degrees(r-1, zs, ind(ind~=i));
        for j = 2:r
          CC = CC + gpcf.sigma2(j).*es(:,:,j-1);
        end
      end
      for j = 1:length(DK)
        DKff{j} = DKff{j} + DK{j}.*CC;
      end
    end
  end

end

function C = gpcf_additive_cov(gpcf, x1, x2)
%GP_ADDITIVE_COV  Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_ADDITIVE_COV(GP, TX, X) takes in covariance function of a
%    Gaussian process GP and two matrixes TX and X that contain
%    input vectors to GP. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i in TX
%    and j in X. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%
%  See also
%    GPCF_ADDITIVE_TRCOV, GPCF_ADDITIVE_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  ncf = length(gpcf.cf);

  if m1~=ncf
    error('input dimension does not match with number of additive kernels')
  end

  r = gpcf.max_deg;

  if r == 1
    % Only one degree
    C = 0;
    for d = 1:ncf
      cf = gpcf.cf{d};
      C = C + cf.fh.cov(cf, x1(:,d), x2(:,d));
    end
    C = C.*gpcf.sigma2(1);

  elseif r == 2
    % Only two degrees
    zs = zeros(n1,n2,ncf);
    for d = 1:ncf
      cf = gpcf.cf{d};
      zs(:,:,d) = cf.fh.cov(cf, x1(:,d), x2(:,d));
    end
    C = 0;
    for i1 = 1:ncf-1
      for i2 = i1+1:ncf
        C = C + zs(:,:,i1).*zs(:,:,i2);
      end
    end
    C = C.*gpcf.sigma2(2);
    C = C + sum(zs,3).*gpcf.sigma2(1);

  elseif r > 2
    % Over two degrees ... use Newton-Girard formulae
    zs = zeros(n1,n2,ncf);
    for d = 1:ncf
      cf = gpcf.cf{d};
      zs(:,:,d) = cf.fh.cov(cf, x1(:,d), x2(:,d));
    end
    sk = zeros(n1,n2,r);
    sk(:,:,1) = sum(zs,3);
    for i = 2:r
      sk(:,:,i) = sum(zs.^i, 3);
    end
    clear zs
    es = zeros(n1,n2,r);
    es(:,:,1) = sk(:,:,1);
    for i = 2:r
      presign = true;
      for k = 1:i-1
        if presign
          es(:,:,i) = es(:,:,i) + es(:,:,i-k).*sk(:,:,k);
          presign = false;
        else
          es(:,:,i) = es(:,:,i) - es(:,:,i-k).*sk(:,:,k);
          presign = true;
        end
      end
      if presign
        es(:,:,i) = es(:,:,i) + sk(:,:,i);
      else
        es(:,:,i) = es(:,:,i) - sk(:,:,i);
      end
      es(:,:,i) = es(:,:,i)./i;
    end
    for i = 1:r
      es(:,:,i) = es(:,:,i).*gpcf.sigma2(i);
    end
    C = sum(es,3);

  else
    error('Invalid max_deg parameter')
  end

end

function C = gpcf_additive_trcov(gpcf, x)
%GP_ADDITIVE_TRCOV     Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_ADDITIVE_TRCOV(GP, TX) takes in covariance function of a
%    Gaussian process GP and matrix TX that contains training
%    input vectors. Returns covariance matrix C. Every element ij
%    of C contains covariance between inputs i and j in TX. This
%    is a mandatory subfunction used for example in prediction and
%    energy computations.
%
%  See also
%    GPCF_ADDITIVE_COV, GPCF_ADDITIVE_TRVAR, GP_COV, GP_TRCOV

  [n,m]=size(x);

  ncf = length(gpcf.cf);

  if m~=ncf
    error('input dimension does not match with number of additive kernels')
  end

  r = gpcf.max_deg;

  if r == 1
    % Only one degree
    C = 0;
    for d = 1:ncf
      cf = gpcf.cf{d};
      C = C + cf.fh.trcov(cf, x(:,d));
    end
    C = C.*gpcf.sigma2(1);

  elseif r == 2
    % Only two degrees
    zs = zeros(n,n,ncf);
    for d = 1:ncf
      cf = gpcf.cf{d};
      zs(:,:,d) = cf.fh.trcov(cf, x(:,d));
    end
    C = 0;
    for i1 = 1:ncf-1
      for i2 = i1+1:ncf
        C = C + zs(:,:,i1).*zs(:,:,i2);
      end
    end
    C = C.*gpcf.sigma2(2);
    C = C + sum(zs,3).*gpcf.sigma2(1);

  elseif r > 2
    % Over two degrees ... use Newton-Girard formulae
    zs = zeros(n,n,ncf);
    for d = 1:ncf
      cf = gpcf.cf{d};
      zs(:,:,d) = cf.fh.trcov(cf, x(:,d));
    end
    sk = zeros(n,n,r);
    sk(:,:,1) = sum(zs,3);
    for i = 2:r
      sk(:,:,i) = sum(zs.^i, 3);
    end
    clear zs
    es = zeros(n,n,r);
    es(:,:,1) = sk(:,:,1);
    for i = 2:r
      presign = true;
      for k = 1:i-1
        if presign
          es(:,:,i) = es(:,:,i) + es(:,:,i-k).*sk(:,:,k);
          presign = false;
        else
          es(:,:,i) = es(:,:,i) - es(:,:,i-k).*sk(:,:,k);
          presign = true;
        end
      end
      if presign
        es(:,:,i) = es(:,:,i) + sk(:,:,i);
      else
        es(:,:,i) = es(:,:,i) - sk(:,:,i);
      end
      es(:,:,i) = es(:,:,i)./i;
    end
    for i = 1:r
      es(:,:,i) = es(:,:,i).*gpcf.sigma2(i);
    end
    C = sum(es,3);

  else
    error('Invalid max_deg parameter')
  end
end

function C = gpcf_additive_trvar(gpcf, x)
% GP_ADDITIVE_TRVAR     Evaluate training variance vector
%
%  Description
%    C = GP_ADDITIVE_TRVAR(GPCF, TX) takes in covariance function of
%    a Gaussian process GPCF and matrix TX that contains training
%    inputs. Returns variance vector C. Every element i of C
%    contains variance of input i in TX. This is a mandatory
%    subfunction used for example in prediction and energy computations.
%
%  See also
%    GPCF_ADDITIVE_COV, GP_COV, GP_TRCOV

  [n,m]=size(x);

  ncf = length(gpcf.cf);

  if m~=ncf
    error('input dimension does not match with number of additive kernels')
  end

  r = gpcf.max_deg;

  if r == 1
    % Only one degree
    C = 0;
    for d = 1:ncf
      cf = gpcf.cf{d};
      C = C + cf.fh.trvar(cf, x(:,d));
    end
    C = C.*gpcf.sigma2(1);

  elseif r == 2
    % Only two degrees
    zs = zeros(n,ncf);
    for d = 1:ncf
      cf = gpcf.cf{d};
      zs(:,d) = cf.fh.trvar(cf, x(:,d));
    end
    C = 0;
    for i1 = 1:ncf-1
      for i2 = i1+1:ncf
        C = C + zs(:,i1).*zs(:,i2);
      end
    end
    C = C.*gpcf.sigma2(2);
    C = C + sum(zs,2).*gpcf.sigma2(1);

  elseif r > 2
    % Over two degrees ... use Newton-Girard formulae
    zs = zeros(n,ncf);
    for d = 1:ncf
      cf = gpcf.cf{d};
      zs(:,d) = cf.fh.trvar(cf, x(:,d));
    end
    sk = zeros(n,r);
    sk(:,1) = sum(zs,2);
    for i = 2:r
      sk(:,i) = sum(zs.^i, 2);
    end
    clear zs
    es = zeros(n,r);
    es(:,1) = sk(:,1);
    for i = 2:r
      presign = true;
      for k = 1:i-1
        if presign
          es(:,i) = es(:,i) + es(:,i-k).*sk(:,k);
          presign = false;
        else
          es(:,i) = es(:,i) - es(:,i-k).*sk(:,k);
          presign = true;
        end
      end
      if presign
        es(:,i) = es(:,i) + sk(:,i);
      else
        es(:,i) = es(:,i) - sk(:,i);
      end
      es(:,i) = es(:,i)./i;
    end
    for i = 1:r
      es(:,i) = es(:,i).*gpcf.sigma2(i);
    end
    C = sum(es,2);

  else
    error('Invalid max_deg parameter')
  end

end

function reccf = gpcf_additive_recappend(reccf, ri, gpcf)
%RECAPPEND  Record append
%
%  Description
%    RECCF = GPCF_ADDITIVE_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains all
%    the old samples and the current samples from GPCF. This
%    subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC, GP_MC->RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_additive';

    % Initialize parameters
    reccf.sigma2=[];

    ncf = length(ri.cf);
    for i=1:ncf
      cf = ri.cf{i};
      reccf.cf{i} = cf.fh.recappend([], ri.cf{i});
    end

    % Set the function handles
    reccf.fh.pak = @gpcf_additive_pak;
    reccf.fh.unpak = @gpcf_additive_unpak;
    reccf.fh.lp = @gpcf_additive_lp;
    reccf.fh.lpg = @gpcf_additive_lpg;
    reccf.fh.cfg = @gpcf_additive_cfg;
    reccf.fh.cov = @gpcf_additive_cov;
    reccf.fh.trcov  = @gpcf_additive_trcov;
    reccf.fh.trvar  = @gpcf_additive_trvar;
    reccf.fh.recappend = @gpcf_additive_recappend;

    % Set other
    reccf.max_deg = ri.max_deg;
    reccf.p=[];
    reccf.p.sigma2=[];
    if isfield(ri.p,'sigma2') && ~isempty(ri.p.sigma2)
      reccf.p.sigma2 = ri.p.sigma2;
    end

  else
    % Append to the record
    gpp = gpcf.p;

    % record sigma2
    reccf.sigma2(ri,:) = gpcf.sigma2;
    if isfield(gpp,'sigma2') && ~isempty(gpp.sigma2)
      reccf.p.sigma2 = gpp.sigma2.fh.recappend(reccf.p.sigma2, ri, gpcf.p.sigma2);
    end

    % Loop over all of the covariance functions
    ncf = length(gpcf.cf);
    for i=1:ncf
      cf = gpcf.cf{i};
      reccf.cf{i} = cf.fh.recappend(reccf.cf{i}, ri, cf);
    end
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
function clusterLabels = gacMerging(graphW, initClusters, groupNumber, strDescr, z)
%% Cluster merging for Graph Agglomerative Clustering
% Implements an agglomerative clustering algorithm based on maiximum graph
%   strcutural affinity of two groups
% Inputs:
%	- graphW: asymmetric weighted adjacency matrix
%   - initClusters: a cell array of clustered vertices
%   - groupNumber: the final number of clusters
%   - strDescr: structural descriptor, 'zeta' or 'path'
%   - z: (I - z*P), default: 0.01
% Outputs:
%   - clusterLabels: 1 x m list whose i-th entry is the group assignment of
%                   the i-th data vector w_i. Groups are indexed
%                   sequentially, starting from 1.
% by Wei Zhang (wzhang009 at gmail.com), June, 8, 2011

%%
numSample = size(graphW,1);
IminuszW = eye(numSample) - z*graphW;
clear graphW
myInf = 1e10;

%% initialization
VERBOSE = true;
switch lower(strDescr)
    case 'zeta'
        complexity_fun = @gacZetaEntropy;
        conditionalComplexity_fun = @gacZetaCondEntropy;
    case 'path'
        complexity_fun = @gacPathEntropy;
        conditionalComplexity_fun = @gacPathCondEntropy;
    otherwise
        error('GAC: Descriptor type is not supported!');
end

numClusters = length(initClusters);
if numClusters <= groupNumber
    error('GAC: too few initial clusters. Do not need merging!');
end

%% compute the structural complexity of each initial cluster
clusterComp = zeros(numClusters,1);
for i = 1 : numClusters
    clusterComp(i) = complexity_fun(IminuszW(initClusters{i}, initClusters{i}));
end

%% compute initial (negative) affinity table (upper trianglar matrix), very slow
if VERBOSE
    disp('   Computing initial table.' );
end
affinityTab = Inf(numClusters);
for j = 1 : numClusters
    for i = 1 : j-1
        affinityTab(i, j) = - conditionalComplexity_fun(IminuszW, initClusters{i}, initClusters{j});
    end
end
affinityTab = bsxfun(@plus, clusterComp, clusterComp') + affinityTab;

if VERBOSE
    disp('   Starting merging process');
end

curGroupNum = numClusters;
while true
    if mod( curGroupNum, 20 ) == 0 && VERBOSE
        disp(['   Group count: ' num2str(curGroupNum)]);
    end
    % Find two clusters with the best affinity
    [minAff, minIndex1] = min(affinityTab(1:curGroupNum, 1:curGroupNum), [], 1);
    [~, minIndex2] = min(minAff);
    minIndex1 = minIndex1(minIndex2);
    if minIndex2 < minIndex1,  [minIndex1, minIndex2] = swap(minIndex1, minIndex2); end

    % merge the two clusters
    new_cluster = unique([initClusters{minIndex1}; initClusters{minIndex2}]);
    % move the second cluster to be merged to the end of the cluster array
    % note that we only need to copy the end cluster's information to
    % the second cluster's position
    if (minIndex2 ~= curGroupNum)
        initClusters{minIndex2} = initClusters{end};
        clusterComp(minIndex2) = clusterComp(curGroupNum);
        % affinityTab is an upper triangular matrix
        affinityTab(1:minIndex2-1, minIndex2) = affinityTab(1:minIndex2-1, curGroupNum);
        affinityTab(minIndex2, minIndex2+1:curGroupNum-1) = affinityTab(minIndex2+1:curGroupNum-1, curGroupNum);
    end

    % update the first cluster and remove the second cluster
    initClusters{minIndex1} = new_cluster;
    initClusters(end) = [];
    clusterComp(minIndex1) = complexity_fun(IminuszW(new_cluster, new_cluster));
    clusterComp(curGroupNum) = myInf;
    affinityTab(:,curGroupNum) = myInf;
    affinityTab(curGroupNum,:) = myInf;
    curGroupNum = curGroupNum - 1;
    if curGroupNum <= groupNumber
        break;
    end

    % update the affinity table for the merged cluster
    for groupIndex1 = 1:minIndex1-1
        affinityTab(groupIndex1, minIndex1) = - conditionalComplexity_fun(IminuszW, initClusters{groupIndex1}, new_cluster);
    end
    for groupIndex1 = minIndex1+1:curGroupNum
        affinityTab(minIndex1, groupIndex1) = - conditionalComplexity_fun(IminuszW, initClusters{groupIndex1}, new_cluster);
    end
    affinityTab(1:minIndex1-1, minIndex1) = clusterComp(1:minIndex1-1) + clusterComp(minIndex1) + affinityTab(1:minIndex1-1, minIndex1);
    affinityTab(minIndex1, minIndex1+1:curGroupNum) = clusterComp(minIndex1+1:curGroupNum)' + clusterComp(minIndex1) + affinityTab(minIndex1, minIndex1+1:curGroupNum);
end

%% generate sample labels
clusterLabels = ones(numSample,1);
for i = 1:length(initClusters)
    clusterLabels(initClusters{i}) = i;
end
if VERBOSE
    disp(['   Final group count: ' num2str(curGroupNum)]);
end

end

function [y, x] = swap (x, y)
end

author: jwyang


text:
%% Lane Marking Identification
% One component of future luxury automobile safety
% systems is warning drivers that they are drifting
% between lanes.  To do this, it is first necessary
% to identify where the lanes are.  In this example,
% we will look at identifying the lanes from a image
% of a road.
%
% Copyright 2007-2013 MathWorks, Inc.
%

%% Setup
clear

%% Initialize Objects
% Video Reader

hvfr = vision.VideoFileReader('viplanedeparture.avi', ...
    'VideoOutputDataType', 'uint8');

%%
% Create a |ColorSpaceConverter| System object to
% convert the RGB image to an intensity image.
hColor = vision.ColorSpaceConverter( ...
               'Conversion', 'RGB to intensity');

%%
% Create Blob Detecter

hBlob = vision.BlobAnalysis(            ...
    'MinimumBlobAreaSource',     'Property', ...
    'MinimumBlobArea',           10        , ...
    'MajorAxisLengthOutputPort', true      , ...
    'MinorAxisLengthOutputPort', true      , ...
    'OrientationOutputPort',     true      , ...
    'CentroidOutputPort',        true      , ...
    'BoundingBoxOutputPort',     true      , ...
    'LabelMatrixOutputPort',     false     , ...
    'AreaOutputPort',            false     );

%% Create Shape Inserter
% Create a |ShapeInserter| System object to draw
% lanes on the original image.
hShapes = vision.ShapeInserter( ...
                'Shape', 'Lines', ...
                'BorderColor', 'Custom', ...
                'CustomBorderColor', [0 255 0], ...
                'Antialiasing', true);

%% Configure Video Players
hVidSource = vision.VideoPlayer;
hVidBlobs  = vision.VideoPlayer;
hVidLanes  = vision.VideoPlayer;

%% Setup Replay Loop
NumLoops = 3;
for loopIdx = 1:NumLoops

    % Loop through video frames
    while ~isDone(hvfr)
        % Read Frame
        frameRGB = step(hvfr);
        % and display
        step(hVidSource, frameRGB);

        % Convert to Intensity
        frameGray = step(hColor, frameRGB);

        % Thresholding image using Otsu's method
        level = graythresh(frameGray) * intmax(class(frameGray));
        binary = frameGray > level;

        % Perform blob detection
        [centroid, bbox, Major, Minor, Orientation] = step(hBlob, binary);
        step(hVidBlobs, binary);

        % Find lanes
        laneIndex = findlanes(bbox, Major, Minor);

        numLanes = sum(laneIndex);

        % Mark Lanes on original image
        laneVertices = getLaneVerticies(bbox(laneIndex, :), Orientation(laneIndex));
        laneVertices = int32(laneVertices);

        laneImage = step(hShapes, frameRGB, laneVertices);
        step(hVidLanes, laneImage);

    % End While
    end

    reset(hvfr);
% end loop
end

%% Release resources
release(hvfr   );
release(hColor );
release(hBlob  );
release(hShapes);


text:
function [header_array, CateNames, CatLengths, preBaseline] = read_sbin_header(filename)

% READ_SBIN_HEADER reads the header information from an EGI segmented simple binary format file
%
% Use as
%   [header_array, CateNames, CatLengths, preBaseline] = read_sbin_header(filename)
% with
%   header_array     - differs between versions, read code for details
%   CateNames        - category names
%   CatLengths       - length of category names
%   preBaseline      - number of samples in the baseline prior to the baseline event
% and
%   filename    - the name of the data file
%
% Since there is no unique event code for the segmentation event, and hence the baseline period,
% the first event code in the list will be assumed to be the segmentation event.
% NetStation itself simply ignores possible baseline information when importing simple binary files.
%_______________________________________________________________________
%
%
% Modified from EGI's readEGLY.m with permission 2008-03-31 Joseph Dien
%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid=fopen_or_error([filename],'r');

version     = fread(fid,1,'int32');
if isempty(version)
    ft_error('ERROR:  This is not a simple binary file.  Note that NetStation does not successfully directly convert EGIS files to simple binary format.');
end;

%check byteorder
[str,maxsize,cEndian]=computer;
if version < 7
    if cEndian == 'B'
        endian = 'ieee-be';
    elseif cEndian == 'L'
        endian = 'ieee-le';
    end;
elseif (version > 6) && ~bitand(version,6)
    if cEndian == 'B'
        endian = 'ieee-le';
    elseif cEndian == 'L'
        endian = 'ieee-be';
    end;
    version = swapbytes(uint32(version));
else
    ft_error('ERROR:  This is not a simple binary file.  Note that NetStation does not successfully directly convert EGIS files to simple binary format.');
end;

if bitand(version,1) == 0
    unsegmented = 1;
else
    unsegmented = 0;
end;

precision = bitand(version,6);
if precision == 0
    ft_error('File precision is not defined.');
end;

%       read header...
year        = fread(fid,1,'int16',endian);
month       = fread(fid,1,'int16',endian);
day         = fread(fid,1,'int16',endian);
hour        = fread(fid,1,'int16',endian);
minute      = fread(fid,1,'int16',endian);
second      = fread(fid,1,'int16',endian);
millisecond = fread(fid,1,'int32',endian);
Samp_Rate   = fread(fid,1,'int16',endian);
NChan       = fread(fid,1,'int16',endian);
Gain        = fread(fid,1,'int16',endian);
Bits        = fread(fid,1,'int16',endian);
Range       = fread(fid,1,'int16',endian);

epochMarked=0;
if unsegmented,
    NumCategors = 0;
    NSegments   = 1;
    NSamples    = fread(fid,1,'int32',endian);
    NEvent      = fread(fid,1,'int16',endian);
    EventCodes='';
    for j = 1:NEvent
        EventCodes(j,1:4) = char(fread(fid,[1,4],'char',endian));
    end
    CateNames   = [];
    CatLengths  = [];
    preBaseline = 0;

    if any(strcmp('epoc',cellstr(EventCodes))) && any(strcmp('tim0',cellstr(EventCodes))) %actually epoch-marked segmented file format

        %Note that epoch-marked simple binary file format loses cell information so just the one "cell"

        switch precision
            case 2
                [temp,count]    = fread(fid,[NChan+NEvent, NSamples],'int16',endian);
            case 4
                [temp,count]    = fread(fid,[NChan+NEvent, NSamples],'single',endian);
            case 6
                [temp,count]    = fread(fid,[NChan+NEvent, NSamples],'double',endian);
        end
        eventData(:,1:NSamples)  = temp( (NChan+1):(NChan+NEvent), 1:NSamples);

        NSegments = length(find(eventData(find(strcmp('epoc',cellstr(EventCodes))),:)));
        epocSamples=find(eventData(find(strcmp('epoc',cellstr(EventCodes))),:));
        segmentLengths=diff(epocSamples);
        segmentLengths=[segmentLengths length(eventData)-epocSamples(end)+1];
        NSamples=max(segmentLengths); % samples per segment, will zero pad out to longest length
        preBaseline=min(find(eventData(find(strcmp('tim0',cellstr(EventCodes))),:)))-min(find(eventData(find(strcmp('epoc',cellstr(EventCodes))),:))); %make assumption all prestimulus durations are the same
    end;
else
    NumCategors = fread(fid,1,'int16',endian);
    for j = 1:NumCategors
        CatLengths(j)   = fread(fid,1,'int8',endian);
        for i = 1:CatLengths(j)
            CateNames{j}(i) = char(fread(fid,1,'char',endian));
        end
    end
    NSegments   = fread(fid,1,'int16',endian);
    NSamples    = fread(fid,1,'int32',endian);          % samples per segment
    NEvent      = fread(fid,1,'int16',endian);          % num events per segment
    EventCodes = '';
    for j = 1:NEvent
        EventCodes(j,1:4)   = char(fread(fid,[1,4],'char',endian));
    end

    preBaseline=0;
    if NEvent > 0

        for j = 1:NSegments
            [segHdr(j,1), count]    = fread(fid, 1,'int16',endian);    %cell
            [segHdr(j,2), count]    = fread(fid, 1,'int32',endian);    %time stamp
            switch precision
                case 2
                    [temp,count]    = fread(fid,[NChan+NEvent, NSamples],'int16',endian);
                case 4
                    [temp,count]    = fread(fid,[NChan+NEvent, NSamples],'single',endian);
                case 6
                    [temp,count]    = fread(fid,[NChan+NEvent, NSamples],'double',endian);
            end
            eventData(:,((j-1)*NSamples+1):j*NSamples)  = temp( (NChan+1):(NChan+NEvent), 1:NSamples);
        end
    end;

    if ~unsegmented
        if NEvent == 0
            theEvent=[];
        elseif NEvent == 1
            %assume this is the segmentation event
            theEvent=find(eventData(1,:)>0);
            theEvent=mod(theEvent(1)-1,NSamples)+1;
        else
            %assume the sample that always has an event is the baseline
            %if more than one, choose the earliest one
            totalEventSamples=mod(find(sum(eventData(:,:))')-1,NSamples)+1;
            baselineCandidates = unique(totalEventSamples);
            counters=hist(totalEventSamples, length(baselineCandidates));
            theEvent=min(baselineCandidates(find(ismember(counters,NSegments))));
        end;

        preBaseline=theEvent-1;
        if (preBaseline == -1)
            preBaseline =0;
        end;
        if isempty(preBaseline)
            preBaseline =0;
        end;
    end
end

fclose(fid);

header_array    = double([version year month day hour minute second millisecond Samp_Rate NChan Gain Bits Range NumCategors, NSegments, NSamples, NEvent]);

author: spm


text:
function [x, y] = intline(x1, x2, y1, y2)
%INTLINE Integer-coordinate line drawing algorithm.
%  [X, Y] = INTLINE(X1, X2, Y1, Y2) computes an
%  approximation to the line segment joining (X1, Y1) and
%  (X2, Y2) with integer coordinates.  X1, X2, Y1, and Y2
%  should be integers.  INTLINE is reversible; that is,
%  INTLINE(X1, X2, Y1, Y2) produces the same results as
%  FLIPUD(INTLINE(X2, X1, Y2, Y1)).

%   Copyright 1993-2002 The MathWorks, Inc.  Used with permission.
%   $Revision: 1.4 $  $Date: 2003/11/21 14:38:20 $

dx = abs(x2 - x1);
dy = abs(y2 - y1);

% Check for degenerate case.
if ((dx == 0) & (dy == 0))
   x = x1;
   y = y1;
  return;
end

flip = 0;
if (dx >= dy)
   if (x1 > x2)
      % Always "draw" from left to right.
      t = x1; x1 = x2; x2 = t;
      t = y1; y1 = y2; y2 = t;
      flip = 1;
   end
   m = (y2 - y1)/(x2 - x1);
   x = (x1:x2).';
   y = round(y1 + m*(x - x1));
else
   if (y1 > y2)
      % Always "draw" from bottom to top.
      t = x1; x1 = x2; x2 = t;
      t = y1; y1 = y2; y2 = t;
      flip = 1;
   end
   m = (x2 - x1)/(y2 - y1);
   y = (y1:y2).';
   x = round(x1 + m*(y - y1));
end

if (flip)
   x = flipud(x);
   y = flipud(y);
end

author: HuangCongQing


text:
function [z,flag] = formationP2z(polytope,zd,range,lc,lm)
max_aw = 1;
min_aw = 0.45;
polynum = length(polytope);
% get the A and b from polytope
Ap = []; bp = [];
for i=1:polynum
    Ap = [Ap;polytope(i).A];
    bp = [bp;polytope(i).b];
end
if length(zd)==12
    % constraints al<=a<=au and thl<=th<=thu just use lb and ub
    lb = [range.lb;-pi;min_aw*ones(4,1);-pi/4*ones(4,1);min_aw]';
    ub = [range.ub;pi;max_aw*ones(4,1);pi/4*ones(4,1);max_aw]';
    % calculate z
    [z,~,flag] = fmincon(@(z) (z-zd)*(z-zd)',zd,[],[],[],[],lb,ub,@(z) mycon(z,Ap,bp,lc,lm));
else
    % constraints al<=a<=au and thl<=th<=thu just use lb and ub
    lb = [range.lb;min_aw*ones(2,1)]';
    ub = [range.ub;max_aw*ones(2,1)]';
    % calculate z
    [z,~,flag] = fmincon(@(z) (z-zd)*(z-zd)',zd,[],[],[],[],lb,ub,@(z) mycon2(z,Ap,bp,lc,lm));
end
end

function [c,ceq] = mycon(z,A,b,lc,lm)
import rvctools.*
max_aw = 1;
min_aw = 0.45;
xr = [z(1:2),0,0,0,z(3)]; a = z(4:7); th = z(8:11); w = z(12);
xm = xr2m(xr,a,th,lc,lm);
for i=1:size(xm,1)
    plat(i) = pkgMechanics.RigidCuboid(1,xm(i,:),lm);
    pm_v(:,:,i) = plat(i).vertices(1:4,:)';
    Rm(:,:,i) = rpy2r(xm(i,4:end));
    pw_m(:,i) = xm(i,1:3)';
    pw_v(:,4*i-3:4*i) = pw_m(:,i)+Rm(:,:,i)*pm_v(:,:,i);
end
for i=1:size(pw_v,2)
    c1(:,i) = A*pw_v(1:2,i)-b;
end
c2 = (a.^2+w.^2-max_aw^2)';
c3 = -(a.^2+w.^2-min_aw^2)';
c = [c1(:);c2(:);c3(:)];
ceq = 0;
end

function [c,ceq] = mycon2(z,A,b,lc,lm)
import rvctools.*
max_aw = 1;
min_aw = 0.45;
xr = [z(1:2),0,0,0,0]; a = z(3)*ones(1,4); th = zeros(1,4); w = z(4);
xm = xr2m(xr,a,th,lc,lm);
for i=1:size(xm,1)
    plat(i) = pkgMechanics.RigidCuboid(1,xm(i,:),lm);
    pm_v(:,:,i) = plat(i).vertices(1:4,:)';
    Rm(:,:,i) = rpy2r(xm(i,4:end));
    pw_m(:,i) = xm(i,1:3)';
    pw_v(:,4*i-3:4*i) = pw_m(:,i)+Rm(:,:,i)*pm_v(:,:,i);
end
for i=1:size(pw_v,2)
    c1(:,i) = A*pw_v(1:2,i)-b;
end
c2 = (a.^2+w.^2-max_aw^2)';
c3 = -(a.^2+w.^2-min_aw^2)';
c = [c1(:);c2(:);c3(:)];
ceq = 0;
end

repo: paper-simulation