text:
%ETS2 Elementary transform sequence in 2D
%
% This class and package allows experimentation with sequences of spatial
% transformations in 2D.
%
%          import ETS2.*
%          a1 = 1; a2 = 1;
%          E = Rz('q1') * Tx(a1) * Rz('q2') * Tx(a2)
%
% Operation methods::
%   fkine      forward kinematics
%
% Information methods::
%   isjoint    test if transform is a joint
%   njoints    the number of joint variables
%   structure  a string listing the joint types
%
% Display methods::
%   display    display value as a string
%   plot       graphically display the sequence as a robot
%   teach      graphically display as robot and allow user control
%
% Conversion methods::
%   char       convert to string
%   string     convert to string with symbolic variables
%
% Operators::
%   *          compound two elementary transforms
%   +          compound two elementary transforms
%
% Notes::
% - The sequence is an array of objects of superclass ETS2, but with
%   distinct subclasses: Rz, Tx, Ty.
% - Use the command 'clear imports' after using ETS3.
%
%
% See also ETS3.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
classdef ETS2
    properties
        what     % type of transform (string): Rx, Ry, etc
        param    % the constant numerical parameter (if not joint)
        qvar     % the integer joint index 1..N (if joint)
        qlim     % for prismatic joint, a 2 vector [min,max]
    end

    methods
        function obj = ETS2(what, x, varargin)
            %ETS2.ETS2  Create an ETS2 object
            %
            % E = ETS2(W, V) is a new ETS2 object that defines an elementary transform where
            % W is 'Rz', 'Tx' or 'Ty' and V is the paramter for the transform.  If V is a string
            % of the form 'qN' where N is an integer then the transform is considered
            % to be a joint.  Otherwise the transform is a constant.
            %
            % E = ETS2(E1) is a new ETS2 object that is a clone of the ETS2 object E1.
            %
            % See also ETS2.Rz, ETS2.Tx, ETS2.Ty.

            assert(nargin > 0, 'RTB:ETS2:ETS2:badarg', 'no arguments given');

            opt.qlim = [];
            opt = tb_optparse(opt, varargin);

            obj.qvar = NaN;
            obj.param = 0;

            if ~isempty(opt.qlim)
                assert(length(opt.qlim) == 2, 'ETS2: qlim must be a 2-vector');
            end
            obj.qlim = opt.qlim;

            if nargin > 1
                if isa(x, 'ETS2')
                    % clone it
                    obj.what = x.what;
                    obj.qvar = x.qvar;
                    obj.param = x.param;
                else
                    % create a new one
                    assert(ismember(what, {'Tx','Ty','R','Rz'}), 'ETS2: invalid transform type given');
                    if strcmp(what, 'R')
                        what = 'Rz';
                    end
                    if ischar(x)
                        obj.qvar = str2num(x(2:end));
                    else
                        obj.param = x;
                    end
                    obj.what = what;
                end
            end
        end

        function r = fkine(ets, q, varargin)
            %ETS2.fkine Forward kinematics
            %
            % ETS.fkine(Q, OPTIONS) is the forward kinematics, the pose of the end of the
            % sequence as an SE2 object.  Q (1xN) is a vector of joint variables.
            %
            % ETS.fkine(Q, N, OPTIONS) as above but process only the first N elements
            % of the transform sequence.
            %
            % Options::
            %  'deg'     Angles are given in degrees.
            r = SE2;

            opt.deg = false;
            [opt,args] = tb_optparse(opt, varargin);

            if opt.deg
                opt.deg = pi/180;
            else
                opt.deg = 1;
            end

            n = length(ets);
            if ~isempty(args) && isreal(args{1})
                n = args{1};
            end
            assert(n>0 && n <= length(ets), 'RTB:ETS2:badarg', 'bad value of n given');

            for i=1:n
                e = ets(i);
                if e.isjoint
                    v = q(e.qvar);
                else
                    v = e.param;
                end
                switch e.what
                    case 'Tx'
                        r = r * SE2(v, 0, 0);
                    case 'Ty'
                        r = r * SE2(0, v, 0);
                    case 'Rz'
                        r = r * SE2(0, 0, v*opt.deg);
                end
            end
            r = r.simplify();  % simplify it if symbolic
        end

        function b = isjoint(ets)
            %ETS2.isjoint  Test if transform is a joint
            %
            % E.isjoint is true if the transform element is a joint, that is, its
            % parameter is of the form 'qN'.
            b = ~isnan(ets.qvar);
        end

        function v = isprismatic(ets)
            %ETS2.isprismatic  Test if transform is prismatic joint
            %
            % E.isprismatic is true if the transform element is a joint, that is, its
            % parameter is of the form 'qN' and it controls a translation.
            v = isjoint(ets) && (ets.what(1) == 'T');
        end

        function k = find(ets, j)
            %ETS2.find  Find joints in transform sequence
            %
            % E.find(J) is the index in the transform sequence ETS (1xN) corresponding
            % to the J'th joint.
            [~,k] = find([ets.qvar] == j);
        end

        function n = njoints(ets)
            %ETS2.njoints  Number of joints in transform sequence
            %
            % E.njoints is the number of joints in the transform sequence.
            %
            % See also ETS2.n.
            n = max([ets.qvar]);
        end

        function v = n(ets)
            %ETS2.n  Number of joints in transform sequence
            %
            % E.njoints is the number of joints in the transform sequence.
            %
            % Notes::
            % - Is a wrapper on njoints, for compatibility with SerialLink object.
            % See also ETS2.n.
            v = ets.njoints;
        end

        function s = string(ets)
            %ETS2.string  Convert to string with symbolic variables
            %
            % E.string is a string representation of the transform sequence where
            % non-joint parameters have symbolic names L1, L2, L3 etc.
            %
            % See also trchain.
            for i = 1:length(ets)
                e = ets(i);
                if e.isjoint
                    term = sprintf('%s(q%d)', e.what, e.qvar);
                else
                    term = sprintf('%s(L%d)', e.what, constant);
                    constant = constant + 1;
                end
                if i == 1
                    s = term;
                else
                    s = [s ' ' term];
                end
            end
        end

        function out = mtimes(ets1, ets2)
            %ETS2.mtimes Compound transforms
            %
            % E1 * E2 is a sequence of two elementary transform.
            %
            % See also ETS2.plus.
            assert( strcmp(superclasses(ets1), superclasses(ets2)), 'ETS2: both operands must have superclass ETS2, perhaps run ''clear import'', and start over');
            out = [ets1 ets2];
        end

        function out = plus(ets1, ets2)
            %ETS2.plus Compound transforms
            %
            % E1 + E2 is a sequence of two elementary transform.
            %
            % See also ETS2.mtimes.
            assert( strcmp(superclasses(ets1), superclasses(ets2)), 'ETS2: both operands must have superclass ETS2, perhaps run ''clear import'', and start over');

            out = [ets1 ets2];
        end

        function s = structure(ets)
            %ETS2.structure  Show joint type structure
            %
            % E.structure is a character array comprising the letters 'R' or 'P' that
            % indicates the types of joints in the elementary transform sequence E.
            %
            % Notes::
            % - The string will be E.njoints long.
            %
            % See also SerialLink.config.
            s = '';
            for e = ets
                if e.qvar > 0
                    switch e.what
                        case {'Tx', 'Ty'}
                            s = [s 'P'];
                        case 'Rz'
                            s = [s 'R'];
                    end
                end
            end
        end

        function display(ets)
            %ETS2.display Display parameters
            %
            % E.display() displays the transform or transform sequence parameters in
            % compact single line format.
            %
            % Notes::
            % - This method is invoked implicitly at the command line when the result
            %   of an expression is an ETS2 object and the command has no trailing
            %   semicolon.
            %
            % See also ETS2.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');

            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(ets) );
        end % display()

        function s = char(ets)
            %ETS2.char Convert to string
            %
            % E.char() is a string showing transform parameters in a compact format.  If E is a transform sequence (1xN) then
            % the string describes each element in sequence in a single line format.
            %
            % See also ETS2.display.
            s = '';

            function s = render(z)
                if isa(z, 'sym')
                    s = char(z);
                else
                    s = sprintf('%g', z);
                end
            end

            for e = ets
                if e.isjoint
                    s = [s sprintf('%s(q%d)', e.what, e.qvar) ];
                else
                    s = [s sprintf('%s(%s)', e.what, render(e.param))];

                end
            end
        end

        function teach(robot, varargin)
            %ETS2.teach Graphical teach pendant
            %
            % Allow the user to "drive" a graphical robot using a graphical slider
            % panel.
            %
            % ETS.teach(OPTIONS) adds a slider panel to a current ETS plot. If no
            % graphical robot exists one is created in a new window.
            %
            % ETS.teach(Q, OPTIONS) as above but the robot joint angles are set to Q (1xN).
            %
            % Options::
            % 'eul'           Display tool orientation in Euler angles (default)
            % 'rpy'           Display tool orientation in roll/pitch/yaw angles
            % 'approach'      Display tool orientation as approach vector (z-axis)
            % '[no]deg'       Display angles in degrees (default true)
            %
            % GUI::
            % - The Quit (red X) button removes the teach panel from the robot plot.
            %
            % Notes::
            % - The currently displayed robots move as the sliders are adjusted.
            % - The slider limits are derived from the joint limit properties.  If not
            %   set then for
            %   - a revolute joint they are assumed to be [-pi, +pi]
            %   - a prismatic joint they are assumed unknown and an error occurs.
            %
            % See also ETS2.plot.

            %-------------------------------
            % parameters for teach panel
            bgcol = [135 206 250]/255;  % background color
            height = 0.06;  % height of slider rows
            %-------------------------------

            %---- handle options
            opt.deg = true;
            opt.orientation = {'rpy', 'eul', 'approach'};
            opt.d_2d = true;
            opt.callback = [];
            [opt,args] = tb_optparse(opt, varargin);

            if nargin == 1
                q = zeros(1,robot.n);
            else
                q = varargin{1};
                varargin = varargin{2:end};
            end
            robot.plot(q, varargin{:})

            RTBPlot.install_teach_panel('ETS2', robot, q, opt);
        end

        function plot(ets, qq, varargin)
            %ETS2.plot Graphical display and animation
            %
            % ETS.plot(Q, options) displays a graphical animation of a robot based on
            % the transform sequence.  Constant translations are represented as pipe segments, rotational joints as cylinder, and
            % prismatic joints as boxes. The robot is displayed at the joint angle Q (1xN), or
            % if a matrix (MxN) it is animated as the robot moves along the M-point trajectory.
            %
            % Options::
            % 'workspace', W    Size of robot 3D workspace, W = [xmn, xmx ymn ymx zmn zmx]
            % 'floorlevel',L    Z-coordinate of floor (default -1)
            %-
            % 'delay',D         Delay betwen frames for animation (s)
            % 'fps',fps         Number of frames per second for display, inverse of 'delay' option
            % '[no]loop'        Loop over the trajectory forever
            % '[no]raise'       Autoraise the figure
            % 'movie',M         Save an animation to the movie M
            % 'trail',L         Draw a line recording the tip path, with line style L
            %-
            % 'scale',S         Annotation scale factor
            % 'zoom',Z          Reduce size of auto-computed workspace by Z, makes
            %                   robot look bigger
            % 'ortho'           Orthographic view
            % 'perspective'     Perspective view (default)
            % 'view',V          Specify view V='x', 'y', 'top' or [az el] for side elevations,
            %                   plan view, or general view by azimuth and elevation
            %                   angle.
            % 'top'             View from the top.
            %-
            % '[no]shading'     Enable Gouraud shading (default true)
            % 'lightpos',L      Position of the light source (default [0 0 20])
            % '[no]name'        Display the robot's name
            %-
            % '[no]wrist'       Enable display of wrist coordinate frame
            % 'xyz'             Wrist axis label is XYZ
            % 'noa'             Wrist axis label is NOA
            % '[no]arrow'       Display wrist frame with 3D arrows
            %-
            % '[no]tiles'       Enable tiled floor (default true)
            % 'tilesize',S      Side length of square tiles on the floor (default 0.2)
            % 'tile1color',C   Color of even tiles [r g b] (default [0.5 1 0.5]  light green)
            % 'tile2color',C   Color of odd tiles [r g b] (default [1 1 1] white)
            %-
            % '[no]shadow'      Enable display of shadow (default true)
            % 'shadowcolor',C   Colorspec of shadow, [r g b]
            % 'shadowwidth',W   Width of shadow line (default 6)
            %-
            % '[no]jaxes'       Enable display of joint axes (default false)
            % '[no]jvec'        Enable display of joint axis vectors (default false)
            % '[no]joints'      Enable display of joints
            % 'jointcolor',C    Colorspec for joint cylinders (default [0.7 0 0])
            % 'jointcolor',C    Colorspec for joint cylinders (default [0.7 0 0])
            % 'jointdiam',D     Diameter of joint cylinder in scale units (default 5)
            %-
            % 'linkcolor',C     Colorspec of links (default 'b')
            %-
            % '[no]base'        Enable display of base 'pedestal'
            % 'basecolor',C     Color of base (default 'k')
            % 'basewidth',W     Width of base (default 3)
            %
            % The options come from 3 sources and are processed in order:
            % - Cell array of options returned by the function PLOTBOTOPT (if it exists)
            % - Cell array of options given by the 'plotopt' option when creating the
            %   SerialLink object.
            % - List of arguments in the command line.
            %
            % Many boolean options can be enabled or disabled with the 'no' prefix.  The
            % various option sources can toggle an option, the last value encountered is used.
            %
            % Graphical annotations and options::
            %
            % The robot is displayed as a basic stick figure robot with annotations
            % such as:
            % - shadow on the floor
            % - XYZ wrist axes and labels
            % - joint cylinders and axes
            % which are controlled by options.
            %
            % The size of the annotations is determined using a simple heuristic from
            % the workspace dimensions.  This dimension can be changed by setting the
            % multiplicative scale factor using the 'mag' option.
            %
            % Figure behaviour::
            %
            % - If no figure exists one will be created and the robot drawn in it.
            % - If no robot of this name is currently displayed then a robot will
            %   be drawn in the current figure.  If hold is enabled (hold on) then the
            %   robot will be added to the current figure.
            % - If the robot already exists then that graphical model will be found
            %   and moved.
            %
            %
            % Notes::
            % - The options are processed when the figure is first drawn, to make different options come
            %   into effect it is neccessary to clear the figure.
            % - Delay betwen frames can be eliminated by setting option 'delay', 0 or
            %   'fps', Inf.
            % - The size of the plot volume is determined by a heuristic for an all-revolute
            %   robot.  If a prismatic joint is present the 'workspace' option is
            %   required.  The 'zoom' option can reduce the size of this workspace.
            %
            % See also ETS2.teach, SerialLink.plot3d.

            % heuristic to figure robot size
            reach = 0;
            for e=ets
                switch e.what
                    case {'Tx', 'Ty', 'Tz'}
                        if isjoint(e)
                            reach = reach + e.qlim(2);
                        else
                            reach = reach + e.param;
                        end
                end
            end

            opt = RTBPlot.plot_options([], [varargin 'reach', 3, 'top']);
            h = draw_ets(ets, qq, opt);

            set(gca, 'Tag', 'RTB.plot');
            set(gcf, 'Units', 'Normalized');
            pf = get(gcf, 'Position');

            if opt.raise
                % note this is a very time consuming operation
                figure(gcf);
            end

            if strcmp(opt.projection, 'perspective')
                set(gca, 'Projection', 'perspective');
            end

            if isstr(opt.view)
                switch opt.view
                    case 'top'
                        view(0, 90);
                    case 'x'
                        view(0, 0);
                    case 'y'
                        view(90, 0)
                    otherwise
                        error('rtb:plot:badarg', 'view must be: x, y, top')
                end
            elseif isnumeric(opt.view) && length(opt.view) == 2
                view(opt.view)
            end

            % enable mouse-based 3D rotation
            rotate3d on

            ets.animate(qq);
        end

        function animate(ets, qq)
            handles = findobj('Tag', 'ETS2');
            h = handles.UserData;
            opt = h.opt;

            ets = h.ets;
            for q = qq'
                for i=1:length(ets)

                    % create the transform for displaying this element (joint cylinder + link)
                    e = ets(i);
                    if i == 1
                        T = SE2;
                    else
                        T = ets.fkine(q, i-1, 'setopt', opt);
                    end

                    % update the pose of the corresponding graphical element (joint cylinder + link)
                    set(h.element(i), 'Matrix', T.SE3.T);

                    if isprismatic(e)
                        % for prismatic joints, scale the box
                        switch e.what
                            case 'Tx'
                                set(h.pjoint(e.qvar), 'Matrix', diag([q(e.qvar) 1 1 1]));
                            case 'Ty'
                                set(h.pjoint(e.qvar), 'Matrix', diag([1 q(e.qvar) 1 1]));
                        end
                    end
                end

                % update the wrist frame
                T = ets.fkine(q, 'setopt', opt);
                if ~isempty(h.wrist)
                    trplot2(T, 'handle', h.wrist);
                end

                % render and pause
                if opt.delay > 0
                    pause(opt.delay);
                    drawnow
                end
            end
        end

   end

   methods (Access=private)

            function h_ = draw_ets(ets, q, opt)

                clf
                disp('creating new ETS plot');

                axis(opt.workspace);

                s = opt.scale;
                % create an axis
                ish = ishold();
                if ~ishold
                    % if hold is off, set the axis dimensions
                    axis(opt.workspace);
                    hold on
                end

                group = hggroup('Tag', 'ETS2');
                h.group = group;

                % create the graphical joint and link elements
                for i=1:length(ets)
                    e = ets(i);

                    if opt.debug
                        fprintf('create graphics for %s\n', e.char );
                    end

                    % create a graphical depiction of the transform element
                    % This is drawn to resemble orthogonal plumbing.

                    if i == 1
                        T = SE2;
                    else
                        T = ets.fkine(q, i-1, 'setopt', opt);
                    end

                    % create the transform for displaying this element (joint cylinder + link)
                    h.element(i) = hgtransform('Tag', sprintf('element%d', i), 'Matrix', T.SE3.T, 'Parent', h.group);

                    if isjoint(e)
                        % it's a joint element: revolute or prismatic
                        switch e.what
                            case 'Tx'
                                h.pjoint(e.qvar) = hgtransform('Tag', 'prismatic', 'Parent', h.element(i), 'Matrix', diag([q(e.qvar) 1 1 1]));
                                RTBPlot.box('x', opt.jointdiam*s, [0 1], opt.pjointcolor, [], 'Parent', h.pjoint(i));
                            case 'Ty'
                                h.pjoint(e.qvar) = hgtransform('Tag', 'prismatic', 'Parent', h.element(i), 'Matrix', diag([1 q(e.qvar) 1 1]));
                                RTBPlot.box('y', opt.jointdiam*s, [0 1], opt.pjointcolor, [], 'Parent', h.pjoint(i));
                            case 'Rz'
                                RTBPlot.cyl('z', opt.jointdiam*s, opt.jointlen*s*[-1 1], opt.jointcolor, [], 'Parent', h.element(i));
                        end
                    else
                        % it's a constant transform
                        switch e.what
                            case 'Tx'
                                RTBPlot.cyl('x', s, [0 e.param], opt.linkcolor, [], 'Parent', h.element(i));
                            case 'Ty'
                                RTBPlot.cyl('y', s, [0 e.param], opt.linkcolor, [], 'Parent', h.element(i));
                            case 'Rz'
                                RTBPlot.cyl('z', opt.jointdiam*s, opt.jointlen*s*[-1 1], opt.linkcolor, [], 'Parent', h.element(i));
                        end
                    end

                    assert( ~(opt.jaxes && opt.jvec), 'RTB:ETS2:plot:badopt', 'Can''t specify ''jaxes'' and ''jvec''')

                    % create the joint axis line
                    if opt.jaxes
                        if e.isjoint
                            line('XData', [0 0], ...
                                'YData', [0 0], ...
                                'ZData', 14*s*[-1 1], ...
                                'LineStyle', ':', 'Parent', h.element(i));

                            % create the joint axis label
                            text(0, 0, 14*s, sprintf('q%d', e.qvar), 'Parent', h.element(i))
                        end
                    end

                    % create the joint axis vector
                    if opt.jvec
                        if e.isjoint
                            daspect([1 1 1]);
                            ha = arrow3([0 0 -12*s], [0 0 15*s], 'c');
                            set(ha, 'Parent', h.element(i));

                            % create the joint axis label
                            text(0, 0, 20*s, sprintf('q%d', e.qvar), 'Parent', h.element(i))
                        end
                    end

                end

                % display the wrist coordinate frame
                if opt.wrist
                    if opt.arrow
                        % compute arrow3 scale factor...
                        d = axis(gca);
                        if length(d) == 4
                            d = norm( d(3:4)-d(1:2) ) / 72;
                        else
                            d = norm( d(4:6)-d(1:3) ) / 72;
                        end
                        extra = {'arrow', 'width', 1.5*s/d};
                    else
                        extra = {};
                    end
                    h.wrist = trplot2(eye(3,3), 'labels', upper(opt.wristlabel), ...
                        'color', 'k', 'length', opt.wristlen*s, extra{:});
                else
                    h.wrist = [];
                end

                xlabel('X')
                ylabel('Y')
                zlabel('Z')
                grid on

                % restore hold setting
                if ~ish
                    hold off
                end

                h.opt = opt;
                h.ets = ets;

                if nargout > 0
                    h_ = h;
                end

                % attach the handle structure to the top graphical element

                h.q = q;
                handles.opt = opt;

                set(group, 'UserData', h);
            end

         end

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function jFrame = gui_edit_channel( ChannelFile )
% GUI_EDIT_CHANNEL: Load, edit and save a channel file.
%
% USAGE: gui_edit_channel( ChannelFile )

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2018

global GlobalData;

% Progress bar
bst_progress('start', 'Channel editor', 'Initialization...');
% Check if a "Channel Editor" panel already exists
panel = bst_get('Panel', 'ChannelEditor');
% If panel exists
if ~isempty(panel)
    % Close it
    gui_hide(panel);
end
% Create new panel
bstPanel = panel_channel_editor('CreatePanel', ChannelFile);
if isempty(bstPanel)
    bst_progress('stop');
    return;
end
% Show panel in a Java window
bstContainer = gui_show(bstPanel, 'JavaWindow', ['Channel editor: ' ChannelFile], [], 0, 0, 1);
if isempty(bstContainer)
    return;
end

% Get current layout
[jBstArea, FigArea, nScreens, jFigArea, jInsets] = gui_layout('GetScreenBrainstormAreas');
% Place this figure on the screen to take all the possible figure space
jFrame = bstContainer.handle{1};
jFrame.setLocation(jFigArea.getX(), jFigArea.getY() - jInsets.bottom);
jFrame.setSize(jFigArea.getWidth(), jFigArea.getHeight());

% Set selected sensors
if ~isempty(GlobalData.DataViewer.SelectedRows)
    panel_channel_editor('SetChannelSelection', GlobalData.DataViewer.SelectedRows, 0);
end
% Progress bar
bst_progress('stop');

author: brainstorm-tools


text:
function [g1, g2] = lfmaXlfmKernGradient(lfmKern1, lfmKern2, t1, t2, covGrad, meanVector)

% LFMAXLFMKERNGRADIENT Compute a cross gradient between a LFMA and a LFM.
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two lfm kernels, one lfm kernel corresponds to the accel. and
% the other corresponds to the position. It is supposed to be used together
% with the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel (acceleration).
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel (position).
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two lfm kernels, one lfm kernel corresponds to the accel. and
% the other corresponds to the position. It is supposed to be used together
% with the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel (acceleration).
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel (position).
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two lfm kernels, one lfm kernel corresponds to the accel. and
% the other corresponds to the position. It is supposed to be used together
% with the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel (acceleration).
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel (position).
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% ARG meanVec : precomputed factor that is used for the switching dynamical
% latent force model.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% SEEALSO : multiKernParamInit, multiKernCompute, lfmvKernParamInit
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

subComponent = false; % This is just a flag that indicates if this kernel is part of a bigger kernel (SDLFM)

if nargin == 4
    covGrad = t2;
    t2 = t1;
elseif nargin == 6
    subComponent = true;
    if numel(meanVector)>1
        if size(meanVector,1) == 1,
            if size(meanVector, 2)~=size(covGrad, 2)
                error('The dimensions of meanVector don''t correspond to the dimensions of covGrad')
            end
        else
            if size((meanVector'), 2)~=size(covGrad,2)
                error('The dimensions of meanVector don''t correspond to the dimensions of covGrad')
            end
        end
    else
        if numel(t1)==1 && numel(t2)>1
            % matGrad will be row vector and so should be covGrad
            dimcovGrad = length(covGrad);
            covGrad = reshape(covGrad, [1 dimcovGrad]);
        elseif numel(t1)>1 && numel(t2)==1
            % matGrad will be column vector and sp should be covGrad
            dimcovGrad = length(covGrad);
            covGrad = reshape(covGrad, [dimcovGrad 1]);
        end
    end
end

if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end
if lfmKern1.inverseWidth ~= lfmKern2.inverseWidth
    error('Kernels cannot be cross combined if they have different inverse widths.')
end

% Parameters of the simulation (in the order providen by kernExtractParam)

m = [lfmKern1.mass lfmKern2.mass];                  % Par. 1
D = [lfmKern1.spring lfmKern2.spring];              % Par. 2
C = [lfmKern1.damper lfmKern2.damper];              % Par. 3
sigma2 = 2/lfmKern1.inverseWidth;                   % Par. 4
sigma = sqrt(sigma2);
S = [lfmKern1.sensitivity lfmKern2.sensitivity];    % Par. 5

alpha = C./(2*m);
omega = sqrt(D./m-alpha.^2);

% Initialization of vectors and matrices

g1 = zeros(1,5);
g2 = zeros(1,5);

% Precomputations
computeH  = cell(4,1);
computeUpsilonMatrix = cell(2,1);
computeUpsilonVector = cell(2,1);
gradientUpsilonMatrix = cell(4,1);
gradientUpsilonVector = cell(4,1);
gamma1_p = alpha(1) + j*omega(1);
gamma1_m = alpha(1) - j*omega(1);
gamma2_p = alpha(2) + j*omega(2);
gamma2_m = alpha(2) - j*omega(2);
preExp1 = zeros(length(t1),2);
preExp2 = zeros(length(t2),2);
preExpg1 = zeros(length(t1),2);
preExpgg1 = zeros(length(t1),2);
preExpt1 = zeros(length(t1),2);
preExpt2 = zeros(length(t2),2);
preGamma(1) = gamma1_p + gamma2_p;
preGamma(2) = gamma1_p + gamma2_m;
preGamma(3) = gamma1_m + gamma2_p;
preGamma(4) = gamma1_m + gamma2_m;
preGamma2 = preGamma.^2;
preConst = 1./preGamma;
preConst2 = 1./(preGamma2);
preFactors(1) = preConst(2) - preConst(1);
preFactors(2) = preConst(3) - preConst(4);
preFactors(3) = preConst(3) - preConst(1);
preFactors(4) = preConst(2) - preConst(4);
preFactors2(1) = -preConst2(2) + preConst2(1);
preFactors2(2) = -preConst2(3) + preConst2(4);
preFactors2(3) = -preConst2(3) + preConst2(1);
preFactors2(4) = -preConst2(2) + preConst2(4);
preExp1(:,1) = exp(-gamma1_p*t1);
preExp1(:,2) = exp(-gamma1_m*t1);
preExp2(:,1) = exp(-gamma2_p*t2);
preExp2(:,2) = exp(-gamma2_m*t2);
preExpg1(:,1) = (2*gamma1_p)*preExp1(:,1);
preExpg1(:,2) = (2*gamma1_m)*preExp1(:,2);
preExpgg1(:,1) = (gamma1_p^2)*preExp1(:,1);
preExpgg1(:,2) = (gamma1_m^2)*preExp1(:,2);
preExpt1(:,1) = t1.*preExpgg1(:,1);
preExpt1(:,2) = t1.*preExpgg1(:,2);
preExpt2(:,1) = t2.*exp(-gamma2_p*t2);
preExpt2(:,2) = t2.*exp(-gamma2_m*t2);
[computeH{1}, computeUpsilonMatrix{1}, computeUpsilonMatrixLocal{1}] =  lfmComputeH3AP(gamma1_p, gamma1_m, sigma2, t1,t2,preFactors([1 2]), 0);
[computeH{2}, computeUpsilonMatrix{2}, computeUpsilonMatrixLocal{2}] =  lfmComputeH3AP(gamma2_p, gamma2_m, sigma2, t2,t1,preFactors([3 4]), 1);
[computeH{3}, computeUpsilonVector{1}, computeUpsilonVectorLocal{1}] =  lfmComputeH4AP(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExp2, 0 );
[computeH{4}, computeUpsilonVector{2}, computeUpsilonVectorLocal{2}] =  lfmComputeH4AP(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExpgg1, 1 );
preKernel = ( computeH{1} + computeH{2}.' + computeH{3} + computeH{4}.');

% Precompute derivatives
gradientUpsilonMatrix{1} = lfmapGradientUpsilonMatrix(gamma1_p,sigma2,t1, t2, 0, computeUpsilonMatrixLocal{1}{1});
gradientUpsilonMatrix{2} = lfmapGradientUpsilonMatrix(gamma1_m,sigma2,t1, t2, 0, computeUpsilonMatrixLocal{1}{2});
gradientUpsilonMatrix{3} = lfmapGradientUpsilonMatrix(gamma2_p,sigma2,t2, t1, 1, computeUpsilonMatrixLocal{2}{1});
gradientUpsilonMatrix{4} = lfmapGradientUpsilonMatrix(gamma2_m,sigma2,t2, t1, 1, computeUpsilonMatrixLocal{2}{2});
gradientUpsilonVector{1} = lfmapGradientUpsilonVector(gamma1_p,sigma2,t1, computeUpsilonVectorLocal{1}{1});
gradientUpsilonVector{2} = lfmapGradientUpsilonVector(gamma1_m,sigma2,t1, computeUpsilonVectorLocal{1}{2});
gradientUpsilonVector{3} = lfmGradientUpsilonVector(gamma2_p,sigma2,t2);
gradientUpsilonVector{4} = lfmGradientUpsilonVector(gamma2_m,sigma2,t2);

if lfmKern1.isNormalised
    K0 =  lfmKern1.sensitivity*lfmKern2.sensitivity/(8*sqrt(2)*lfmKern1.mass*lfmKern2.mass*prod(omega));
    K02 = 1/(8*sqrt(2)*prod(m)*prod(omega));
else
    K0 =  sigma*sqrt(pi)*lfmKern1.sensitivity*lfmKern2.sensitivity/(8*lfmKern1.mass*lfmKern2.mass*prod(omega));
    K02 = sigma*sqrt(pi)/(8*prod(m)*prod(omega));
end

% Gradient with respect to m, D and C
for ind_theta = 1:3 % Parameter (m, D or C)
    for ind_par = 0:1 % System (1 or 2)
        % Choosing the right gradients for m, omega, gamma1 and gamma2
        switch ind_theta
            case 1  % Gradient wrt m
                gradThetaM = [1-ind_par ind_par];
                gradThetaAlpha = -C./(2*(m.^2));
                gradThetaOmega = (C.^2-2*m.*D)./(2*(m.^2).*sqrt(4*m.*D-C.^2));
            case 2  % Gradient wrt D
                gradThetaM = zeros(1,2);
                gradThetaAlpha = zeros(1,2);
                gradThetaOmega = 1./sqrt(4*m.*D-C.^2);
            case 3  % Gradient wrt C
                gradThetaM = zeros(1,2);
                gradThetaAlpha = 1./(2*m);
                gradThetaOmega = -C./(2*m.*sqrt(4*m.*D-C.^2));
        end
        gradThetaGamma1 = gradThetaAlpha + j*gradThetaOmega;
        gradThetaGamma2 = gradThetaAlpha - j*gradThetaOmega;
        % Gradient evaluation
        gradThetaGamma11 = [gradThetaGamma1(1) gradThetaGamma2(1)];
        gradThetaGamma2 = [gradThetaGamma1(2) gradThetaGamma2(2)];
        gradThetaGamma1 = gradThetaGamma11;

        if ~ind_par % ind_par = k or d
            matGrad = K0 ...
                * ( lfmGradientH31( preFactors([1 2]), preFactors2([1 2]), gradThetaGamma1, ...
                gradientUpsilonMatrix{1}, gradientUpsilonMatrix{2}, computeUpsilonMatrix{1}{1}, ...
                computeUpsilonMatrix{1}{2}, 1) + ...
                lfmGradientH32( preGamma2,  gradThetaGamma1, computeUpsilonMatrix{2}{1}, ...
                computeUpsilonMatrix{2}{2}, 1).' + ...
                lfmGradientH41( preGamma, preGamma2, gradThetaGamma1, preExp2, ...
                gradientUpsilonVector{1}, gradientUpsilonVector{2}, computeUpsilonVector{1}{1},...
                computeUpsilonVector{1}{2}, 1) + ...
                lfmGradientH42AP(preGamma, preGamma2, gradThetaGamma1, preExpg1, preExpgg1, preExpt1, ...
                computeUpsilonVector{2}{1}, computeUpsilonVector{2}{2}).'...
                - (gradThetaM(1+ind_par)/m(1+ind_par) ...
                + gradThetaOmega(1+ind_par)/omega(1+ind_par)) ...
                *preKernel);
        else % ind_par = r or d'
            matGrad = K0 ...
                * ( lfmGradientH31( preFactors([3 4]), preFactors2([3 4]), gradThetaGamma2, ...
                gradientUpsilonMatrix{3}, gradientUpsilonMatrix{4}, computeUpsilonMatrix{2}{1}, ...
                computeUpsilonMatrix{2}{2}, 1).' + ...
                lfmGradientH32( preGamma2([1 3 2 4]),  gradThetaGamma2, computeUpsilonMatrix{1}{1}, ...
                computeUpsilonMatrix{1}{2}, 1) + ...
                lfmGradientH41( preGamma([1 3 2 4]), preGamma2([1 3 2 4]), gradThetaGamma2, preExpgg1, ...
                    gradientUpsilonVector{3}, gradientUpsilonVector{4}, computeUpsilonVector{2}{1},...
                    computeUpsilonVector{2}{2}, 1).' + ...
                lfmGradientH42(preGamma([1 3 2 4]), preGamma2([1 3 2 4]), gradThetaGamma2, preExp2, preExpt2, ...
                computeUpsilonVector{1}{1}, computeUpsilonVector{1}{2}, 1)...
                - (gradThetaM(1+ind_par)/m(1+ind_par) ...
                + gradThetaOmega(1+ind_par)/omega(1+ind_par)) ...
                *preKernel);
        end

        if subComponent
            if size(meanVector,1) ==1,
                matGrad = matGrad*meanVector;
            else
                matGrad = (meanVector*matGrad).';
            end
        end

        % Check the parameter to assign the derivative
        if ind_par == 0
            g1(ind_theta) = sum(sum(matGrad.*covGrad));
        else
            g2(ind_theta) = sum(sum(matGrad.*covGrad));
        end
    end
end

% Gradients with respect to sigma
if lfmKern1.isNormalised
    matGrad = K0*(lfmGradientSigmaH3AP(gamma1_p, gamma1_m, sigma2, t1, t2, preFactors([1 2]), 0)...
        +  lfmGradientSigmaH3AP(gamma2_p, gamma2_m, sigma2, t2, t1, preFactors([3 4]), 1).'...
        +  lfmGradientSigmaH4AP(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExp2, 0 )...
        +  lfmGradientSigmaH4AP(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExpgg1, 1 ).');
else
    matGrad = (prod(S)*sqrt(pi)/(8*prod(m)*prod(omega))) ...
        * (preKernel ...
        + sigma*(lfmGradientSigmaH3AP(gamma1_p, gamma1_m, sigma2, t1, t2, preFactors([1 2]), 0)...
        +  lfmGradientSigmaH3AP(gamma2_p, gamma2_m, sigma2, t2, t1, preFactors([3 4]), 1).'...
        +  lfmGradientSigmaH4AP(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExp2, 0 )...
        +  lfmGradientSigmaH4AP(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExpgg1, 1 ).' ));
end

if subComponent
    if size(meanVector,1) ==1,
        matGrad = matGrad*meanVector;
    else
        matGrad = (meanVector*matGrad).';
    end
end
g1(4) = sum(sum(matGrad.*covGrad))*(-(sigma^3)/4);
g2(4) = g1(4);

% Gradients with respect to S

matGrad = K02*preKernel;

if subComponent
    if size(meanVector,1) ==1,
        matGrad = matGrad*meanVector;
    else
        matGrad = (meanVector*matGrad).';
    end
end

g1(5) = sum(sum(S(2)*matGrad.*covGrad));
g2(5) = sum(sum(S(1)*matGrad.*covGrad));

g2(4) = 0; % Otherwise is counted twice, temporarly changed by Mauricio Alvarez

g1 = real(g1);
g2 = real(g2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [tempimg ]=displayout(Polyg,w,h,img)

r=[255 0 0 0 255 255 0];
g=[0 255 0 255 255 0 0];
b=[0 0 255 255 0 255 0];

rl=[0  0 0 255 255 255 0];
gl=[255 255 0 0 255 0 0];
bl=[0 255 255 0 0 255 0];

fields=zeros(h,w);
%1 floor
polyg=Polyg{1};
if numel(polyg)>0
    tempimg1=poly2mask(polyg(:,1),polyg(:,2),h,w);
else
    tempimg1=zeros(h,w);
end
fields=fields.*(fields~=0)+1*tempimg1.*(fields==0);

%2 middlewall

polyg=Polyg{2};
if numel(polyg)>0
    tempimg1=poly2mask(polyg(:,1),polyg(:,2),h,w);
else
    tempimg1=zeros(h,w);
end
fields=fields.*(fields~=0)+2*tempimg1.*(fields==0);
%3 right wall

polyg=Polyg{3};
if numel(polyg)>0
    tempimg1=poly2mask(polyg(:,1),polyg(:,2),h,w);
else
    tempimg1=zeros(h,w);
end
fields=fields.*(fields~=0)+3*tempimg1.*(fields==0);
%4 left wall

polyg=Polyg{4};
if numel(polyg)>0
    tempimg1=poly2mask(polyg(:,1),polyg(:,2),h,w);
else
    tempimg1=zeros(h,w);
end
fields=fields.*(fields~=0)+4*tempimg1.*(fields==0);
%
%5 ceiling

polyg=Polyg{5};
if numel(polyg)>0
    tempimg1=poly2mask(polyg(:,1),polyg(:,2),h,w);
else
    tempimg1=zeros(h,w);
end
fields=fields.*(fields~=0)+5*tempimg1.*(fields==0);

fields = fields.*(fields~=0) + 6*(fields==0);

mask_r = r(fields);
mask_g = g(fields);
mask_b = b(fields);
mask_color(:,:,1) = mask_r;
mask_color(:,:,2) = mask_g;
mask_color(:,:,3) = mask_b;

tempimg = double(img)*0.5 + mask_color*0.5;
% figure;
% imshow(uint8(tempimg),[]);

return

author: zouchuhang

repo: LayoutNet


text:
function [err,elemErr] = getL2errorRT0(node,elem,exactSigma,sigmah,d)
%% GETL2ERRORRT0 L2 norm of RT0 element.
%
%  The input exactSigma can be a function bundle or vector array.
%  When exactSigma is a vector array, it should be size NT*2, and
%  exactSigma = \nabla u (or exactSigma= K\nabla u), vector array of size
%  NT by 2.
%
%  [err,elemErr] = getL2errorRT0(node,elem,exactSigma,sigmah).
%  err gives the L2 norm between exact flux and approximate one from RT0,
%  elemErr, an NT by 1 array, gives the square of L2 norm between exact
%  flux and approximate one from RT0 on each element.
%
%  [err,elemErr] = getL2errorRT0(node,elem,exactSigma,sigmah,d).
%  If the coefficient 'd' is inlcuded, then the input should be in the way
%  sigmah     = K\nabla u_h  (including K),
%  exactSigma =  \nabla u    (function bundle, don't including K),    or
%  exactSigma = K\nabla u    (vector array, including K).
%  err=||K exactSigma - sigma_h||_{K^(-1)} = |||u-u_h||| is the energy norm.
%  elemErr, an NT by 1 array, gives the square of L2 norm between exact
%  flux and approximate one from RT0 on each element
%
% Example
%
%     maxIt = 5;
%     node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
%     elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
%     bdEdge = setboundary(node,elem,'Dirichlet');
%     pde = mixBCdata;
%     err = zeros(maxIt,1); N = zeros(maxIt,1);
%     for i =1:maxIt
%         [node,elem,bdEdge] = uniformrefine(node,elem,bdEdge);
%         [u,sigma,NULL] = PoissonRT0(node,elem,pde,bdEdge);
%         err(i) = getL2errorRT0(node,elem,pde.Du,sigma);
%         N(i) = size(u,1);
%     end
%     r1 = showrate(N,err,2);
%     legend('||\sigma - \sigma_h||',['N^{' num2str(r1) '}'],...
%            'LOCATION','Best');
%
% See also getHdiverrorRT0, getL2error3RT0.
%
% Created by Ming Wang at Jan 17, 2011, M-lint modified at May 15, 2011.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if nargin >=5 && ~isempty(d)
    if isreal(d)
        K = d;                  % d is an array
    else
    center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
              node(elem(:,3),:))/3;
    K = d(center);              % d is a function
    end
else
    K = [];
end

%% Construct Data Structure
[elem2dof,~,elem2edgeSign] = dofedge(elem);
NT = size(elem,1);% Ndof = max(elem2dof(:)); %N = size(node,1);
[Dlambda,area] = gradbasis(node,elem);
locEdge = [2 3; 3 1; 1 2];

%% Compute square of the L2 error element-wise
[lambda,w] = quadpts(3); % quadrature order is 3
nQuad = size(lambda,1);
err = zeros(NT,1);
rotMat = [0 -1; 1 0]; % rotation matrix for computing rotLambda.
for p = 1:nQuad
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:);
    if ~isnumeric(exactSigma)
        sigmap = exactSigma(pxy);
        if(nargin >=5)&& ~isempty(K) % multiply coeff. if fun bundle flux.
            sigmap = repmat(K,1,2).*sigmap;
        end
    else
        sigmap = exactSigma;         %  K*\nabla u_h.
    end
    sigmahp = zeros(NT,2);
    for k = 1:3 % for each basis
        i = locEdge(k,1); j = locEdge(k,2);
        % phi_k = lambda_iRot_j - lambda_jRot_i;
        sigmahp = sigmahp + repmat(elem2edgeSign(:,k).*sigmah(elem2dof(:,k)),1,2).*...
                   (lambda(p,i)*Dlambda(:,:,j)*rotMat-lambda(p,j)*Dlambda(:,:,i)*rotMat);
    end
    err = err + w(p)*sum((sigmap - sigmahp).^2,2);
end
err = err.*area;               % ||sigma - K\nabla u_h||^2
if(nargin >=5)&& ~isempty(K)   % ||sigma - K\nabla u_h||^2_{K^(-1)}
    err = err./K;
end
elemErr = err;           % ||sigma - K\nabla u_h||^2_{K^(-1)}
% modify the error
err(isnan(err)) = 0;
err = sqrt(sum(err));

repo: Aero-matlab


text:
function [IPPEPoses,refinedPoses] = perspectiveIPPE(U,Q,hEstMethod,opts)
%perspectiveIPPE: The solution to Perspective IPPE with point correspondences computed
%between points in world coordinates on the plane z=0, and normalised points in the
%camera's image.
%
%Inputs:
%
%U: 2xN or 3xN matrix holding the model points in world coordinates. If U
%is 2xN then the points are defined in world coordinates on the plane z=0
%
%Q: 2xN matrix holding the points in the image. These are in normalised
%pixel coordinates. That is, the effects of the camera's intrinsic matrix
%and lens distortion are corrected, so that the Q projects with a perfect
%pinhole model.
%
%hEstMethod: the homography estimation method, either 'Harker' (default) or 'DLT'. If
%'Harker' then the method of Harker and O'Leary is used, from the paper
%"Computation of Homographies" in the 2005 British Computer Vision
%Conference. Otherwise the Direct Linear Transform is used, from Peter Kovesi's implementation at http://www.csse.uwa.edu.au/~pk.
%If hEstMethod=[] then the default is used.
%
%opts is an optional datastructure holding the fields:
%setup IPPE options:
%opts.withPoseRefinement = true; %If this is set, we also refine the camera pose from IPPE by gradient-based optimisation (Levenberg–Marquardt). This is designed to give the statistically optimal pose assuming zero mean I.I.D nose for the feature's 2D positions, but requires initialisation (from IPPE). Note that it about two orders of magnitude slower than IPPE.
%opts.displayTiming = true; %set this to true if you want to see how long IPPE takes in Matlab.
%
%Outputs:
%IPPEPoses: structure containing the two pose solutions from IPPE.
%   IPPEPoses.R1 is the first pose's rotation
%   IPPEPoses.t1 is the first pose's translation

%   IPPEPoses.R2 is the second pose's rotation
%   IPPEPoses.t2 is the second pose's translation

%   IPPEPoses.reprojError1 is the first pose's reprojection error
%   IPPEPoses.reprojError2 is the second pose's reprojection error
%
%refinedPoses: structure containing the two refined pose solutions using
%Levenberg–Marquardt (same format as IPPEPoses). This is only outputted if opts.withPoseRefinement = true
%
%Note that the poses are transforms from world
%coordinates to camera coordinates. If you
%want to get the camera's pose in world coordinates you would get it with
%M1 = inv([[R1,t1];[0,0,0,1]]) and M2 = inv([[R2,t2];[0,0,0,1]])
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of the IPPE package for very fast plane-based pose
% estimation from a single perspective image, from the paper "Infinitesimal Plane-based Pose Estimation"
% by Toby Collins and Adrien Bartoli, published in the International Journal of Computer Vision, September 2014.
% A copy of the author's pre-print version can be found here: http://isit.u-clermont1.fr/~ab/Publications/Collins_Bartoli_IJCV14.pdf
%
% Feel free to contact Toby (toby.collins@gmail.com) if you have any
% questions about the paper and IPPE.
%
% This package is free and covered by the BSD licence without any warranty. We hope you find this code useful and please cite our paper in your work:
% (c) Toby Collins 2015
%
%@article{
%year={2014},
%issn={0920-5691},
%journal={International Journal of Computer Vision},
%volume={109},
%number={3},
%doi={10.1007/s11263-014-0725-5},
%title={Infinitesimal Plane-Based Pose Estimation},
%url={http://dx.doi.org/10.1007/s11263-014-0725-5},
%publisher={Springer US},
%keywords={Plane; Pose; SfM; PnP; Homography},
%author={Collins, Toby and Bartoli, Adrien},
%pages={252-286},
%language={English}
%}

opts = parseOpts(opts);

measureTiming = opts.measureTiming;
refineSolutions = opts.withPoseRefinement;

U = double(U);
Q = double(Q);

%set defaults for hEstMethod:
if isempty(hEstMethod)
    hEstMethod = 'Harker';
else
    if ~(strcmp(hEstMethod,'Harker')||strcmp(hEstMethod,'DLT'))
        error('You must specify a valid homography estimation method in perspectiveIPPE. The options are Harker or DLT');
    end
end

%argument sanity check:
assert(size(U,1)==2|size(U,1)==3);
assert(size(U,2)==size(Q,2));
assert(size(Q,1) == 2);

%n is the number of points:
n = size(U,2);
modelDims = size(U,1);

Uuncentred = U;

if modelDims==2
   %we zero-center the model points:
    Pbar = [mean(U,2);0];
    U(1,:) = U(1,:)-Pbar(1);
    U(2,:) = U(2,:)-Pbar(2);
else
    %we rotate the model points onto the plane z=0 and zero center them:
    Pbar = mean(U,2);
    MCenter = eye(4);
    MCenter(1:3,end) = -Pbar;
    U_ = MCenter(1:3,:)*[U;ones(1,size(U,2))];
    %compute the rotation that rotates the model points onto the plane z=0:
    [modelRotation,sigs,~] = svd(U_*U_'); modelRotation=modelRotation';
    if det(modelRotation)<0
        modelRotation(3,:)  = -modelRotation(3,:); %this ensures modelRotation is a member of SO3
    end
    if (sigs(3,3)/sigs(2,2)>1e-5)
        error('IPPE requires the model points to be planar!');
    end

    modelRotation(4,4) = 1;
    Mcorrective = modelRotation*MCenter;
    U = Mcorrective(1:2,:)*[U;ones(1,size(U,2))];

end

if measureTiming
    tic;
end
%compute the model-to-image homography:
switch hEstMethod
    case 'DLT'
        H = homography2d([U;ones(1,n)],[Q;ones(1,n)]);
    case 'Harker'
        H = homographyHarker([U;ones(1,n)],[Q;ones(1,n)]);
end

%[ px, py, 1, 0, 0, 0, -px*qx, -py*qx]h = qx
%[ 0, 0, 0, px, py, 1, -px*qy, -py*qy]h = qy

%Compute the Jacobian J of the homography at (0,0):
H = H./H(3,3);
J(1,1) = H(1,1)-H(3,1)*H(1,3);
J(1,2) = H(1,2)-H(3,2)*H(1,3);
J(2,1) = H(2,1)-H(3,1)*H(2,3);
J(2,2) = H(2,2)-H(3,2)*H(2,3);

%Compute the two rotation solutions:
v = [H(1,3);H(2,3)];
[R1,R2] = IPPE_dec(v,J);

%Compute the two translation solutions:
t1_ = estT(R1,U,Q);
t2_ = estT(R2,U,Q);

if modelDims==2
    t1 = [R1,t1_]*[-Pbar;1];
    t2 = [R2,t2_]*[-Pbar;1];
else
    M1 = [R1,t1_];
    M1(4,4) = 1;
    M2 = [R2,t2_];
    M2(4,4) = 1;
    M1 = M1*Mcorrective;
    M2 = M2*Mcorrective;
    R1 = M1(1:3,1:3);
    R2 = M2(1:3,1:3);
    t1 = M1(1:3,end);
    t2 = M2(1:3,end);
end

[reprojErr1,reprojErr2] = computeReprojErrs(R1,R2,t1,t2,Uuncentred,Q);
%sort solutions:
if reprojErr1>reprojErr2
    [R1,R2,t1,t2,reprojErr1,reprojErr2] = swapSolutions(R1,R2,t1,t2,reprojErr1,reprojErr2);
end
IPPEPoses.R1 = R1;
IPPEPoses.t1 = t1;
IPPEPoses.R2 = R2;
IPPEPoses.t2 = t2;
IPPEPoses.reprojError1 = reprojErr1;
IPPEPoses.reprojError2 = reprojErr2;

%%with timing?
if measureTiming
    disp('Computing time to estimate homography for IPPE. Warning this slows things down, so turn it off when you are not benchmarking!');
    tic;
    %average timing over 100 runs:
    for i=1:100
        switch hEstMethod
            case 'DLT'
                H = homography2d([U;ones(1,n)],[Q;ones(1,n)]);
            case 'Harker'
                H = homographyHarker([U;ones(1,n)],[Q;ones(1,n)]);
        end
    end

    tm = toc;
    disp(['Time to estimate homography: ' num2str(1000*tm/100) 'ms']);

    disp('Computing time to solve IPPE given homography...');
    tic;
    %average timing over 100 runs:
    for i=1:100
        %Compute the Jacobian J of the homography at (0,0):
        H = H./H(3,3);
        J(1,1) = H(1,1)-H(3,1)*H(1,3);
        J(1,2) = H(1,2)-H(3,2)*H(1,3);
        J(2,1) = H(2,1)-H(3,1)*H(2,3);
        J(2,2) = H(2,2)-H(3,2)*H(2,3);
        %Compute the two rotation solutions:
        v = [H(1,3);H(2,3)];
        [R1,R2] = IPPE_dec(v,J);
        %Compute the two translation solutions:
        t1_ = estT(R1,U,Q);
        t2_ = estT(R2,U,Q);
        t1 = [R1,t1_]*[-Pbar;1];
        t2 = [R2,t2_]*[-Pbar;1];
    end
    tm = toc;
    disp(['Time to solve IPPE given homography: ' num2str(1000*tm/100) 'ms']);
end

if modelDims==2
    Uuncentred(3,:) = 0;
end

if refineSolutions
    poseEst= pnpLMRefinement(Uuncentred,Q,R1,t1,'MLE');
    R1 = poseEst.R;
    t1 = poseEst.t;

    poseEst= pnpLMRefinement(Uuncentred,Q,R2,t2,'MLE');
    R2 = poseEst.R;
    t2 = poseEst.t;

    [reprojErr1,reprojErr2] = computeReprojErrs(R1,R2,t1,t2,Uuncentred,Q);
    %Sort the solutions:
    if reprojErr1>reprojErr2
        [R1,R2,t1,t2,reprojErr1,reprojErr2] = swapSolutions(R1,R2,t1,t2,reprojErr1,reprojErr2);
    end
    refinedPoses.R1 = R1;
    refinedPoses.t1 = t1;
    refinedPoses.R2 = R2;
    refinedPoses.t2 = t2;
    refinedPoses.reprojError1 = reprojErr1;
    refinedPoses.reprojError2 = reprojErr2;

    %%with timing?
    if measureTiming
        disp('Computing time to refine IPPE solutions with Levenberg Marquardt...');
        tic;
        %average timing over 100 runs:
        for i=1:100
            pnpLMRefinement(Uuncentred,Q,IPPEPoses.R1,IPPEPoses.t1,'MLE');
            pnpLMRefinement(Uuncentred,Q,IPPEPoses.R2,IPPEPoses.t2,'MLE');
        end
        tm = toc;
        disp(['Time to refine IPPE solutions with Levenberg Marquardt: ' num2str(1000*tm/100) 'ms']);
    end

end

function [reprojErr1,reprojErr2] = computeReprojErrs(R1,R2,t1,t2,U,Q)
%computeReprojErrs: Computes the reprojection errors for the two solutions
%generated by IPPE.

%transform model points to camera coordinates and project them onto the
%image:
if size(U,1)==2
    PCam1 = R1(:,1:2)*U;
    PCam2 = R2(:,1:2)*U(1:2,:);
else
    PCam1 = R1*U;
    PCam2 = R2*U;
end

PCam1(1,:) = PCam1(1,:) + t1(1);
PCam1(2,:) = PCam1(2,:) + t1(2);
PCam1(3,:) = PCam1(3,:) + t1(3);

PCam2(1,:) = PCam2(1,:) + t2(1);
PCam2(2,:) = PCam2(2,:) + t2(2);
PCam2(3,:) = PCam2(3,:) + t2(3);

Qest_1 = PCam1./[PCam1(3,:);PCam1(3,:);PCam1(3,:)];
Qest_2 = PCam2./[PCam2(3,:);PCam2(3,:);PCam2(3,:)];

%compute reprojection errors:
reprojErr1 = norm(Qest_1(1:2,:)-Q);
reprojErr2 = norm(Qest_2(1:2,:)-Q);

function t = estT(R,psPlane,Q)
%Computes the least squares estimate of translation given the rotation solution.
if size(psPlane,1)==2
    psPlane(3,:) =0;
end
%qq = homoMult(Kinv,pp')';
Ps = R*psPlane;

numPts = size(psPlane,2);
Ax = zeros(numPts,3);
bx = zeros(numPts,1);

Ay = zeros(numPts,3);
by = zeros(numPts,1);

Ax(:,1) = 1;
Ax(:,3) = -Q(1,:);
bx(:) = Q(1,:).*Ps(3,:) -  Ps(1,:);

Ay(:,2) = 1;
Ay(:,3) = -Q(2,:);
by(:) = Q(2,:).*Ps(3,:) -  Ps(2,:);

A = [Ax;Ay];
b = [bx;by];

AtA = A'*A;
Atb = A'*b;

Ainv = IPPE_inv33(AtA);
t = Ainv*Atb;

function [R1,R2,t1,t2,reprojErr1,reprojErr2] = swapSolutions(R1_,R2_,t1_,t2_,reprojErr1_,reprojErr2_)
%swaps the solutions:
R1 = R2_;
t1 = t2_;
reprojErr1 = reprojErr2_;

R2 = R1_;
t2 = t1_;
reprojErr2 = reprojErr1_;

function Ainv = IPPE_inv33(A)
%computes the inverse of a 3x3 matrix, assuming it is full-rank.
a11 = A(1,1);
a12 = A(1,2);
a13 = A(1,3);

a21 = A(2,1);
a22 = A(2,2);
a23 = A(2,3);

a31 = A(3,1);
a32 = A(3,2);
a33 = A(3,3);

Ainv = [[ a22*a33 - a23*a32, a13*a32 - a12*a33, a12*a23 - a13*a22];
    [ a23*a31 - a21*a33, a11*a33 - a13*a31, a13*a21 - a11*a23];
    [ a21*a32 - a22*a31, a12*a31 - a11*a32, a11*a22 - a12*a21]];
Ainv = Ainv./(a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);

function opts = parseOpts(opts)
if ~isfield(opts,'measureTiming')
    opts.measureTiming = false;
end
if ~isfield(opts,'withPoseRefinement')
    opts.withPoseRefinement = true;
end

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
function stats = rmanova2(data,alpha,doplot,ttst)
%
%Repeated-measures two-way ANOVA
%
% :Usage:
% ::
%
%     stats = rmanova2(data,[alpha],[doplot],[ttst]);
%
% :Inputs:
%
%   **data:**
%        can be one of two formats:
%        1. Cell array - each row represents a level of factor 1, and
%           each column represents a level of factor 2. Each cell contains a
%           vector of values of the dependent variable for each subject.
%        2. Matrix - each row represents a trial, with the following
%           columns:
%             - column1 - dependent variable
%             - column2 - grouping variable for subject
%             - column3 - grouping variable for factor 1
%             - column4 - grouping variable for factor 2
%
%   **alpha:**
%        (optional) p-value threshold (default: 0.05)
%
%   **doplot:**
%        (optional) if 1, will produce a line plot.
%                       Works only for cell input data (default: 1)
%
%   **ttst:**
%        (optional) if 1, will perform pairwise t-tests (default: 0)
%
% ..
%    Aaron Schurger (2005.02.04)
%    Derived from Keppel & Wickens (2004) "Design and Analysis" ch. 18
%    Modified by Sam Gershman (2006.11.16)
% ..

if nargin < 2; alpha = 0.05; doplot = 0; ttst = 0; end
if nargin < 3; doplot = 1;  ttst = 0; end
if nargin < 4; ttst = 0; end

if iscell(data);
    Y = []; S = []; F1 = []; F2 = [];
    [m n] = size(data);
    for i = 1:m;
        for j = 1:n;
            Y = cat(1,Y,data{i,j});
            nsubs = length(data{i,j});
            S = cat(2,S,1:nsubs);
            F1 = cat(1,F1,zeros(nsubs,1)+i);
            F2 = cat(1,F2,zeros(nsubs,1)+j);
        end
    end
    S = S';
else
    Y = data(:,1);
    S = data(:,2);
    F1 = data(:,3); F2 = data(:,4);
    clear data
    su = unique(S); f1u = unique(F1); f2u = unique(F2);
    for a = 1:length(f1u);
        for b = 1:length(f2u);
            ii = intersect(find(F1==f1u(a)),find(F2==f2u(b)));
            d = [Y(ii) S(ii)];
            d = sortrows(d,2);
            data{a,b} = d(:,1);
        end
    end
end

F1_lvls = unique(F1);
F2_lvls = unique(F2);
Subjs = unique(S);

a = length(F1_lvls); % # of levels in factor 1
b = length(F2_lvls); % # of levels in factor 2
n = length(Subjs); % # of subjects

INDS = cell(a,b,n); % this will hold arrays of indices
CELLS = cell(a,b,n); % this will hold the data for each subject X condition
MEANS = zeros(a,b,n); % this will hold the means for each subj X condition

% Calculate means for each subject X condition.
% Keep data in CELLS, because in future we may want to allow options for
% how to compute the means (e.g. leaving out outliers > 3stdev, etc...).
for i=1:a % F1
    for j=1:b % F2
        for k=1:n % Subjs
            INDS{i,j,k} = find(F1==F1_lvls(i) & F2==F2_lvls(j) & S==Subjs(k));
            CELLS{i,j,k} = Y(INDS{i,j,k});
            MEANS(i,j,k) = mean(CELLS{i,j,k});
        end
    end
end

% make tables (see table 18.1, p. 402)
AB = reshape(sum(MEANS,3),a,b); % across subjects
AS = reshape(sum(MEANS,2),a,n); % across factor 2
BS = reshape(sum(MEANS,1),b,n); % across factor 1

A = sum(AB,2); % sum across columns, so result is ax1 column vector
B = sum(AB,1); % sum across rows, so result is 1xb row vector
S = sum(AS,1); % sum across columns, so result is 1xs row vector
T = sum(sum(A)); % could sum either A or B or S, choice is arbitrary

% degrees of freedom
dfA = a-1;
dfB = b-1;
dfAB = (a-1)*(b-1);
dfS = n-1;
dfAS = (a-1)*(n-1);
dfBS = (b-1)*(n-1);
dfABS = (a-1)*(b-1)*(n-1);

% bracket terms (expected value)
expA = sum(A.^2)./(b*n);
expB = sum(B.^2)./(a*n);
expAB = sum(sum(AB.^2))./n;
expS = sum(S.^2)./(a*b);
expAS = sum(sum(AS.^2))./b;
expBS = sum(sum(BS.^2))./a;
expY = sum(Y.^2);
expT = T^2 / (a*b*n);

% sums of squares
ssA = expA - expT;
ssB = expB - expT;
ssAB = expAB - expA - expB + expT;
ssS = expS - expT;
ssAS = expAS - expA - expS + expT;
ssBS = expBS - expB - expS + expT;
ssABS = expY - expAB - expAS - expBS + expA + expB + expS - expT;
ssTot = expY - expT;

% mean squares
msA = ssA / dfA;
msB = ssB / dfB;
msAB = ssAB / dfAB;
msS = ssS / dfS;
msAS = ssAS / dfAS;
msBS = ssBS / dfBS;
msABS = ssABS / dfABS;

% f statistic
fA = msA / msAS;
fB = msB / msBS;
fAB = msAB / msABS;

% p values
pA = 1-fcdf(fA,dfA,dfAS);
pB = 1-fcdf(fB,dfB,dfBS);
pAB = 1-fcdf(fAB,dfAB,dfABS);

% return values
 stats.SS1 = ssA;
 stats.df1 = dfA;
 stats.MS1 = msA;
 stats.F1 = fA;
 stats.P1 = pA;
 stats.SS2 = ssB;
 stats.df2 = dfB;
 stats.MS2 = msB;
 stats.F2 = fB;
 stats.P2 = pB;
 stats.SS12 = ssAB;
 stats.df12 = dfAB;
 stats.MS12 = msAB;
 stats.F12 = fAB;
 stats.P12 = pAB;

 stats.alpha = alpha;

 %decision rule
if stats.P1 < alpha; stats.significant1 = 'Yes'; else; stats.significant1 = 'No'; end;
if stats.P2 < alpha; stats.significant2 = 'Yes'; else; stats.significant2 = 'No'; end;
if stats.P12 < alpha; stats.significant12 = 'Yes'; else; stats.significant12 = 'No'; end;

%make line plots
if doplot
    anova_line_plot(data);

%     [m n] = size(data); plotcounter = 0;
%     figure;
%     for i = 1:m;
%         for j = 1:n;
%             plotcounter = plotcounter + 1;
%             subplot(m,n,plotcounter);
%             boxplot(data{i,j});
%             xlabel([num2str(i),',',num2str(j)]);
%             set(gca,'XTickLabel',[]);
%         end
%     end
end

%pairwise t-tests
if ttst;
    disp('Performing post-hoc t-tests...');
    tcounter = 0; a = 0; b = 0;
    [m n] = size(data);
    for t1 = 1:m;
        for t2 = 1:n;
            a = a + 1; b = 0;
            for t3 = 1:m;
                for t4 = 1:n;
                    b = b + 1;
                    if b > a;
                        tcounter = tcounter + 1;
                        [h,p,ci,stat] = ttest(data{t1,t2},data{t3,t4});
                        stats.ttests(tcounter).P = p;
                        stats.ttests(tcounter).comparison = [num2str(t1),',',num2str(t2),' > ',num2str(t3),',',num2str(t4)];
                        stats.ttests(tcounter).means = [mean(data{t1,t2}) mean(data{t3,t4})];
                        stats.ttests(tcounter).ci = ci';
                    end
                end
            end
        end
    end
end

end

function anova_line_plot(anova_dat)

    if ~iscell(anova_dat)
        error('For line plots, enter cells rather than matrix input!');
    end

for i = 1:size(anova_dat, 1)
for j = 1:size(anova_dat, 2)
linedat(i, j) = nanmean(anova_dat{i, j});
lineste(i, j) = ste(anova_dat{i, j});
end
end

% Transpose so F2 is on X-axis, F1 is lines
linedat = linedat';
lineste = lineste';

create_figure('Lines');

plot(linedat, 'o-', 'LineWidth', 3, 'MarkerSize', 12, 'MarkerFaceColor', [.5 .5 .5])

set(gca, 'XLim', [.5 size(anova_dat, 2)+.5], 'XTick', [1:size(anova_dat, 2)]); xlabel('Factor 2');

for j = 1:size(anova_dat, 2), f1names{j} = ['F1: Level ' num2str(j)]; end

legend(f1names)

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function c = tapas_hgf_binary_mab_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF) in a multi-armded bandit
% situation for binary inputs in the absence of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading exactly to the model introduced in
% Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_hgf_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_binary_mab_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.rho        row vector of rhos (representing drift; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.sa          sigma (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels, 3rd dim: bandits)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_hgf_binary_mab_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013-2017 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'hgf_binary_mab';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Number of bandits
c.n_bandits = 3;

% Coupling
% This may only be set to true if c.n_bandits is set to 2 above. If
% true, it means that the two bandits' winning probabilities are
% coupled in the sense that they add to 1 and are both updated on
% each trial even though only the outcome for one of them is observed.
c.coupled = false;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 0, 0];

c.logsa_0mu = [NaN,   log(0.1), log(1)];
c.logsa_0sa = [NaN,          0,      0];

% Rhos
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero to turn off drift.
c.rhomu = [NaN, 0, 0];
c.rhosa = [NaN, 0, 0];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,      0];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,  -2,  -6];
c.omsa = [NaN, 4^2, 4^2];

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.rhomu,...
    c.logkamu,...
    c.ommu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.rhosa,...
    c.logkasa,...
    c.omsa,...
         ];

% Check whether we have the right number of priors
expectedLength = 3*c.n_levels+2*(c.n_levels-1)+1;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_hgf_binary_mab;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_binary_mab_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function A_xyz = computeVelCoeffsMtx(sectorOrder)
%COMPUTEVELCOEFFSMTX Compute the generating matrices of products of
%beampatterns with dipoles
%
%   This routine computes the matrices that generate the coefficients of
%   the beampattern of order (sectorOrder+1) that is essentially the
%   product of a pattern of order=sectorOrder and a dipole. It is used in
%   beamWeightsVelocityPatterns.m. For the derivation of the matrices see:
%
%       Politis, A. and Pulkki, V., 2016.
%       Acoustic intensity, energy-density and diffuseness estimation in a directionally-constrained region.
%       arXiv preprint arXiv:1609.03409.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% COMPUTEVELCOEFFSMTX.M - 7/06/2015
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% orders of sector and resulting velocity patterns
Ns = sectorOrder;
Nxyz = Ns+1;

% Complex SH coefficients of velocity x,y,z
% X = cos\phi sin\theta
% Y = sin\phi sin\theta
% Z = cos\theta
% q = 1 -> Yq = Y_{1(-1)}
% q = 2 -> Yq = Y_{10}
% q = 3 -> Yq = Y_{11}
% x_q != 0 only for q=1,3
% y_q != 0 only for q=1,3
% z_q != 0 only for q=2
x1 = sqrt(2*pi/3);
x3 = -sqrt(2*pi/3);
y1 = 1i*sqrt(2*pi/3);
y3 = 1i*sqrt(2*pi/3);
z2  = sqrt(4*pi/3);

A_xyz = zeros((Nxyz+1)^2, (Ns+1)^2, 3);

G_mtx = gaunt_mtx(Ns, 1, Nxyz);

for ii=1:(Nxyz+1)^2
    for jj=1:(Ns+1)^2

        A_xyz(ii,jj,1) = x1*G_mtx(jj,2,ii) + x3*G_mtx(jj,4,ii);
        A_xyz(ii,jj,2) = y1*G_mtx(jj,2,ii) + y3*G_mtx(jj,4,ii);
        A_xyz(ii,jj,3) = z2*G_mtx(jj,3,ii);
    end
end

author: polarch

repo: Spherical-Array-Processing


text:
function h = times(f, g)
%.*   Times of two BALLFUNV objects.
%   F.*G if F is a BALLFUNV and G is double returns the BALLFUNV after
%   componentwise multiplication.
%
%   F.*G if F is a BALLFUNV and G is a BALLFUNV returns the BALLFUNV
%   after multiplication of F by each component of G.
%
% See also MTIMES.

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

if isempty( f ) || isempty(g)
    h = ballfunv();
    return
end

fIsBallfun = isa(f, 'ballfun');
gIsBallfun = isa(g, 'ballfun');

if fIsBallfun
    G = g.comp;
    h = [f.*G{1} ; f.*G{2} ; f.*G{3}];
elseif gIsBallfun
    F = f.comp;
    h = [F{1}.*g ; F{2}.*g ; F{3}.*g];
else
    F = f.comp;
    G = g.comp;
    h = [F{1}.*G{1} ; F{2}.*G{2} ; F{3}.*G{3}];
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function prtm(hpop)
    % function to print output form GenAS either to printer or postscript file
    % (see call from gogenas.m)
    %

    report_this_filefun(mfilename('fullpath'));

    val = get(hpop,'Value');
    if val == 2
        print
    elseif val == 3
        print -dpsc genas.ps
        disp('  Plot saved as poscript file genas.ps ...');
    end

author: CelsoReyes


text:
function y = geo_mean( x, dim, w )
error( nargchk( 1, 3, nargin ) );

%GEO_MEAN   Internal cvx version.

%
% Basic argument check
%

sx = size( x );
if nargin < 2 || isempty( dim ),
    dim = cvx_default_dimension( sx );
elseif ~cvx_check_dimension( dim ),
    error( 'Second argument must be a positive integer.' );
end

%
% Determine sizes, quick exit for empty arrays
%

sx = [ sx, ones( 1, dim - length( sx ) ) ];
nx = sx( dim );
sy = sx;
sy( dim ) = 1;
if any( sx == 0 ),
    y = ones( sy );
    return
end

%
% Third and fourth argument check
%

if nargin < 3 || isempty( w ),
    w = [];
elseif numel( w ) ~= length( w ) || ~isnumeric( w ) || ~isreal( w ) || any( w < 0 ) || any( w ~= floor( w ) ),
    error( 'Third argument must be a vector of nonnegative integers.' );
elseif length( w ) ~= nx,
    error( 'Third argument must be a vector of length %d', nx );
elseif ~any( w ),
    y = ones( sy );
    return
end

%
% Type check
%

persistent remap_1 remap_2 remap_3 remap_4
if isempty( remap_4 ),
    % Constant (postive or negative)
    remap_1 = cvx_remap( 'real' );
    remap_2 = cvx_remap( 'concave' );
    remap_3 = cvx_remap( 'log-convex' );
    remap_4 = cvx_remap( 'log-concave' );
end
vx = cvx_reshape( cvx_classify( x ), sx );
t1 = all( reshape( remap_1( vx ), sx ), dim );
t2 = all( reshape( remap_2( vx ), sx ), dim );
t3 = all( reshape( remap_3( vx ), sx ), dim ) | ...
     all( reshape( remap_4( vx ), sx ), dim );
% Valid combinations with zero or negative entries can be treated as constants
t1 = t1 | ( ( t2 | t3 ) & any( vx == 1 | vx == 9, dim ) );
ta = t1 + ( 2 * t2 + 3 * t3 ) .* ~t1;
nu = sort( ta(:) );
nu = nu([true;diff(nu)~=0]);
nk = length( nu );

%
% Permute and reshape, if needed
%

perm = [];
if nk > 1 || ( any( nu > 1 ) && nx > 1 ),
    if dim > 1 && any( sx( 1 : dim - 1 ) > 1 ),
        perm = [ dim, 1 : dim - 1, dim + 1 : length( sx ) ];
        x    = permute( x,  perm );
        sx   = sx( perm ); %#ok
        sy   = sy( perm );
        ta   = permute( ta, perm );
        dim  = 1;
    else
        perm = [];
    end
    nv = prod( sy );
    x  = reshape( x, nx, nv );
    ta = reshape( ta, 1, nv );
end

%
% Perform the computations
%

if nk > 1,
    y = cvx( [ 1, nv ], [] );
end
for k = 1 : nk,

    if nk == 1,
        xt = x;
        sz = sy; %#ok
    else
        tt = ta == nu( k );
        xt = cvx_subsref( x, ':', tt );
        nv = nnz( tt );
        sz = [ 1, nv ]; %#ok
    end

    switch nu( k ),
        case 0,
            error( 'Disciplined convex programming error:\n   Invalid computation: geo_mean( {%s} )', cvx_class( xt, true, true ) );
        case 1,
            yt = cvx( geo_mean( cvx_constant( xt ), dim, w ) );
        case 2,
            cvx_begin
                hypograph variable yt(sz);
                { cvx_accept_concave(xt), yt } == geo_mean_cone( size(xt), dim,  w, 'func' );
            cvx_end
        case 3,
            if nx == 1,
                yt = xt;
            elseif isempty( w ),
                yt = exp( sum( log( xt ), 1 ) * ( 1 / nx ) );
            else
                yt = exp( ( w / sum( w ) ) * log( xt ) );
            end
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    if nk == 1,
        y = yt;
    else
        y = cvx_subsasgn( y, tt, yt );
    end

end

%
% Reverse the reshaping and permutation steps
%

y = reshape( y, sy );
if ~isempty( perm ),
    y = ipermute( y, perm );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
% maximally flat delay LC filter pp 323 A. I. Zveref
% Copyright 2004-2013 The MathWorks, Inc.
R=8;      % load resistance
fc=80e3;  % bandpass corner

w=2*pi*fc;

c1=1.5012; l2=0.978; c3=.612;l4=.211;  % Rs=inf,  Rl = 1
% scale for freq and impedance ...
C1=c1/(R*w)
L2=l2*R/w
C3=c3/(R*w)
L4=l4*R/w

author: Sable

repo: mcbench-benchmarks


text:
function [xVals, yVals, zVals] = getUniformizedXYZVals(planC)
%"getUniformizedXYZVals"
%   Return the X Y and Z vals of the uniformized dataset.
%
%   These are arrays, with a value for each row, col and slice in the
%   uniformized dataset.
%
%   In case of any error accessing the data 0 is returned for all vals.
%
% JRA 11/14/03
%
% Usage: [xVals, yVals, zVals] = getUniformizedXYZVals(planC)

indexS = planC{end};

[xVals, yVals, junk] = getScanXYZVals(planC{indexS.scan}(1));

scanInfo = planC{indexS.scan}(1).scanInfo(1);
sizeArray = getUniformizedSize(planC);

uniformScanInfo = planC{indexS.scan}(1).uniformScanInfo;
nZSlices = sizeArray(3);
zVals = uniformScanInfo.firstZValue : uniformScanInfo.sliceThickness : uniformScanInfo.sliceThickness * (nZSlices-1) + uniformScanInfo.firstZValue;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function d =  testing(a,d)

 n= get_dim(d);
 r = zeros(n,1);
 for i=1:a.kmax
 	r =  r + a.alpha(i).*get_x(test(a.child{i},d));
 end

 r=sign(r);
%  if(c > 1)
%  	L=-ones(c)+2*eye(c);
% 	for i=1:n
% 		[mi ind]=min(sum(abs(repmat(r(i,:),c,1)-L),2));
% 		r(i,:)=L(ind,:);
% 	end
%  end

 d=set_x(d,r);
 d=set_name(d,[get_name(d) ' -> ' get_name(a)]);


text:
% [MIN, MAX] = range2(MTX)
%
% Compute minimum and maximum values of MTX, returning them as a 2-vector.

% Eero Simoncelli, 3/97.

function [mn, mx] = range2(mtx)

%% NOTE: THIS CODE IS NOT ACTUALLY USED! (MEX FILE IS CALLED INSTEAD)

fprintf(1,'WARNING: You should compile the MEX code for "range2", found in the MEX subdirectory.  It is MUCH faster.\n');

if (~isreal(mtx))
  error('MTX must be real-valued');
end

mn = min(min(mtx));
mx = max(max(mtx));

author: ricedsp


text:
function [t,u,v,idx]=raytrace(p0,v0,node,face)
%
% [t,u,v,idx]=raytrace(p0,v0,node,face)
%
% perform a Havel-styled ray tracing for a triangular surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   p0: starting point coordinate of the ray
%   v0: directional vector of the ray
%   node: a list of node coordinates (nn x 3)
%   face: a surface mesh triangle list (ne x 3)
%
% output:
%   t: signed distance from p to the intersection point for each surface
%      triangle, if ray is parallel to the triangle, t is set to Inf
%   u: bary-centric coordinate 1 of all intersection points
%   v: bary-centric coordinate 2 of all intersection points
%      the final bary-centric triplet is [u,v,1-u-v]
%   idx: optional output, if requested, idx lists the IDs of the face
%      elements that intersects the ray; users can manually calc idx by
%
%      idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
%
% Reference:
%  [1] J. Havel and A. Herout, "Yet faster ray-triangle intersection (using
%          SSE4)," IEEE Trans. on Visualization and Computer Graphics,
%          16(3):434-438 (2010)
%  [2] Q. Fang, "Comment on 'A study on tetrahedron-based inhomogeneous
%          Monte-Carlo optical simulation'," Biomed. Opt. Express, (in
%          press)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

p0=p0(:)';
v0=v0(:)';

AB=node(face(:,2),1:3)-node(face(:,1),1:3);
AC=node(face(:,3),1:3)-node(face(:,1),1:3);

N=cross(AB',AC')';
d=-dot(N',node(face(:,1),1:3)')';

Rn2=1./sum((N.*N)')';

N1=cross(AC',N')'.*repmat(Rn2,1,3);
d1=-dot(N1',node(face(:,1),1:3)')';

N2=cross(N',AB')'.*repmat(Rn2,1,3);
d2=-dot(N2',node(face(:,1),1:3)')';

den=(v0*N')';
t=-(d+(p0*N')');
P=(p0'*den'+v0'*t')';
u=dot(P',N1')'+den.*d1;
v=dot(P',N2')'+den.*d2;

idx=find(den);
den(idx)=1./den(idx);

t=t.*den;
u=u.*den;
v=v.*den;

 % if den==0, ray is parallel to triangle, set t to infinity
t(find(den==0))=Inf;

if(nargout>=4)
    idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ip = inplace(a)
% INPLACE Convert to a type that supports inplace modification
% ipa = inplace(v) creates an inplace double object from any matrix v.
%
% Example:
%    ipa = inplace(ones(5));

if nargin == 0
    error('inplace:invalidParameter', ['ipdouble must be created with an initial .']);
elseif isa(a,'ipdouble')
    % make a copy
    ipd = inplace(a.a);

    return;

end

    function out = ip_get_a()
        out = a;
    end

    function ip_assign(in)
        a = in;
    end

    function ip_subsasgn(S,B)
        subsasgn(a,S,B);
    end

    ip.get_a = @ip_get_a;
    ip.assign = @ip_assign;
    ip.subsasgn = @ip_subsasgn;

    ip = class(ip,'inplace');

% ***** end inplace *****
end


text:
function [X, Y, B] = ginput_datetime(ax, varargin)
    % standard ginput, but can handle date_time. Axes must be specified
    % [x, y, button] = ginput_datetime(ax, ...)
    %
    % because, for large date ranges, there is an artificial level of precision,
    % the addition of a 'nearest_...' parameter (as the last parameter) will do rounding
    %
    % [x, y, button] = ginput_datetime(ax, ... , 'nearest_day')
    % [x, y, button] = ginput_datetime(ax, ... , 'nearest_hour')
    %
    %
    % see also ginput
    %
    % Celso G Reyes, 2017

    selector = '';
    if ~isempty(varargin) && ischar(varargin{end}) && startsWith(varargin{end},'nearest_')
        selector = varargin{end};
        varargin(end) = [];
    end

    axes(ax);
    [X, Y, B] = ginput(varargin{:});

    X = num2ruler(X, ax.XAxis);
    Y = num2ruler(Y, ax.YAxis);
    if isa(X,'datetime')
        X = round_time(X,selector);
    end
    if isa(Y,'datetime')
        Y = round_time(Y,selector);
    end

end

function dates = round_time(dates, selector)
    switch selector
        case 'nearest_day'
            dates = datetime(dates.Year, dates.Month, dates.Day + round(dates.Hour ./ 24));
        case 'nearest_hour'
            dates = datetime(dates.Year, dates.Month, dates.Day, dates.Hour + round(dates.Minute/60),0,0);
        otherwise
            % do nothing
    end
end


text:
%CODEGENERATOR.CONSTRUCTHEADERSTRINGC Creates common toolbox header string.
%
% HFSTRING = CGen.constructheaderstringc(HSTRUCT) is the formatted header
% string.
% HSTRUCT  is the header content struct
%
% Notes::
% The contents of the header are coded in a struct that is the input
% parameter and has the following self-explaining fields:
% Fieldname              Datatype
% - funName              'string'
% - shortDescription     'string'
% - detailedDescription  {'line1','line2',...}
% - inputs               {'input1: description','input2: description',...}
% - outputs              {'output1: description','output2: description',...}
% - examples             {'line1','line2',...}
% - knownBugs            {'line1','line2',...}
% - toDO                 {'line1','line2',...}
% - references           {'line1','line2',...}
% - authors              {'line1','line2',...}
% - seeAlso              {'function1','function2',...}
%
% Example::
% hStruct.funName = 'myFirstFunction';
% hStruct.shortDescription = ['Just an example!'];
% hStruct.calls = {'result = myFirstFunction(A,B)'};
% constructheaderstringc(hStruct)
%
% Authors::
%  Joern Malzahn   (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.constructheaderstring, sprintf.

% Copyright (C) 1993-2014, by Peter I. Corke
% Copyright (C) 2014, by Joern Malzahn
%
% 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
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

% function [hFString] = constructheaderstringc(CGen,hStruct)
%
% % reuse existing header construction routine,
% hFString = CGen.constructheaderstring(hStruct);
%
% % replace comment characters for C-compliance
% hFString = regexprep(hFString, '%', '//');

%CODEGENERATOR.CONSTRUCTHEADERSTRING Creates common toolbox header string.
%
% HFSTRING = CGen.constructheaderstring(HSTRUCT) is the formatted header
% string.
% HSTRUCT  is the geader content struct
%
% Notes::
% The contents of the header are coded in a struct that is the input
% parameter and has the following self-explaining fields:
% Fieldname              Datatype
% - funName              'string'
% - shortDescription     'string'
% - detailedDescription  {'line1','line2',...}
% - inputs               {'input1: description','input2: description',...}
% - outputs              {'output1: description','output2: description',...}
% - examples             {'line1','line2',...}
% - knownBugs            {'line1','line2',...}
% - toDO                 {'line1','line2',...}
% - references           {'line1','line2',...}
% - authors              {'line1','line2',...}
% - seeAlso              {'function1','function2',...}
%
% Example::
% hStruct.funName = 'myFirstFunction';
% hStruct.shortDescription = ['Just an example!'];
% hStruct.calls = {'result = myFirstFunction(A,B)'};
% constructheaderstring(hStruct)
%
% Authors::
%  Joern Malzahn   (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.replaceheader, sprintf.

% Copyright (C) 1993-2014, by Peter I. Corke
% Copyright (C) 2014, by Joern Malzahn
%
% 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
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [hFString] = constructheaderstringc(CGen,hStruct)

hFString = [];
hFString = [hFString,sprintf('%s\n',['/*! \file ',hStruct.funName,'.h'])];

if isfield(hStruct,'funName') && isfield(hStruct,'shortDescription')

    hFString = [hFString, sprintf('%s \n',['\brief ', hStruct.shortDescription, ''])];
    hFString = [hFString, sprintf('%s \n', '   ')]; % Empty line
else
    error('Header must include the function name and provide a short description!')
end

%% Detailed description
hFString = [hFString, sprintf('%s \n', ' ')];
if isfield(hStruct,'detailedDescription')
    nDescription = length(hStruct.detailedDescription);
    for iDescription = 1:nDescription
        hFString = [hFString, sprintf('%s \n', [hStruct.detailedDescription{iDescription}])];
    end
else
    hFString = [hFString, sprintf('%s \n', ' ')];
end
hFString = [hFString, sprintf('%s \n', '   ')]; % Empty line

%% Examples
hFString = [hFString, sprintf('%s \n', '__Example__:
')];
if isfield(hStruct,'example')
    nExamples = length(hStruct.example);
    for iExamples = 1:nExamples
        hFString = [hFString, sprintf('%s \n', [hStruct.example{iExamples}])];
    end
else
    hFString = [hFString, sprintf('%s \n', '    ')];
end

hFString = [hFString, sprintf('%s \n', '   ')];

%% Known Bugs
hFString = [hFString, sprintf('%s \n', '__Known Bugs__:
')];
if isfield(hStruct,'knownBugs')
    nBugs = length(hStruct.knownBugs);
    for iBugs = 1:nBugs
        hFString = [hFString, sprintf('%s \n', [hStruct.knownBugs{iBugs}])];
    end
else
    hFString = [hFString, sprintf('%s \n', ' ')];
end

hFString = [hFString, sprintf('%s \n', '   ')];

%% TODO list
hFString = [hFString, sprintf('%s \n', '__TODO__:
')];
if isfield(hStruct,'toDo')
    nToDo = length(hStruct.toDo);
    for iToDo = 1:nToDo
        hFString = [hFString, sprintf('%s \n', [hStruct.toDo{iToDo}])];
    end
else
    hFString = [hFString, sprintf('%s \n', ' ')];
end

hFString = [hFString, sprintf('%s \n', '   ')];

%% References
hFString = [hFString, sprintf('%s \n', '__References__:
')];
if isfield(hStruct,'references')
    nRef = length(hStruct.references);
    for iRef = 1:nRef
        hFString = [hFString, sprintf('+ %s \n', [hStruct.references{iRef}])];
    end
else
    hFString = [hFString, sprintf('%s \n', ' ')];
end

hFString = [hFString, sprintf('%s \n', '   ')];

%% Authors
hFString = [hFString, sprintf('%s \n', '__Authors__:
')];
if isfield(hStruct,'authors')
    nAuthors = length(hStruct.authors);
    for iAuthors = 1:nAuthors
        hFString = [hFString, sprintf('%s', [hStruct.authors{iAuthors}])];
    end
else
    hFString = [hFString, sprintf('%s \n', ' ')];
end

hFString = [hFString, sprintf('%s \n', '   ')];
hFString = [hFString, sprintf('%s \n', '

')];

%% Copyright note
hFString = [hFString, sprintf('%s \n', '__Copyright Note__:')];
cprNote = CGen.generatecopyrightnote;
cprNote = regexprep(cprNote, '%', ' ');
cprNote = regexprep(cprNote, 'Copyright', '
Copyright ');

hFString = [hFString,cprNote];

hFString = [hFString,sprintf('%s\n','*/')];
hFString = [hFString, sprintf('%s \n', '   ')]; % Empty line
hFString = [hFString,sprintf('%s\n',['/*!\fn ',hStruct.calls])];

%% Explanation of the Inputs
if isfield(hStruct,'inputs')
    nInputs = length(hStruct.inputs);
    for iInputs = 1:nInputs
        hFString = [hFString, sprintf('\\param %s\n', hStruct.inputs{iInputs})];
    end
else
    hFString = [hFString, sprintf('%s \n', '')];
end
%% Explanation of the Outputs
if isfield(hStruct,'outputs')
    nOutputs = length(hStruct.outputs);
    for iOutputs = 1:nOutputs
        hFString = [hFString, sprintf('\\param %s\n', hStruct.outputs{iOutputs})];
    end
else
    hFString = [hFString, sprintf('%s \n', ' ')];
end

hFString = [hFString,sprintf('%s\n','*/')];


text:
function yhat=spence21(y)
% function for Spencer's 21-point graduation rule.
% set out following Spencer's hand-calculation method,
% which isn't the shortest computer program!

oten = ones(10,1);
n = length(y);
y = [ oten*y(1); y; oten*y(n) ];

n = length(y);
k = 4:(n-3);
y1 = -y(k-3) + y(k-1) + 2*y(k) + y(k+1) -y(k+3);

n = length(y1);
k = 4:(n-3);
y2 = y1(k-3)+y1(k-2)+y1(k-1)+y1(k)+y1(k+1)+y1(k+2)+y1(k+3);

n = length(y2);
k = 3:(n-2);
y3 = y2(k-2)+y2(k-1)+y2(k)+y2(k+1)+y2(k+2);

n = length(y3);
k = 3:(n-2);
y4 = y3(k-2)+y3(k-1)+y3(k)+y3(k+1)+y3(k+2);

yhat = y4/350;
return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function T = ful(T)
%FUL Convert formatted data set to an array.
%   T = ful(T) converts the formatted data set T into a MATLAB array. If T
%   is incomplete, unknown entries are represented by the value NaN.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if isnumeric(T), return; end
if ~isstruct(T)
    error('ful:T','T must be a data set formatted by fmt.');
end

val = T.val;
size_tens = T.size;
if 8*prod(size_tens) > 8e9
    error('ful:size','T is too large to be stored as an array.');
end
if ~isfield(T,'ind')
    ind = sub2ind(size_tens,T.sub{:});
else
    ind = T.ind;
end

if isfield(T,'sparse') && T.sparse
    T = zeros(size_tens);
elseif isfield(T,'incomplete') && T.incomplete
    T = nan(size_tens);
else
    error('ful')
end
T(ind) = val;

author: andrewssobral

repo: lrslibrary