text:
function plot_varN2(params,nMode,const)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: plot_parvar1(params, 2,200,[],100)
%
% This function plots z, beta, and the probability of them as a function of
% input parameter like sampling size, window length and bin size. For each
% plot the other two parameters are kept fixed.
% The plot shows the distribution of values over all the grid point, the
% mean and the error as a standard deviation.
%
% Input
% params        structural array from
% nMode         parameter to plot z, beta, prob's against. Choose
%                      1 : Sampling Volume, N
%                      2 : Window Length, Tw
%                      3 : Bin Size, Tbin
% vN            Vector with range of sample sizes
% vTw           Vector with range of window length
% vTbin         Vector with range of bin size
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author:
% van Stiphout, Thomas, vanstiphout@sed.ethz.ch
%
% Created on 16.08.2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% prepare x-labels and ...
%         sXLabel='Sampling Volume, N';
switch nMode
    case 1
        %prepare calculation
        nPos1=1;
        vVariable2=params.vTw; nPos2=2;
        vVariable3=params.vTbin(const); nPos3=3;
        % perform calculation
        calc_plotval
        % prepare plotting
        sXLabel='Sampling Volume, N';
        sName=sprintf('z(N,Tw) - const. Tbin=%4.1f',vVariable3);
        sPrint=sprintf('RatesfN_Tw_const_Tbin-%03.0f.png',vVariable3);
        plot_values
    case 2
        %prepare calculation
        nPos1=1;
        vVariable2=params.vTbin;nPos2=3;
        vVariable3=params.vTw(const); nPos3=2;
        % perform calculation
        calc_plotval;
        % prepare plotting
        sXLabel='Sampling Volume, N';
        sName=sprintf('z(N,Tbin) - const. Tw=%4.1f',vVariable3);
        sPrint=sprintf('RatesfN_Tbin_const_Tw-%03.0f.png',vVariable3);
        plot_values;
    case 3
        %prepare calculation
        nPos1=2;
        vVariable2=params.vN;nPos2=1;
        vVariable3=params.vTbin(const); nPos3=3;
        % perform calculation
        calc_plotval;
        % prepare plotting
        sXLabel='Sampling Volume, Tw';
        sName=sprintf('z(Tw,N) - const. Tbin=%4.1f',vVariable3);
        sPrint=sprintf('RatesfTw_N_const_Tbin-%03.0f.png',vVariable3);
        plot_values;
    case 4
        %prepare calculation
        nPos1=2;
        vVariable2=params.vTbin;nPos2=3;
        vVariable3=params.vN(const); nPos3=1;
        % perform calculation
        calc_plotval;
        % prepare plotting
        sXLabel='Sampling Volume, Tw';
        sName=sprintf('z(Tw,Tbin) - const. N=%4.1f',vVariable3);
        sPrint=sprintf('RatesfTw_Tbin_const_N-%03.0f.png',vVariable3);
        plot_values;
    case 5
        %prepare calculation
        nPos1=3;
        vVariable2=params.vN; nPos2=1;
        vVariable3=params.vTw(const); nPos3=2;
        % perform calculation
        calc_plotval;
        % prepare plotting
        sXLabel='Sampling Volume, Tbin';
        sName=sprintf('z(Tbin,N) - const. Tw=%4.1f',vVariable3);
        sPrint=sprintf('RatesfTbin_N_const_Tw-%03.0f.png',vVariable3);
        plot_values;
    case 6
        %prepare calculation
        nPos1=3;
        vVariable2=params.vTw; nPos2=2;
        vVariable3=params.vN(const); nPos3=1;
        % perform calculation
        calc_plotval;
        % prepare plotting
        sXLabel='Sampling Volume, Tbin';
        sName=sprintf('z(Tbin,Tw) - const. N=%4.1f',vVariable3);
        sPrint=sprintf('RatesfTbin_Tw_const_N-%03.0f.png',vVariable3);
        plot_values;
end  % end switch

    function calc_plotval
        for ii=1:size(vVariable2,1)
            vSel=logical((params.mVar(:,nPos2) == vVariable2(ii)).*(params.mVar(:,nPos3)==vVariable3));
            params.mVar(vSel,:);
            mN.V1=params.mVar(vSel,nPos1);
            for i=1:size(params.mResult_,1)
                mN.V2.Z{ii}(:,i)=squeeze(params.mResult_(i,1,vSel));
                mN.V2.B{ii}(:,i)=squeeze(params.mResult_(i,2,vSel));
                mN.V2.pZ{ii}(:,i)=squeeze(params.mResult_(i,3,vSel));
                mN.V2.pB{ii}(:,i)=squeeze(params.mResult_(i,4,vSel));
            end  % end for i
            mN.V2.Z_mean(:,ii)=mean(mN.V2.Z{ii},2);
            mN.V2.B_mean(:,ii)=mean(mN.V2.B{ii},2);
            mN.V2.pZ_mean(:,ii)=mean(mN.V2.pZ{ii},2);
            mN.V2.pB_mean(:,ii)=mean(mN.V2.pB{ii},2);
            mN.V2.Z_std(:,ii)=std(mN.V2.Z{ii},0,2);
            mN.V2.B_std(:,ii)=std(mN.V2.B{ii},0,2);
            mN.V2.pZ_std(:,ii)=std(mN.V2.pZ{ii},0,2);
            mN.V2.pB_std(:,ii)=std(mN.V2.pB{ii},0,2);
        end % end for ii
    end  % end function calc_plotval

    function plot_values
        figure_w_normalized_uicontrolunits('Position',[100 25 400 600],'Name',sName);
        mPlot{1}.X1=mN.V1
        mPlot{1}.Y1=mN.V2.Z_mean;
        mPlot{1}.Y2=mN.V2.Z_mean-mN.V2.Z_std;
        mPlot{1}.Y3=mN.V2.Z_mean+mN.V2.Z_std;
        mPlot{2}.X1=mN.V1;
        mPlot{2}.Y1=mN.V2.B_mean;
        mPlot{2}.Y2=mN.V2.B_mean-mN.V2.B_std;
        mPlot{2}.Y3=mN.V2.B_mean+mN.V2.B_std;
        mPlot{3}.X1=mN.V1;
        mPlot{3}.Y1=mN.V2.pZ_mean;
        mPlot{3}.Y2=mN.V2.pZ_mean-mN.V2.pZ_std;
        mPlot{3}.Y3=mN.V2.pZ_mean+mN.V2.pZ_std;
        mPlot{4}.X1=mN.V1;
        mPlot{4}.Y1=mN.V2.pB_mean;
        mPlot{4}.Y2=mN.V2.pB_mean-mN.V2.pB_std;
        mPlot{4}.Y3=mN.V2.pB_mean+mN.V2.pB_std;

        sTitle=[cellstr(char('z(lta)-variation')), cellstr(char('p(Z)-variation')),...
            cellstr(char('\beta-variation')),  cellstr(char('p(\beta)-variation'))];
        sYLabel=[cellstr(char('z(lta)')),cellstr(char('p(Z)')),cellstr(char('\beta')), cellstr(char('p(\beta)'))];

        for pp=1:4
            subplot(2,2,pp);
            hold on;plot(mPlot{pp}.X1,mPlot{pp}.Y1,'x--','LineWidth',2,'MarkerSize',12)
            if pp==1  legend(num2str(vVariable2)); end
            hold on;plot(mPlot{pp}.X1,mPlot{pp}.Y2  ,'+:','LineWidth',1)
%             hold on;plot(mPlot{pp}.X1,mPlot{pp}.Y2,'+')
            hold on;plot(mPlot{pp}.X1,mPlot{pp}.Y3,':+','LineWidth',1)
%             hold on;plot(mPlot{pp}.X1,mPlot{pp}.Y3,'+')
            title(sTitle(pp));
            ylabel(sYLabel(pp));
            xlabel(sXLabel);
        end % end for pp
    end % end plot_values

    mN;
    sPrint=sprintf('print -dpng %s',sPrint);
    eval(sPrint)
end % end function

%
%     xlabel(sXLabel);
%     switch nSubplot
%         case 1
%             title('z(lta)-variation');
%             ylabel('z(lta)');
%         case 2
%             title('p(z)-variation');
%             ylabel('p(z)');
%         case 3
%             title('\beta-variation');
%             ylabel('\beta');
%         case 4
%             title('p(\beta)-variation');
%             ylabel('p(\beta)');
%     end
%     subplot(2,2,nSubplot);
%     Y=mean(squeeze(params.mResult_(:,nSubplot,vSel)));
%     E=std(squeeze(params.mResult_(:,nSubplot,vSel)));
%     errorbar(params.mVar(vSel,nMode),Y,E,'xr','MarkerSize',10,'LineWidth',1);
% end

author: CelsoReyes


text:
function a = lt( x, y )

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '<' );
if nargout, a = b; end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function data = readframeSUNRGBD(thispath,dataRoot,cls,bbmode)
    % example code to read annotation from ".json" file.
    % thispath: full path to the data folder.
    % dataRoot: root directory of all data folder.
    % cls : object category of ground truth to load. If not speficfy, the code will load all ground truth.
    if ~exist('cls','var')
        cls  =[];
    end
    if ~exist('bbmode','var')
        bbmode  ='2Dbb';
    end

    if ~exist('dataRoot','var')||isempty(dataRoot)
        dataRoot = '/n/fs/sun3d/data/';
    end
    sequenceName  = getSequenceName(thispath,dataRoot);
    if ~exist(thispath,'dir')
        data.sequenceName = sequenceName;
        data.valid = 0;
        return;
    end
    indd = find(sequenceName=='/');
    sensorType = sequenceName(indd(1)+1:indd(2)-1);
    % get K
    fID = fopen([thispath '/intrinsics.txt'],'r');
    K = reshape(fscanf(fID,'%f'),[3,3])';
    fclose(fID);

     % get image and depth path
    depthpath = dir([thispath '/depth/' '/*.png']);
    depthname = depthpath(1).name;
    depthpath = [thispath '/depth/' depthpath(1).name];

    rgbpath = dir([thispath '/image/' '/*.jpg']);
    rgbname = rgbpath(1).name;
    rgbpath = [thispath '/image/' rgbpath(1).name];

    if exist(sprintf('%s/annotation3Dfinal/index.json',thispath),'file')
        annoteImage =loadjson(sprintf('%s/annotation3Dfinal/index.json',thispath));
        % get Box
        filename = dir([fullfile(thispath,'extrinsics') '/*.txt']);
        Rtilt = dlmread([fullfile(thispath,'extrinsics') '/' filename(end).name]);
        Rtilt = Rtilt(1:3,1:3);
        anno_extrinsics = Rtilt;
        % convert it into matlab coordinate
        Rtilt = [1 0 0; 0 0 1 ;0 -1 0]*Rtilt*[1 0 0; 0 0 -1 ;0 1 0];

        cnt =1;
        for obji =1:length(annoteImage.objects)
                annoteobject =annoteImage.objects(obji);
                if ~isempty(annoteobject)&&~isempty(annoteobject{1})&&~isempty(annoteobject{1}.polygon)
                    annoteobject =annoteobject{1};
                    box = annoteobject.polygon{1};

                    % class name and label
                    ind = find(annoteobject.name==':');
                    if isempty(ind)
                        classname  = annoteobject.name;
                        labelname ='';
                    else
                        if ismember(annoteobject.name(ind-1),{'_',' '}),
                            clname = annoteobject.name(1:ind-2);
                        else
                            clname = annoteobject.name(1:ind-1);
                        end
                        %[~,classId]= ismember(clname,classNames);
                        classname  = clname;
                        labelname = annoteobject.name(ind+2:end);
                        %[~,label]= ismember(Labelname,labelNames);
                    end
                    if ismember(classname,{'wall','floor','ceiling'})||(~isempty(cls)&&~(sum(ismember(cls,{classname}))>0)),
                        continue;
                    end

                    x =box.X;
                    y =box.Z;
                    vector1 =[x(2)-x(1),y(2)-y(1),0];
                    coeff1 =norm(vector1);
                    vector1 =vector1/norm(vector1);
                    vector2 =[x(3)-x(2),y(3)-y(2),0];
                    coeff2 = norm(vector2);
                    vector2 =vector2/norm(vector2);
                    up = cross(vector1,vector2);
                    vector1 = vector1*up(3)/up(3);
                    vector2 = vector2*up(3)/up(3);
                    zmax =-box.Ymax;
                    zmin =-box.Ymin;
                    centroid2D = [0.5*(x(1)+x(3)); 0.5*(y(1)+y(3))];

                    thisbb.basis = [vector1;vector2; 0 0 1]; % one row is one basis
                    thisbb.coeffs = abs([coeff1, coeff2, zmax-zmin])/2;
                    thisbb.centroid = [centroid2D(1), centroid2D(2), 0.5*(zmin+zmax)];
                    thisbb.classname = classname;
                    thisbb.labelname = labelname;
                    thisbb.sequenceName = sequenceName;
                    orientation = [([0.5*(x(2)+x(1)),0.5*(y(2)+y(1))] - centroid2D(:)'), 0];
                    thisbb.orientation = orientation/norm(orientation);

                    if strcmp(bbmode,'2Dbb'),
                        [bb2d,bb2dDraw] = projectStructBbsTo2d(thisbb,Rtilt,[],K);
                        %gtBb2D = crop2DBB(gtBb2D,427,561);
                        thisbb.gtBb2D = bb2d(1:4);
                    end
                    groundtruth3DBB(cnt) =thisbb;
                    cnt=cnt+1;
                end
        end
        if cnt==1,groundtruth3DBB =[];end
    else
        groundtruth3DBB =[];
        filename = dir([fullfile(thispath,'extrinsics') '/*.txt']);
        Rtilt = dlmread([fullfile(thispath,'extrinsics') '/' filename(end).name]);
        Rtilt = Rtilt(1:3,1:3);
        anno_extrinsics = Rtilt;
        Rtilt = [1 0 0; 0 0 1 ;0 -1 0]*Rtilt*[1 0 0; 0 0 -1 ;0 1 0];

    end
    % read in room
    gtCorner3D =[];
    if exist([thispath '/annotation3Dlayout/index.json'],'file')
       json=loadjson([thispath '/annotation3Dlayout/index.json']);
       for objectID=1:length(json.objects)
           try
                groundTruth = json.objects{objectID}.polygon{1};
                numCorners = length(groundTruth.X);

                gtCorner3D(1,:) = [groundTruth.X groundTruth.X];
                gtCorner3D(2,:) = [repmat(groundTruth.Ymin,[1 numCorners]) repmat(groundTruth.Ymax,[1 numCorners])];
                gtCorner3D(3,:) = [groundTruth.Z groundTruth.Z];
                gtCorner3D = anno_extrinsics'*gtCorner3D;
                gtCorner3D = gtCorner3D([1,3,2],:);
                gtCorner3D(3,:) = -1*gtCorner3D(3,:);
                gtCorner3D = Rtilt*gtCorner3D;
                break;
           catch
           end
       end

    end

    data =struct('sequenceName',sequenceName,'groundtruth3DBB',...
         groundtruth3DBB,'Rtilt',Rtilt,'K',K,...
         'depthpath',depthpath,'rgbpath',rgbpath,'anno_extrinsics',anno_extrinsics,'depthname',depthname,...
         'rgbname',rgbname,'sensorType',sensorType,'valid',1,'gtCorner3D',gtCorner3D);

end

author: thusiyuan


text:
function varargout = msg(Flag, varargin)
% MSG Report process informations to the workspace.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% With this function informations about processes and variables can be printed
% to the workspace.
% ------------------------------------------------------------------------------
% INPUT
% 1 [Flag]
%   Flag describing the scope of the report. Any Flag can be used, however,
%   the following flags have a special purpose:
%   S ... to define the start of a process
%   E ... to define the end of a process
%   V ... to report the value of a variable
%   T ... to report a text string
%   D ... to report a table (of data)
%   O ... to set persistant options for this function (e.g. log level)
%   Conventions for arbitrary flags:
%   I ... for output of informations
%   X ... for error output
%   W ... for warnings
%
% 2 [ProcHierarchy] (not required for flags V, T, D and O)
%   Cell describing the hierarchy of the process.
%
% 3 [Info] (required for all flags, except S, E, D and O)
%   String containing report information.
%
% 4 ['Prec', Prec] (optional for flag V)
%   Precision of variable output.
%
% 5 ['MaxPrec', MaxPrec] (optional for flag V)
%   Maximal precision output for variables (see examples section).
%
% 6 ['FieldWidth', FieldWidth] (optional for flag V)
%   Field width for variable output (see examples section).
%
% 7 ['LogLevel', LogLevel] (optional for all flags, except O)
%   Define the log level of the message. Possible values:
%   * 'basic' -> for output of basic information.
%   * 'debug' -> for output of detailed information.
%
% 8 ['SetLogLevel', SetLogLevel] (optional for flag O)
%   Define which messages are reported to the workspace. Possible values:
%   * 'basic'  -> output of messages with log level 'basic'.
%   * 'debug'  -> output of messages with log level 'basic' and 'debug'.
%   * 'off'    -> no message output.
%
% 9 'GetLogLevel' (optional for flag O)
%   Get actual log level.
%
% 10 ['SetLogfile', SetLogfile] (optional for flag O)
%    Define a logfile instead of displaying the output information on the
%    workspace.
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 Report the start of a process and begin time measurement for this process.
%   msg('S', {'POINTCLOUD' 'ESTIMATE NORMALS'});
%
% 2 Report an info update for the same process.
%   msg('I', {'POINTCLOUD' 'ESTIMATE NORMALS'}, 'Normal estimation by PCA');
%
% 3 Report the end of the process.
%   msg('E', {'POINTCLOUD' 'ESTIMATE NORMALS'});
%
% 4 Print a variable.
%   msg('V', 6378137, 'Semi-major axis (GRS80)');
%
% 5 Print variable with specified precision.
%   msg('V', pi, 'pi', 'Prec', 7);
%
% 6 Print variable with definition of MaxPrec and FieldWidth.
%   msg('V', pi, 'pi', 'Prec', 7, 'MaxPrec', 10, 'FieldWidth', 20);
%
%   V :         3.1415927    = pi
%       12345678901234567890
%      >     FieldWidth     <
%                > MaxPrec  <
%                > Prec  <
%
% 8 Print a text string.
%   msg('T', 'A text information');
%
% 9 Print a data table.
%
%   tabularData.Caption = 'CAPTION OF DATA TABLE';
%   tabularData.ColumnCaptions     = {'value' 'usage'  'file'    % first row
%                                     ''      'yes/no' ''     }; % second row
%   tabularData.ColumnWidths       = [8        10       0     ];
%   tabularData.ColumnPrecs        = [2        NaN      NaN   ]; % NaN for char columns
%   tabularData.ColumnRightMargins = [1        3        0     ]; % optional
%
%   tabularData.Data = {7.654  'yes' 'file1.txt'
%                       6.3245 'no'  'file2_with_long_name.txt'
%                       8.12   'yes' 'file3.txt'};
%
%   tabularData.Notes = {'Note: files are located in folder ''D:\data\'''}; % optional
%
%   msg('D', tabularData);
%
% 10 Set log level to 'basic', so that only messages with log level equal to
%    'basic' are reported.
%    msg('O', 'SetLogLevel', 'basic');
%
% 11 Get actual log level.
%    origLogLevel = msg('O', 'GetLogLevel');
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Preparations -----------------------------------------------------------------

persistent logLevel % log level of output
if isempty(logLevel), logLevel = 'debug'; end % initialize

persistent logFile
if isempty(logFile), logFile = ''; end % initialize

% Set options
if strcmpi(Flag, 'O')

    % Set log level
    idx = find(strcmpi('SetLogLevel', varargin));
    if ~isempty(idx)
        logLevel = varargin{idx+1};
    end

    % Get log level
    idx = find(strcmpi('GetLogLevel', varargin));
    if ~isempty(idx)
        varargout{1} = logLevel;
    end

    % Set log file
    idx = find(strcmpi('SetLogFile', varargin));
    if ~isempty(idx)
        logFile = varargin{idx+1};
    end

    return;

% Retrieve log level of message
else

    % Check if log level of message is defined (without parsing all input arguments)
    idx = find(strcmpi('LogLevel', varargin));
    if ~isempty(idx) % if log level of message is defined
        logLevelMsg = varargin{idx+1};
        varargin(idx:idx+1) = []; % delete log level arguments (no parsing)
    else % if log level of message is NOT defined
        logLevelMsg = 'debug'; % default value
    end

end

% Check log level and eventually stop function ---------------------------------
% In order to save time, this section has to be placed at the beginning of the
% function

% If logLevel is set to off leave the function
if strcmpi(logLevel, 'off')
    return;

% If logLevel is set to basic ...
elseif strcmpi(logLevel, 'basic')

    % ... and log level of message is set to debug, leave the function
    if strcmpi(logLevelMsg, 'debug'), return; end

end

% ... continue with normal parsing of input arguments

% Input parsing ----------------------------------------------------------------

% Global parameters
p = inputParser;

if strcmpi(Flag, 'V')
    p.addRequired( 'Var'              , @(x) isnumeric(x) || islogical(x));
    p.addRequired( 'Info'             , @ischar);
    p.addParameter('Prec'         , 3 , @isnumeric);
    p.addParameter('MaxPrec'      , 9 , @isnumeric);
    p.addParameter('FieldWidth'   , 20, @isnumeric);

elseif any(strcmpi(Flag, {'S' 'E'}))
    p.addRequired(  'ProcHierarchy'    , @iscell);

elseif strcmpi(Flag, 'T')
    p.addRequired(  'Info'             , @ischar);

elseif strcmpi(Flag, 'D')
    p.addRequired(  'TabularData'      , @isstruct);

else
    p.addRequired(  'ProcHierarchy'    , @iscell);
    p.addRequired(  'Info'             , @ischar);
end

p.parse(varargin{:});
p = p.Results;

% Variable declaration ---------------------------------------------------------

persistent tictoc   % for time measurement

% Output string
out = [];

% Add time info to output ------------------------------------------------------

if any(strcmpi(Flag, {'V' 'T' 'D'}))
    % No time information, just empty string with same length
    tInfo = char( ones(1,numel(timeinfo)) *32 );

elseif strcmpi(Flag, 'S')
    tInfo = timeinfo;
    % Start timer for specific process
    tictoc.(genvarname([p.ProcHierarchy{:}])).tic = tic;

else
    % Elapsed time for specific process
    t = toc(tictoc.(genvarname([p.ProcHierarchy{:}])).tic);
    % No time information, just empty string with same length
    tInfo = char( ones(1,numel(timeinfo)) *32 );
    % Write elapsed time right-aligned into string
    tString = sprintf('+%.3fs', t);
    tInfo(end-numel(tString)+1:end) = tString;
end

out = [out tInfo];

% Add separator between two columns to output ----------------------------------

out = [out ' | '];

% Add flag to output -----------------------------------------------------------

out = [out Flag ' : '];

% Add process hierarchy string to output ---------------------------------------

if all(~strcmpi(Flag, {'V' 'T' 'D'}))
    for i = 1:numel(p.ProcHierarchy)
        switch i
            case 1
                procStr = upper(p.ProcHierarchy{i});
            otherwise
                procStr = [procStr ' > ' upper(p.ProcHierarchy{i})];
        end
    end
    out = [out procStr];
end

% Add info string to output ----------------------------------------------------

if ~any(strcmpi(Flag, {'S' 'E' 'V' 'T' 'D'}))
    out = [out ' : '];

    % If p.Info contains \n, replace them by left part of output
    logIdxSep = out == '|';
    leftPart = char( ones(1,numel(out)) *32 );
    leftPart(logIdxSep) = '|';
    p.Info = strrep(p.Info, sprintf('\n'), sprintf('\n%s', leftPart));

    out = [out p.Info];
end

% Add variable info to output --------------------------------------------------

if strcmpi(Flag, 'V')
    if p.Prec == 0, p.MaxPrec = p.MaxPrec+1; end
    out = [out sprintf('%*.*f%*s = %s', p.FieldWidth-p.MaxPrec+p.Prec, p.Prec, p.Var, p.MaxPrec-p.Prec, '', p.Info)];
end

% Add text info to output ------------------------------------------------------

if strcmpi(Flag, 'T')
    out = [out p.Info];
end

% Add tabular data to output ---------------------------------------------------

if strcmpi(Flag, 'D')

    % First: delete all columns containing only NaNs
    idx2del = [];
    for i = 1:size(p.TabularData.Data,2)
        if all(isnan([p.TabularData.Data{:,i}]))
            idx2del = [idx2del; i];
        end
    end
    p.TabularData.Data(:,idx2del) = [];
    p.TabularData.ColumnCaptions(:,idx2del) = [];
    p.TabularData.ColumnWidths(idx2del) = [];
    p.TabularData.ColumnPrecs(idx2del) = [];

    % Add column margins if not defined
    if ~isfield(p.TabularData, 'ColumnRightMargins')
        p.TabularData.ColumnRightMargins = [ones(1, size(p.TabularData.Data,2)-1) 0];
    end

    % Save prefix
    prefix = out;
    clear out

    % Separation line
    lengthSeparationLine = sum(p.TabularData.ColumnWidths) + sum(p.TabularData.ColumnRightMargins);
    if lengthSeparationLine < 80, lengthSeparationLine = 80; end
    out{1,1} = [prefix repmat('-', 1, lengthSeparationLine)];

    % Report table caption
    out{2,1} = [prefix p.TabularData.Caption];

    % Report column captions
    formatSpec = [];
    for col = 1:size(p.TabularData.ColumnCaptions,2)
        formatSpec = [formatSpec '%' num2str(p.TabularData.ColumnWidths(col)) 's' repmat(' ', 1,p.TabularData.ColumnRightMargins(col))];
    end
    formatSpec = strtrim(formatSpec); % remove trailing spaces
    for row = 1:size(p.TabularData.ColumnCaptions,1)
        out{end+1,1} = [prefix sprintf(formatSpec, p.TabularData.ColumnCaptions{row, :})];
    end

    % Report data
    formatSpec = [];
    for col = 1:size(p.TabularData.Data,2)
        if isnan(p.TabularData.ColumnPrecs(col))
            columnFormatSpec = ['%' num2str(p.TabularData.ColumnWidths(col)) 's'];
        else
            columnFormatSpec = ['%' num2str(p.TabularData.ColumnWidths(col)) '.' num2str(p.TabularData.ColumnPrecs(col)) 'f'];
        end
        formatSpec = [formatSpec columnFormatSpec repmat(' ', 1,p.TabularData.ColumnRightMargins(col))];
    end
    formatSpec = strtrim(formatSpec); % remove trailing spaces
    for row = 1:size(p.TabularData.Data,1)
        out{end+1,1} = [prefix sprintf(formatSpec, p.TabularData.Data{row, :})];
    end

    % Notes?
    if isfield(p.TabularData, 'Notes')
        for i = 1:numel(p.TabularData.Notes)
           out{end+1,1} = [prefix p.TabularData.Notes{i}];
        end
    end

    % Separation line
    out{end+1,1} = [prefix repmat('-', 1, lengthSeparationLine)];

end

% Write output to command window -----------------------------------------------

% Output to log file
if ~isempty(logFile)
    fid = fopen(logFile, 'at');
    if ~iscell(out)
        fprintf(fid, '%s\n', out);
    else % only in case of Flag = 'D'
        for i = 1:numel(out)
            fprintf(fid, '%s\n', out{i});
        end
    end
    fclose(fid);
% Output to work space
else
    if any(strcmpi(Flag, {'X'})) % standard error output
        fprintf(2, '%s\n', out);
    else % standard output
        if ~iscell(out)
            fprintf(1, '%s\n', out);
        else % only in case of Flag = 'D'
            for i = 1:numel(out)
                fprintf(1, '%s\n', out{i});
            end
        end
    end
end

% wk 2015-04-25 Wolfgang Rieger criticized that when redirecting the output of opalsICP to a file, the output stream seems to be buffered, and hence he cannot closely follow the progress of opalsICP when watching the file in a text editor.
% In fact, this seems to be the case for stdout and stderr, see: http://undocumentedmatlab.com/blog/improving-fwrite-performance
% Let's check if stdout and stderr can be flushed explicitly by the following command, see: http://stackoverflow.com/questions/2633019/how-can-i-flush-the-output-of-disp-in-matlab-or-octave
% wk 2015-04-26 this seems to be unnecessary, because the output of ICP4opals is directed to the opals logging system via the handler installed by opalsICP. Thus, flush the screen log stream (std::cout) in opalsICP instead.
% drawnow('update')

end % function end

% ------------------------------------------------------------------------------

function tInfo = timeinfo
tInfo = [datestr(now,'yyyy-mm-dd HH:MM:SS.FFF') 's'];
end

author: pglira


text:
function [X,E,obj,err,iter] = lrr(A,B,lambda,opts)

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

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

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

[d,na] = size(A);
[~,nb] = size(B);

X = zeros(nb,na);
E = zeros(d,na);
J = X;

Y1 = E;
Y2 = X;
BtB = B'*B;
BtA = B'*A;
I = eye(nb);
invBtBI = (BtB+I)\I;

iter = 0;
for iter = 1 : max_iter
    Xk = X;
    Ek = E;
    Jk = J;
    % first super block {J,E}
    [J,nuclearnormJ] = prox_nuclear(X+Y2/mu,1/mu);
    if strcmp(loss,'l1')
        E = prox_l1(A-B*X+Y1/mu,lambda/mu);
    elseif strcmp(loss,'l21')
        E = prox_l21(A-B*X+Y1/mu,lambda/mu);
    elseif strcmp(loss,'l2')
        E = mu*(A-B*X+Y1/mu)/(lambda+mu);
    else
        error('not supported loss function');
    end
    % second  super block {X}
    X = invBtBI*(B'*(Y1/mu-E)+BtA-Y2/mu+J);

    dY1 = A-B*X-E;
    dY2 = X-J;
    chgX = max(max(abs(Xk-X)));
    chgE = max(max(abs(Ek-E)));
    chgJ = max(max(abs(Jk-J)));
    chg = max([chgX chgE chgJ max(abs(dY1(:))) max(abs(dY2(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = nuclearnormJ+lambda*comp_loss(E,loss);
            err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

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

function out = comp_loss(E,loss)

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

author: canyilu

repo: LibADMM-toolbox


text:
function Y=domain(varargin)
%DOMAIN  Extract domain for PWA variable

% Fast version for plotting simple PWA objects
if nargin == 1
    X = varargin{1};
    if isa(varargin{1},'sdpvar')
        if length(X) == 1
            if length(getvariables(X))==1
                extstruct = yalmip('extstruct',getvariables(X));
                if ~isempty(extstruct)
                    if isequal(extstruct.fcn,'pwa_yalmip') | isequal(extstruct.fcn,'pwq_yalmip')
                        Y = [];
                        for i = 1:length(extstruct.arg{1})
                            Y = [Y extstruct.arg{1}{i}.Pn];
                        end
                        return
                    end
                end
            end
        end
    end
end
error('DOMAIN can only be applied to simple MPT related PWA objects')

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function ApToolupdatePHD_axesLabels(handles)
%update PHD axes labels

xlim = get(handles.phd,'xlim');
ylim = get(handles.phd,'ylim');

%setup axes labels
if isfield(handles,'phdXaxis')
    dxlim = diff(xlim);
    xticks = round(xlim(1))+[0 ceil(dxlim/4-1) ceil(dxlim/2-1) ...
        ceil(dxlim/4*3-1) dxlim];
    xticks(xticks < xlim(1)) = ceil(xlim(1));
    xticks(xticks > xlim(2)) = floor(xlim(2));
    xticklabels = cellfun(@(x) sprintf('%0.2g',x),num2cell(handles.phdXaxis(xticks)),'uniformoutput',false);
    xticklabels{3} = '0'; % Should be very close to zero, so this makes display cleaner.
    set(handles.phd,'xtick',xticks,'xticklabel',xticklabels);
    xlabel(handles.phd,'Polar Angle (degrees)','fontweight','bold')
else
    set(handles.phd, 'xtick', []);
end
if isfield(handles,'phdYaxis')
    dylim = diff(ylim);
    yticks = round(ylim(1))+[0 ceil(dylim/4-1) ceil(dylim/2-1) ...
        ceil(dylim/4*3-1) dylim];
    yticks(yticks < ylim(1)) = ceil(ylim(1));
    yticks(yticks > ylim(2)) = floor(ylim(2));
    set(handles.phd,'ytick',yticks,'yticklabel',...
        cellfun(@(x) sprintf('%0.4g',x),num2cell(handles.phdYaxis(yticks)),'uniformoutput',false));
    ylabel(handles.phd,sprintf('Frequency (%s)',handles.freqUnits),'fontweight','bold')
else
    set(handles.phd, 'ytick', []);
end

author: ngageoint


text:
%[2019]-"A new meta-heuristic optimizer: Pathfinder algorithm"

% (8/12/2020)

function PFA = jPathFinderAlgorithm(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;

if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = zeros(N,dim);
for i = 1:N
	for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
	end
end
% Fitness
fit  = zeros(1,N);
fitP = inf;
for i = 1:N
  fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  % Pathfinder update
  if fit(i) < fitP
    fitP = fit(i);
    Xpf  = X(i,:);
  end
end
% Set previous pathfiner
Xpf_old = Xpf;
% Pre
Xpf_new = zeros(1,dim);
Xnew    = zeros(N,dim);

curve = zeros(1,max_Iter);
curve(1) = fitP;
t = 2;
% Iterations
while t <= max_Iter
  % Alpha & beta in [1,2]
  alpha = 1 + rand();
  beta  = 1 + rand();
  for d = 1:dim
    % Define u2 in [-1,1]
    u2 = -1 + 2 * rand();
    % Compute A (2.6)
    A  = u2 * exp(-(2 * t) / max_Iter);
    % Update pathfinder (2.4)
    r3 = rand();
    Xpf_new(d) = Xpf(d) + 2 * r3 * (Xpf(d) - Xpf_old(d)) + A;
  end
  % Boundary
  Xpf_new(Xpf_new > ub) = ub; Xpf_new(Xpf_new <lb) = lb;
  % Update previous path
  Xpf_old = Xpf;
  % Fitness
  Fnew = fun(feat,label,(Xpf_new > thres),opts);
  % Greedy selection
  if Fnew < fitP
    fitP = Fnew;
    Xpf  = Xpf_new;
  end
  % Sort member
  [fit, idx] = sort(fit,'ascend');
  X          = X(idx,:);
  % Update first solution
  if Fnew < fit(1)
    fit(1) = Fnew;
    X(1,:) = Xpf_new;
  end
  % Update
  for i = 2:N
    % Distance (2.5)
    Dij = norm(X(i,:) - X(i-1,:));
    for d = 1:dim
      % Define u1 in [-1,1]
      u1  = -1 + 2 * rand();
      % Compute epsilon (2.5)
      eps = (1 - (t / max_Iter)) * u1 * Dij;
      % Define R1, R2
      r1  = rand();
      r2  = rand();
      R1  = alpha * r1;
      R2  = beta * r2;
      % Update member (2.3)
      Xnew(i,d) = X(i,d) + R1 * (X(i-1, d) - X(i,d)) + ...
        R2 * (Xpf(d) - X(i,d)) + eps;
    end
    % Boundary
    XB = Xnew(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    Xnew(i,:) = XB;
  end
  % Fitness
  for i = 2:N
    % Fitness
    Fnew = fun(feat,label,(Xnew(i,:) > thres),opts);
    % Selection
    if Fnew < fit(i)
      fit(i) = Fnew;
      X(i,:) = Xnew(i,:);
    end
    % Pathfinder update
    if fit(i) < fitP
      fitP = fit(i);
      Xpf  = X(i,:);
    end
  end
  curve(t) = fitP;
  fprintf('\nIteration %d Best (PFA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xpf > thres) == 1);
sFeat = feat(:,Sf);
% Store results
PFA.sf = Sf;
PFA.ff = sFeat;
PFA.nf = length(Sf);
PFA.c  = curve;
PFA.f  = feat;
PFA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
% caliper() - Measure a set of spatial components of a given data epoch relative to
%           a reference epoch and decomposition.
% Usage:
%   >> [amp,window]=caliper(newepoch,refepoch,weights,compnums,filtnums,times,'noplot');
%
% Inputs:
%     newepoch = (nchannels,ntimes) new data epoch
%     refepoch = a (nchannels,ntimes) reference data epoch
%     weights  = (nchannels,ncomponents) unmixing matrix (e.g., ICA weights*sphere)
%     compnums = vector of component numbers to return amplitudes for {def|0: all}
%     filtnums = [srate highpass lowpass] filter limits for refepoch {def|0: allpass}
%     times    = [start_ms end_ms] epoch latency limits, else latencies vector {def|0: 0:EEG.pnts-1}
%     'noplot' = produce no plots {default: plots windows for the first <=3 components}
%
% Outputs:
%     amps = (1,length(compnums)) vector of mean signed component rms amplitudes
%     windows = (length(compnums)),length(times)) matrix of tapering windows used
%
% Notes:
%   Function caliper() works as follows: First the reference epoch is decomposed using
%   the given weight matrix (may be ICA, PCA, or etc). Next, the time course of the
%   main lobe of the activation in the reference epoch (from max to 1st min, forward
%   and backward in time from abs max, optionally after bandpass filtering) is used
%   to window the new epoch. Then, the windowed new epoch is decomposed by the same
%   weight matrix, and signed rms amplitude (across the channels) is returned of the
%   projection of each of the specified component numbers integrated across the windowed
%   epoch. If not otherwise specified, plots the windows for the first <= 3 components listed.
%
% Example: Given a grand mean response epoch and weight matrix, it can be used
%   to measure amps of grand mean components in single-subject averages.
%
% Authors: Scott Makeig & Marissa Westerfield, SCCN/INC/UCSD, La Jolla, 11/2000

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

% Edit History:
% 12/05/00 -- added fig showing data, ref activation, and window vector -mw
% 12/19/00 -- adjusted new icaproj() args -sm
% 01-25-02 reformated help & license -ad

function [amps,windows] = caliper(newepoch,refepoch,weights,compnums,filtnums,times,noplot)

if nargin < 3
  help caliper
  return
end

if nargin<4
  compnums = 0;
end

nchans = size(newepoch,1);
ntimes = size(newepoch,2);

if compnums(1) == 0 | isempty(compnums(1))
   compnums = 1:nchans;
end

if min(compnums) < 1 | max(compnums) > size(weights,2)
   help caliper
   return
end

if nargin<5
  filtnums = [];
else
  if isempty(filtnums)
      filtnums = [];
  elseif length(filtnums)==1 & filtnums(1)==0
     filtnums = [];
  elseif length(filtnums) ~= 3
     fprintf('\ncaliper(): filter parameters (filtnums) must have length 3.\n')
     return
  end
end

if nargin< 6  | isempty(times) | (length(times)==1 & times(1)==0)
  times = 0:ntimes-1;
else
  if length(times) ~= ntimes
    if length(times) ~= 2
      fprintf('caliper(): times argument should be [startms endms] or vector.\n')
      return
    else
      times = times(1):(times(2)-times(1))/(ntimes-1):times(2);
      times = times(1:ntimes);
    end
  end
end

if nargin < 7
   noplot = 0;
else
   noplot = 1;
end

refact = weights(compnums,:)*refepoch; % size (length(compnums),ntimes)
newact = weights(compnums,:)*newepoch;

if length(filtnums) == 3
  if ~exist('eegfilt')
    fprintf('caliper(): function eegfilt() not found - not filtering refepoch.\n');
  else
   try
      refact = eegfilt(refact,filtnums(1),filtnums(2),filtnums(3));
   catch
      fprintf('\n%s',lasterr)
      return
   end
  end
end

if size(weights,1) == size(weights,2)
   winv = inv(weights);
else
   winv = pinv(weights);
end
winvrms = sqrt(mean(winv.*winv)); % map rms amplitudes

amps = [];
windows = [];
n = 1; % component index
for c=compnums
  if floor(c) ~= c
    help caliper
    fprintf('\ncaliper(): component numbers must be integers.\n')
    return
  end

  [lobemax i] = max(abs(refact(n,:)));
  f = i+1;
  oldact = lobemax;
  while f>0 & f<=ntimes & abs(refact(n,f)) < oldact
    oldact = abs(refact(n,f));
    f = f+1;
  end % f is now one past end of "main lobe"
  lobeend = f-1;

  f = i-1;
  oldact = lobemax;
  while f>0 & f<=ntimes & abs(refact(n,f)) < oldact
    oldact = abs(refact(n,f));
    f = f-1;
  end % f is now one past start of "main lobe"
  lobestart = f+1;

  refact(n,1:lobestart) = zeros(1,length(1:lobestart));
  refact(n,lobeend:end) = zeros(1,length(lobeend:ntimes));
  windows = [windows; refact(n,:)];

  refnorm = sum(refact(n,:));
  if abs(refnorm)<1e-25
    fprintf('caliper(): near-zero activation for component %d - returning NaN amp.\n',c)
    amps = [amps NaN];
  else
    refact(n,:) = refact(n,:)/refnorm; % make reference epoch window sum to 1
    amps = [amps winvrms(c)*sum(refact(n,:).*newact(n,:))];
  end
  n = n+1;
  if ~noplot & n <= 4  %%% only plot out at most the first 3 components
     refproj = icaproj(refepoch,weights,c);
     refproj = env(refproj);
     windproj = winv(:,c)*(refact(n-1,:)*refnorm);
     windproj = env(windproj);
     figure; plot(times,newepoch(2:nchans,:),'g');
     hold on;h=plot(times,newepoch(1,:),'g');
     hold on;h=plot(times,newepoch(1,:),'g',...
                    times,refproj(1,:),'b',...
                    times,windproj(1,:),'r','LineWidth',2);
     hold on;plot(times,refproj(2,:),'b',times,windproj(2,:),'r','LineWidth',2);
     set(h(1),'linewidth',1);
     legend(h,'new data','comp. act.','window');
     title(['Component ',int2str(c),';  rms amplitude = ',num2str(amps(n-1))],...
            'FontSize',14);
     ylabel('Potential (uV)');
  end;             %% endif
end % c

author: PatternRecognition

repo: OpenBMI


text:
classdef RollAngleConstraint < AbstractConstraint
    %RollAngleConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = RollAngleConstraint(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            ut = stateLogEntry.time;
            rVect = stateLogEntry.position;
            vVect = stateLogEntry.velocity;
            bodyInfo = stateLogEntry.centralBody;

            [rollAngle, ~, ~] = stateLogEntry.attitude.getEulerAngles(ut, rVect, vVect, bodyInfo);
            value = rad2deg(rollAngle);

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent).deepCopy();

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent).deepCopy();

                    otherwise
                        error('Unknown event node.');
                end

                cartElem = stateLogEntryStateComp.getCartesianElementSetRepresentation();
                cartElem = cartElem.convertToFrame(stateLogEntry.centralBody.getBodyCenteredInertialFrame());
                stateLogEntryStateComp.setCartesianElementSet(cartElem);

                ut = stateLogEntryStateComp.time;
                rVect = stateLogEntryStateComp.position;
                vVect = stateLogEntryStateComp.velocity;
                bodyInfo = stateLogEntryStateComp.centralBody;

                [rollAngle, ~, ~] = stateLogEntryStateComp.attitude.getEulerAngles(ut, rVect, vVect, bodyInfo);
                valueStateComp = rad2deg(rollAngle);
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Roll Angle';
        end

%         function name = getName(obj)
%             name = sprintf('%s - Event %i', obj.getConstraintType(), obj.event.getEventNum());
%         end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'deg';
            lbLim = -360;
            ubLim = 360;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
            addConstraintTf = output.output{1};
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = RollAngleConstraint(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%DEMOPOLYGONCENTROIDS  One-line description here, please.
%
%   output = demoPolygonCentroids(input)
%
%   Example
%   demoPolygonCentroids
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-02-16,    using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2022 INRAE.

%% Demo polygon

% vertices of a "U"-shaped polygon
poly = [0 0;30 0;30 20;20 20;20 11;10 11;10 20;0 20];

% display the filled polygon and the open polyline
figure; hold on; axis equal; axis([-5 35 -5 30]);
fillPolygon(poly, 'c');
hPoly = drawPolygon(poly, 'linewidth', 2, 'color', 'k');
drawVertices(poly);

%% Computes various centroids

% computes the centroids
vxCentroid = centroid(poly);
lsCentroid = polylineCentroid(poly, 'closed');
pgCentroid = polygonCentroid(poly);

% display the centroids
hVxc = drawPoint(vxCentroid, 'marker', 'o', 'color', [0.8 0 0], 'linewidth', 2, 'MarkerFaceColor', 'w');
hLsc = drawPoint(lsCentroid, 'marker', 'o', 'color', [0 0.8 0], 'linewidth', 2, 'MarkerFaceColor', 'w');
hPgc = drawPoint(pgCentroid, 'marker', 'o', 'color', [0 0 0.8], 'linewidth', 2, 'MarkerFaceColor', 'w');

% decorate graph
legend([hPoly hVxc hLsc hPgc], ...
    {'Polygon', 'Vertex Centroid', 'Polyline Centroid', 'Polygon Centroid'}, ...
    'Location', 'NorthEast');

print(gcf, 'polygon_centroids.png', '-dpng');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Y = HealpixNestedInvLinearTrans(X, A, m)
%
% Description
% This function calculates Y = WX where
% W = kron( eye(12), W_0 * W_1 * W_2 * ... * W_{m-2} * W_{m-1} )
% W_i = kron( eye( 4^{m-i-1} ), V_i )
% V_i = kron( eye(4), diag([0 1 1 ... 1]) ) + kron( A, diag([1 0 0 ...0]) )
%
% Parameters
% X : input row-vector
% Y : output row-vector (same size as X)
% A : linear transform applied to the adjacent elements
% m : log4(length(X) / 12)

function Y = HealpixNestedInvLinearTrans(X, A, m)

for itr = (m - 1):(-1):0
    num_trans = 12 * 4^(m - itr - 1);
    interval = 4^itr;
    begin_idx = 1;
    %itr
    for trans = 0:(num_trans - 1)
        IDX = begin_idx + interval * [0:3];
        %IDX
        X(IDX) = A * X(IDX);
        begin_idx = begin_idx + 4 * interval;
    end
end

Y = X;

author: mtex-toolbox

repo: mtex

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


text:
function Y = symNormalRnd(M, S, sz)

% Generates independent symmetric matrices
% according to a matrix-variate symmetric normal distribution.
%
% Y = symNormalRnd(M, S, sz)
%
%
% Input:
%   M       pxp mean symmetric matrix
%   S       qxq covariance matrix
%   sz      size of the array (e.g. number of samples)
%
% Output:
%   Y       pxpx[sz] array of random symmetric matrices.
%
% Copyright: Armin Schwartzman, 2009
%

% HISTORY:
%   2008.12.30 ASH (armins@hsph.harvard.edu) wrote it.

% Check inputs
if (size(M,1) ~= size(M,2)),
    error('Wrong input format');
end
if (size(S,1) ~= size(S,2)),
    error('Wrong input format');
end
p = size(M,1);
q = p*(p+1)/2;
if (size(S,1) ~= q),
    error('Wrong input format');
end
if ~exist('sz'),
    sz = 1;
end

%-----------------------------------------------------------------
% Generate multivariate normal
N = prod(sz);       % number of samples
z = mvnrnd(zeros(1,q), S, N);   % size Nxq
z = reshape(z', [q sz]);        % size qxN
Y = repmat(M, [1 1 sz]) + vecdinv(z);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [A_full, E_full, t_l1f] = l1_filter(D_full, column_seed, row_seed, A_seed)

t_l1f = 0;

[m, n] = size(D_full);
if nargin < 4
    D_seed = D_full(row_seed, column_seed);
    %[A_seed, ~] = ialm_rpca(D_seed);
    [A_seed, ~] = inexact_alm_rpca(D_seed);
end

rA = rank(A_seed);

[U,S,V] = svd(A_seed);
Ur = U(:,1:rA);
Vr = V(:,1:rA);
Sr = diag(S);
Sr = diag(Sr(1:rA));

%ts = tic;

As_row = Ur*Sr;
As_column = Sr*Vr';

A_row = zeros(rA, n);
A_column = zeros(rA, m);

column_comp = setdiff(1:n, column_seed);
pinvUr = Ur';
[~, Ac_row] = solve_ml1(D_full(row_seed, column_comp), Ur, pinvUr);
A_row(:, column_comp) = Ac_row;
A_row(:, column_seed) = As_column;

row_comp = setdiff(1:m, row_seed);
pinvVr = Vr';
[~, Ac_column] = solve_ml1(D_full(row_comp, column_seed)', Vr, pinvVr);
A_column(:, row_comp) = Ac_column;
A_column(:, row_seed) = As_row';

A_full = A_column'*inv(Sr)*A_row;
E_full = D_full - A_full;

%t_l1f = toc(ts);

author: andrewssobral

repo: lrslibrary


text:
function ci = cspselect_fixedNumber(score, ~, nComps, mode)
%CSPSELECT_FIXEDNUMBER - Select a fixed number of components
%
%Synopsis:
% CI = select_fixedNumber(SCORE,~,NCOMPS,MODE)
%
%Arguments:
% SCORE  - score of components
% NCOMPS - number of components (total number depends on MODE)
% MODE   - 'absolutemax' : chooses NCOMPS components corresponding to the
%          NCOMPS maximal absolute scores)
%          'equalperclass' : chooses 2*NCOMPS components corresponding to
%          the NCOMPS lowest and the NCOMPS highest scores (default)
%          'onlyclass1' : chooses NCOMPS components corresponding to the
%          NCOMPS lowest scores
%          'onlyclass2' : chooses NCOMPS components corresponding to the
%          NCOMPS highest scores
%
%Returns:
% CI     - index of components
%
%See also processing/proc_csp

nChans = length(score);

switch mode
    case 'absolutemax'
        [~,ix] = sort(-abs(score));
        ci = ix(1:nComps);
    case 'equalperclass'
        [~,ix] = sort(score);
        ci = [ix(1:nComps); ix(end:-1:nChans-nComps+1)];
    case 'onlyclass1'
        [~,ix] = sort(score);
        ci = ix(1:nComps);
    case 'onlyclass2'
        [~,ix] = sort(score);
        ci = ix(end:-1:nChans-nComps+1);
    otherwise
        error('Unknown mode')
end

author: bbci


text:
function failed_old_halfspace_bug243

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY

% this tests for bug 243

% create vol conductor for halfspace (in mm)
geom.pnt = [[0 2 0];[0 0 0];[1 2 0]];
vol_hs = ft_headmodel_halfspace(geom,[-1 -1 -1],'conductivity',1);

% FIXME: this solution takes into account only dipole projected on surface
%   choose adequate electrodes set

  % create vol conductor for sphere (radius = 100 m)
  r = 1e5; % in mm
  pnt = r*mesh_sphere(162);
  pnt(:,3) = pnt(:,3) - r;
  % figure,ft_plot_mesh(pnt),axis on
  vol_sph = ft_headmodel_singlesphere(pnt,'conductivity',1);

% FIXME: test infinite vs halfspace when both electrodes and dipoles are
% very far away from plane

% FIXME: test halfpsace versus BEM when the surface is very big

% create infinite volume conductor
vol = [];
vol.type = 'infinite';

% create electrodes positions
elec.pnt = [[0 0 -0.2];[0 0 -0.4];[0 0 -0.6];[0 0 -0.8];[0 0 -1];[0 0 -1.2];[0 0 -1.4];[0 0 -1.6]];
for i=1:8
  elec.label{i} = sprintf('ch%0.2d',i);
end

% apply the necessary sensors/volume transformations
[vol_hs1, sens1]  = ft_prepare_vol_sens(vol_hs, elec);
% [vol_sph1, sens2] = ft_prepare_vol_sens(vol_sph, elec);

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

% 1 first try, calculate lf for certain dipoles positions
pos = [[0 -1 -1];[0 -1 -2];[-1 -2.5 -1]];Ndip = size(pos,1);
ori = [[0 0 1];[0 0 1];[0 0 1]]; ori = ori'; mori = ori(:);
[lf]  = ft_compute_leadfield(pos, sens1, vol);
lf  = lf *mori;
[lf1] = ft_compute_leadfield(pos, sens1, vol_hs1);
lf1 = lf1 *mori;

figure,plot(lf)
figure,plot(lf1)

% 2 near to surface
pos = [[0 -1 -0.1];[0 -1 -0.2];[-1 -2.5 -0.5]];Ndip = size(pos,1);
[lf]  = ft_compute_leadfield(pos, sens1, vol);
[lf1] = ft_compute_leadfield(pos, sens1, vol_hs1);
figure,plot(lf)
figure,plot(lf1)

% 3 far away from surface
pos = [[0 -1 -10];[0 -1 -2];[-1 -2.5 -5]];Ndip = size(pos,1);
[lf]  = ft_compute_leadfield(pos, sens1, vol); lf = lf *mori;
[lf1] = ft_compute_leadfield(pos, sens1, vol_hs1);lf1 = lf1 *mori;
figure,plot(lf)
figure,plot(lf1)

% 4 far away and electrode-dipoles are near
elec.pnt = 1e6*[[0 0.1 -0.2];[0 0 -0.3];[0 0.2 -0.2];[0 0.3 -0.3];[0 0.7 -.2];[0 1 -.3];[0 1.1 -.2];[0 1.4 -.3]];
pos = 1e7*[[0 -1 -1];[0 -1 -2];[0 -2.5 -2]];Ndip = size(pos,1);
[lf]  = ft_compute_leadfield(pos, elec, vol); %lf = lf *mori;
[lf1] = ft_compute_leadfield(pos, elec, vol_hs1); %lf1 = lf1 *mori;
figure,plot(lf)
figure,plot(lf1)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Copyright (C) Daphne Koller, Stanford University, 2012

function [MEU OptimalDecisionRule] = OptimizeWithJointUtility( I )
  % Inputs: An influence diagram I with a single decision node and one or more utility nodes.
  %         I.RandomFactors = list of factors for each random variable.  These are CPDs, with
  %              the child variable = D.var(1)
  %         I.DecisionFactors = factor for the decision node.
  %         I.UtilityFactors = list of factors representing conditional utilities.
  % Return value: the maximum expected utility of I and an optimal decision rule
  % (represented again as a factor) that yields that expected utility.
  % You may assume that there is a unique optimal decision.

  % This is similar to OptimizeMEU except that we must find a way to
  % combine the multiple utility factors.  Note: This can be done with very
  % little code.

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % YOUR CODE HERE
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Uti = I.UtilityFactors(1);
 for i = 2:length(I.UtilityFactors)
	 Uti = FactorSum(Uti,I.UtilityFactors(i));
 end
 I.UtilityFactors = Uti;
[MEU OptimalDecisionRule] = OptimizeMEU(I);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [f,P,prob] = lomb(t,h,ofac,hifac)
% LOMB(T,H,OFAC,HIFAC) computes the Lomb normalized periodogram (spectral
% power as a function of frequency) of a sequence of N data points H,
% sampled at times T, which are not necessarily evenly spaced. T and H must
% be vectors of equal size. The routine will calculate the spectral power
% for an increasing sequence of frequencies (in reciprocal units of the
% time array T) up to HIFAC times the average Nyquist frequency, with an
% oversampling factor of OFAC (typically >= 4).
%
% The returned values are arrays of frequencies considered (f), the
% associated spectral power (P) and estimated significance of the power
% values (prob).  Note: the significance returned is the false alarm
% probability of the null hypothesis, i.e. that the data is composed of
% independent gaussian random variables.  Low probability values indicate a
% high degree of significance in the associated periodic signal.
%
% Although this implementation is based on that described in Press,
% Teukolsky, et al. Numerical Recipes  In C, section 13.8, rather than using
% trigonometric rercurrences, this takes advantage of MATALB's array
% operators to calculate the exact spectral power as defined in equation
% 13.8.4 on page 577.  This may cause memory issues for large data sets and
% frequency ranges.
%
% Example
%    [f,P,prob] = lomb(t,h,4,1);
%    plot(f,P)
%    [Pmax,jmax] = max(P)
%    disp(['Most significant period is ',num2str(1/f(jmax)),...
%         ' with FAP of ',num2str(prob(jmax))])
%
% Written by Dmitry Savransky 21 May 2008

%sample length and time span
N = length(h);
T = max(t) - min(t);

%mean and variance
mu = mean(h);
s2 = var(h);

%calculate sampling frequencies
f = (1/(T*ofac):1/(T*ofac):hifac*N/(2*T)).';

%angular frequencies and constant offsets
w = 2*pi*f;
tau = atan2(sum(sin(2*w*t.'),2),sum(cos(2*w*t.'),2))./(2*w);

%spectral power
cterm = cos(w*t.' - repmat(w.*tau,1,length(t)));
sterm = sin(w*t.' - repmat(w.*tau,1,length(t)));
P = (sum(cterm*diag(h-mu),2).^2./sum(cterm.^2,2) + ...
     sum(sterm*diag(h-mu),2).^2./sum(sterm.^2,2))/(2*s2);

%estimate of the number of independent frequencies
M=2*length(f)/ofac;

%statistical significane of power
prob = M*exp(-P);
inds = prob > 0.01;
prob(inds) = 1-(1-exp(-P(inds))).^M;

author: Sable

repo: mcbench-benchmarks


text:
function A = genWattsStrogotzGraph(N,K,rewireProb)
%%GENWATTSSTROGOTZGRAPH Generates a random instance of the Watts-Strogotz
%                       model as described in [1].
%
%INPUT:
% N: The number of nodes to be used in generating the graph.
% K: The number of connected neighbors for each node. This must be an even
%    number.
% rewireProb: The probability that any single edge connecting a given node
%             to its neighbor will be changed to connect to a different
%             node chosen uniformly at random from the set of nodes which
%             do not at that time share an edge with the given node.
%
%OUTPUT:
% A: The adjacency matrix for the final graph.
%
%EXAMPLE: Generates an instance of the Watts-Strogotz model.
% N = 30;
% K = 4;
% rewireProb = 0.3;
% A = genWattsStrogotzGraph(N,K,rewireProb);
% g = graph(A);
% plot(g,"MarkerSize",degree(g),"Layout","circle")
%
%REFERENCES:
%[1] D. J. Watts and S. H. Strogatz, "Collective dynamics of 'small-world'
%    networks," Nature, vol. 393, no. 6684, pp. 440-442, 1998.
%
%August 2022 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if mod(K,2)==1
    error('K must be even')
end
if rewireProb>1 || rewireProb<0
    error('rewireProb must be in the interval [0,1]')
end

A = zeros(N);
for i = 1:N-1
    for j = i+1:N
        A(i,j) = mod(abs(i-1-(j-1)),N-K/2)<=K/2;
    end
end
for i = 1:N
    for j = i+1:i+K/2
        jNode = mod(j-1,N)+1;
        if rand()<rewireProb
            A(i,jNode) = 0;
            free = find(A(i,:)==0);
            newNode = i;
            while newNode==i
                newNode = free(randi(length(free)));
            end
            if newNode<i
                A(newNode,i) = 1;
            else
                A(i,newNode) = 1;
            end
        end
    end
end
A = A+A';
end
%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function prob_test039 ( )

%*****************************************************************************80
%
%% TEST039 tests COSINE_MEAN, COSINE_SAMPLE, COSINE_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST039\n' );
  fprintf ( 1, '  For the Cosine PDF:\n' );
  fprintf ( 1, '  COSINE_MEAN computes the mean;\n' );
  fprintf ( 1, '  COSINE_SAMPLE samples;\n' );
  fprintf ( 1, '  COSINE_VARIANCE computes the variance.\n' );

  a = 2.0;
  b = 1.0;

  check = cosine_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST039 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = cosine_mean ( a, b );
  variance = cosine_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );
  fprintf ( 1, '  PDF mean =        %14f\n', mean );
  fprintf ( 1, '  PDF variance =    %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = cosine_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%NLABELD Return numeric labels of classified dataset (c
%
% 	NLABELS = NLABELD(Z)
% 	NLABELS = Z*NLABELD
%	  NLABELS = NLABELD(A,W)
% 	NLABELS = A*W*NLABELD
%
% INPUT
%		Z        Classified dataset, or
%		A,W      Dataset and classifier mapping
%
% OUTPUT
%		NLABELS	 Column vector of numeric labels)
%
% DESCRIPTION
% Returns the numberic labels of the classified dataset Z (typically the
% result of a mapping or classification A*W). For each object in Z (i.e.
% each row) the feature label or class label (i.e. the column label) of the
% maximum column value is returned. This corresponds with the classes
% stored in W, which can be found by GETLABELS(W).
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, TESTC, PLOTC, GETLABELS

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function labels = nlabeld(a,w)

		if (nargin == 0)

		% Untrained mapping.
		labels = prmapping(mfilename,'fixed');

	elseif (nargin == 1)

		% In a classified dataset, the feature labels contain the output
		% of the classifier.
		[m,k] = size(a); featlist = getfeatlab(a);

		if (k == 1)
			% If there is one output, assume it's a 2-class discriminant:
			% decision boundary = 0.
			J = 2 - (double(a) >= 0);
		else
			% Otherwise, pick the column containing the maximum output.
			[dummy,J] = max(+a,[],2);
		end
		labels = J;
	elseif (nargin == 2)

		% Just construct classified dataset and call again.
		labels = feval(mfilename,a*w);

	else
		error ('too many arguments');
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function x = merl_quadprog(varargin)
  % MERL_QUADPROG Optimize problems of the form:
  %
  % min_x 1/2 x' Q x - x' h
  %
  % s.t. x>=0
  %
  % According to the method described in "PARALLEL QUADRATIC PROGRAMMING FOR
  % IMAGE PROCESSING" [Brand & Chen 11]
  %
  % x = merl_quadprog(Q,h,'ParamterName',ParameterValue)
  %
  % Inputs:
  %   Q  n by n symmetric positive semi-definite quadratic coefficients matrix
  %   h  n by 1 linear coefficents vector (note sign)
  % Outputs:
  %   x  n by 1 solution vector
  %
  %

  Q = varargin{1};
  h = varargin{2};

  n = size(Q,1);
  assert(size(Q,2) == n);
  assert(numel(h) == n);
  h = h(:);

  max_iter = inf;
  tol = 1e-7;

  % initial guess
  x0 = ones(n,1);

  v = 3;
  while v <= numel(varargin)
    switch varargin{v}
    case 'MaxIter'
      assert((v+1)<=numel(varargin));
      v = v+1;
      max_iter = varargin{v};
    case 'X0'
      assert((v+1)<=numel(varargin));
      v = v+1;
      x0 = varargin{v};
    otherwise
      error(['Unknown parameter: ' varargin{v}]);
    end
    v = v+1;
  end

  %% "some r ∈ R^n ≥ 0"
  %r = sparse(n,1);
  % ri = max(Qii , ∑ j Q−i j )
  r = max(diag(Q),max(-Q,[],2));

  % "similarly" ... "si>0"
  s = ones(n,1);

  Qp = max(Q,0) + diag(r);
  Qm = max(-Q,0) + diag(r);
  assert(issparse(Qp));
  hp = max(h,0) + s;
  hm = max(-h,0) + s;

  x = x0;
  it = 0;
  while true
    x_prev = x;
    x = x.*(hp+Qm*x)./(hm+Qp*x);
    it = it + 1;
    diff = max(abs(x-x_prev));
    %fprintf('%d: %g\n',it,diff);
    if it >= max_iter
      return;
    end
    if diff<=tol
      return;
    end
  end

end

author: alecjacobson

repo: gptoolbox


text:
function net = cnnsetup(net, x, y)
    assert(~isOctave() || compare_versions(OCTAVE_VERSION, '3.8.0', '>='), ['Octave 3.8.0 or greater is required for CNNs as there is a bug in convolution in previous versions. See http://savannah.gnu.org/bugs/?39314. Your version is ' myOctaveVersion]);
    inputmaps = 1;
    mapsize = size(squeeze(x(:, :, 1)));

    for l = 1 : numel(net.layers)   %  layer
        if strcmp(net.layers{l}.type, 's')
            mapsize = mapsize / net.layers{l}.scale;
            assert(all(floor(mapsize)==mapsize), ['Layer ' num2str(l) ' size must be integer. Actual: ' num2str(mapsize)]);
            for j = 1 : inputmaps
                net.layers{l}.b{j} = 0;
            end
        end
        if strcmp(net.layers{l}.type, 'c')
            mapsize = mapsize - net.layers{l}.kernelsize + 1;
            fan_out = net.layers{l}.outputmaps * net.layers{l}.kernelsize ^ 2;
            for j = 1 : net.layers{l}.outputmaps  %  output map
                fan_in = inputmaps * net.layers{l}.kernelsize ^ 2;
                for i = 1 : inputmaps  %  input map
                    net.layers{l}.k{i}{j} = (rand(net.layers{l}.kernelsize) - 0.5) * 2 * sqrt(6 / (fan_in + fan_out));
                end
                net.layers{l}.b{j} = 0;
            end
            inputmaps = net.layers{l}.outputmaps;
        end
    end
    % 'onum' is the number of labels, that's why it is calculated using size(y, 1). If you have 20 labels so the output of the network will be 20 neurons.
    % 'fvnum' is the number of output neurons at the last layer, the layer just before the output layer.
    % 'ffb' is the biases of the output neurons.
    % 'ffW' is the weights between the last layer and the output neurons. Note that the last layer is fully connected to the output layer, that's why the size of the weights is (onum * fvnum)
    fvnum = prod(mapsize) * inputmaps;
    onum = size(y, 1);

    net.ffb = zeros(onum, 1);
    net.ffW = (rand(onum, fvnum) - 0.5) * 2 * sqrt(6 / (onum + fvnum));
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function ds_info=my_gen_ds_info_egohands(ds_config)

ds_dir=fullfile('../datasets', 'EgoHands');

train_idx_file=fullfile(ds_dir, 'ImageSets/Segmentation/train.txt');
fid=fopen(train_idx_file);
train_file_names=textscan(fid, '%s');
train_file_names=train_file_names{1};
fclose(fid);

val_idx_file=fullfile(ds_dir, 'ImageSets/Segmentation/val.txt');
fid=fopen(val_idx_file);
val_file_names=textscan(fid, '%s');
val_file_names=val_file_names{1};
fclose(fid);

train_num=length(train_file_names);
img_names=cat(1, train_file_names, val_file_names);
img_num=length(img_names);

img_files=cell(img_num, 1);
mask_files=cell(img_num, 1);

for t_idx=1:img_num
    file_name=img_names{t_idx};
    img_files{t_idx}=[file_name(1:end-4) '.jpg'];
    mask_file_name = strrep(file_name,'JPEGImages','SegmentationClass');
    mask_files{t_idx}=[mask_file_name(1:end-4) '.png'];

end

train_idxes=1:train_num;
val_idxes=train_num+1:img_num;

ds_info=[];

ds_info.img_names=img_names;
ds_info.img_files=img_files;
ds_info.mask_files=mask_files;

ds_info.train_idxes=uint32(train_idxes');
ds_info.test_idxes=uint32(val_idxes');

disp('here::')
img_dir=fullfile(ds_dir, 'JPEGImages')
mask_dir=fullfile(ds_dir, 'SegmentationClass')

data_dirs=[];
data_dirs{1}=img_dir;
data_dirs{2}=mask_dir;

data_dir_idxes_img=zeros([img_num 1], 'uint8')+1;
data_dir_idxes_mask=zeros([img_num 1], 'uint8')+2;

ds_info.data_dir_idxes_img=data_dir_idxes_img;
ds_info.data_dir_idxes_mask=data_dir_idxes_mask;
ds_info.data_dirs=data_dirs;
ds_info.ds_dir=ds_dir;

ds_info.class_info=gen_class_info_ego();

ds_info.ds_name='EgoHands';

ds_info=process_ds_info_classification(ds_info, ds_config);

end

author: aurooj

repo: Hand-Segmentation-in-the-Wild


function [ im_data ] = prepare_img( im_data, mean_data, im_size )
%PREPARE_IMG Summary of this function goes here
%   Detailed explanation goes here
IMAGE_DIM = im_size;
% Convert an image returned by Matlab's imread to im_data in caffe's data
% format: W x H x C with BGR channels
im_data = im_data(:, :, [3, 2, 1]);  % permute channels from RGB to BGR
im_data = permute(im_data, [2, 1, 3]);  % flip width and height
im_data = single(im_data);  % convert from uint8 to single
im_data = imresize(im_data, [IMAGE_DIM IMAGE_DIM], 'bilinear');  % resize im_data
im_data = im_data - mean_data;  % subtract mean_data (already in W x H x C, BGR)
end


text:
% SOSDEMO6 --- MAX CUT
% Section 3.6 of SOSTOOLS User's Manual
%

clear; echo on;
syms x1 x2 x3 x4 x5;
vartable = [x1; x2; x3; x4; x5];

% Number of cuts
f = 2.5 - 0.5*x1*x2 - 0.5*x2*x3 - 0.5*x3*x4 - 0.5*x4*x5 - 0.5*x5*x1;

% Boolean constraints
bc{1} = x1^2 - 1 ;
bc{2} = x2^2 - 1 ;
bc{3} = x3^2 - 1 ;
bc{4} = x4^2 - 1 ;
bc{5} = x5^2 - 1 ;

% =============================================
% First, initialize the sum of squares program
prog = sosprogram(vartable);

% =============================================
% Then define SOSP variables

% -- p1(x) -- : sum of squares
% Monomial vector: 5 independent variables, degree <= 1
Z = monomials(vartable,[0 1]);
[prog,p{1}] = sossosvar(prog,Z);

% -- p2(x) ... p6(x) : polynomials
% Monomial vector: 5 independent variables, degree <= 2
Z = monomials(vartable,0:2);
for i = 1:5
    [prog,p{1+i}] = sospolyvar(prog,Z);
end;

% =============================================
% Next, define SOSP constraints

% Constraint : p1(x)*(gamma - f(x)) +  p2(x)*bc1(x)
%               + ... + p6(x)*bc5(x) - (gamma-f(x))^2 >= 0
gamma = 4;

expr = p{1}*(gamma-f);
for i = 2:6
    expr = expr + p{i}*bc{i-1};
end;
expr = expr - (gamma-f)^2;

prog = sosineq(prog,expr);

% =============================================
% And call solver
prog = sossolve(prog);

% =============================================
% Program is feasible, thus 4 is an upper bound for the cut.
echo off

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function out=sech(x)

out=1./cosh(x);

author: opencobra

repo: cobratoolbox


text:
%kcparallel_signature 'Create a signature for a class for Parallelepiped classification '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cparallel_signature.pane file
%
% Parameters:
% InputFile: isample 'Input Sample', required: 'input image file with the sample'
% OutputFile: osig 'Output Signature', required: 'output file with signature'
% Toggle: useasis 'Use data as-is for bounds extraction', default: 0: 'if TRUE will use data as is for bounds extraction'
% Integer: usethreshold 'Use a threshold for bounds extraction', default: 0: 'if chosen will use a threshold for bounds extraction'
%
% Example: osig = kcparallel_signature(isample, {'isample','';'osig','';'useasis',0;'usethreshold',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cparallel_signature - Create a signature for a class for Parallelepiped classification
%
%  DESCRIPTION
% This routine creates a signature for the Parallelepiped classifier from a ROI or sample image (passed with the [-isample] parameter). Output will be written to the file specified with the [-osig] parameter.
% Signatures for the Parallelepiped classifier consist of two values per feature denoting a lower and upper threshold that determine if a pixel belongs or not to that class.
% Determination of these thresholds can be done in three modes:
% 1- With the parameter [-useasis], in this case the region will be scanned and for each element the smallest and larger values for each element will be used as signatures;
% 2- With the parameter [-usethreshold], a threshold will be used to filter the minimum values of each element - this is usually done to avoid using minimum values that are not representatives of the classes.
% 3- With the parameter [-usestddevs], a number of standard deviations will be used with the averages for each element to determine the minimum and maximum values of the signature for that element.
% The dimensions of the signature for the Parallelepiped classifier will be Ex2x1x1x1 where E is the number of element dimensions in the input data object. Since the signature has a fixed size, it should be appended in blocks with the "kappend" operator for classification with the \fIcparallel_classify\fP operator.
% Signature contents can be printed in tables with the "cparallel_printsig" operator.
%
%
%
%  EXAMPLES
% All examples for the Classify toolbox are listed on the Classify Toolbox Manual. For an example of this operator, please see the example workspace Classify:workspaces:PARALLEL.
%
%  "SEE ALSO"
% cparallel_printsig, cparallel_classify, kappend. See also the Classify Toolbox Manual for details about the Parallelepiped Classification Algorithm.
%
%  RESTRICTIONS
%
%  REFERENCES
% All references for the Classify toolbox are listed on the Classify Toolbox Manual.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

function varargout = kcparallel_signature(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kcparallel_signature(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'isample', '__input';'osig', '__output';'useasis', 0;'usethreshold', 0};
maxval={0,0,0,0};
minval={0,0,0,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Toggle','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'cparallel_signature"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
function MOT_copy_images

% test sequences
seqs = {'TUD-Crossing', 'PETS09-S2L2', 'ETH-Jelmoli', ...
    'ETH-Linthescher', 'ETH-Crossing', 'AVG-TownCentre', 'ADL-Rundle-1', ...
    'ADL-Rundle-3', 'KITTI-16', 'KITTI-19', 'Venice-1'};

indexes = {{31, 107}, {68, 111}, {82, 215}, {51, 78}, {97}, {52, 220} ...
    {232}, {183}, {90}, {281}, {235}};

for i = 1:numel(seqs)
    seq_name = seqs{i};
    index = indexes{i};

    for j = 1:numel(index)
        ind = index{j};
        file_src = sprintf('results_MOT/results_MOT_1/%s/%06d.png', seq_name, ind);
        disp(file_src);
        file_dst = sprintf('figures/%s_%06d.png', seq_name, ind);
        copyfile(file_src, file_dst);
    end

end

author: yuxng


text:
clear; clc;

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

% Add dependencies
dir_lib_caffe = './fast-rcnn/caffe-fast-rcnn/matlab/caffe/';
dir_lib_frcn = './fast-rcnn/matlab/';
dir_lib_edgebox = './edges/';
dir_lib_ptoolbox = './toolbox/';

addpath(genpath(dir_lib_caffe));
addpath(genpath(dir_lib_frcn));
addpath(genpath(dir_lib_edgebox));
addpath(genpath(dir_lib_ptoolbox));

% Set directories
dir_data = './data/';
file_input_img = [dir_data, 'input_img.jpg'];

file_model_edgebox = [dir_lib_edgebox, 'models/forest/modelBsds.mat'];
file_def_frcn = './models/fashion_detector.prototxt';
file_net_frcn = './models/fashion_detector.caffemodel';

dir_results = './results/';
file_output_bbox = [dir_results, 'bbox_img.mat'];

% Set hyper-parameters
flag_visualize = true;
size_img_max = 227;
id_category = 1; % '1' represents upper-body clothes, '2' represents lower-body clothes, '3' represents full-body clothes

% Initialize EdgeBox model
model_edgebox = load(file_model_edgebox);
model_edgebox = model_edgebox.model;
model_edgebox.opts.multiscale = 0;
model_edgebox.opts.sharpen = 2;
model_edgebox.opts.nThreads = 4;

opts_edgebox = edgeBoxes;
opts_edgebox.alpha = .65;      % step size of sliding window search
opts_edgebox.beta  = .65;      % nms threshold for object proposals
opts_edgebox.minScore = .01;   % min score of boxes to detect
opts_edgebox.maxBoxes = 1e4;   % max number of boxes to detect

% Initialize fast R-CNN model
tic;
use_gpu = true;
model_frcn = fast_rcnn_load_net(file_def_frcn, file_net_frcn, use_gpu);
disp('### Initialization completed: ');
toc;

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

% Read image
img_cur = imread(file_input_img);
height_img = size(img_cur, 1);
width_img = size(img_cur, 2);
channel_img = size(img_cur, 3);

% Resize image for speed-up
if max([height_img, width_img]) > size_img_max
	ratio = size_img_max ./ max([height_img, width_img]);
	img_resized = imresize(img_cur, ratio);
else
	ratio = 1;
	img_resized = img_cur;
end

% Generate object proposals using EdgeBox
tic;
try
	proposals_cur = edgeBoxes(img_resized, model_edgebox, opts_edgebox);
catch
	continue;
	warning('!!! error generating object proposals');
end
proposals_cur = proposals_cur(:, 1:4); % only keep coordinates
proposals_cur(:, 3:4) = proposals_cur(:, 1:2) + proposals_cur(:, 3:4); % convert [x, y, w, h] to [x1, y1, x2, y2]
proposals_cur = single(proposals_cur) + 1; % account for the 0-based indexing in EdgeBox
disp('### EdgeBox completed : ');
toc;

% Sanity check proposals
if isempty(proposals_cur)
	continue;
	warning('!!! no proposal generated');
end

% Detecting clothes using fast R-CNN
tic;
try
	bbox_pred = fast_rcnn_im_detect(model_frcn, img_resized, proposals_cur);
catch
	continue;
	warning('!!! error running fast R-CNN');
end
disp('### Fast R-CNN completed: ');
toc;

% Transform bounding box back into original coordinates w.r.t input image
bbox_img = cell(3, 1);
for id_category = 1:3
	bbox_img{id_category}(:, 1:4) = bbox_pred{id_category}(:, 1:4) ./ ratio;
	bbox_img{id_category}(:, 5) = bbox_pred{id_category}(:, 5);
end

% Visualize the most salient detection results for upper-body clothes, lower-body clothes and full-body clothes respectively
if flag_visualize
	figure(1);
	subplot(1, 3, 1); showboxes(img_cur, bbox_img{1}(1, :)); title('Upper-body Clothes');
	subplot(1, 3, 2); showboxes(img_cur, bbox_img{2}(1, :)); title('Lower-body Clothes');
	subplot(1, 3, 3); showboxes(img_cur, bbox_img{3}(1, :)); title('Full-body Clothes');
end

% Save bounding box outputs
save(file_output_bbox, 'img_cur', 'bbox_img', '-v7.3');

repo: fashion-detection


text:
function [newpt elemid weight]=proj2mesh(v,f,pt,nv,cn,radmax)
%  [newpt elemid weight]=proj2mesh(v,f,pt,nv,cn)
%
%  project a point cloud on to the surface mesh (surface can only be triangular)
%
%  author: Qianqian Fang <q.fang at neu.edu>
%  date: 12/12/2008
%
% parameters:
%      v: node coordinate of the surface mesh (nn x 3)
%      f: element list of the surface mesh (3 columns for
%            triangular mesh, 4 columns for cubic surface mesh)
%      pt: points to be projected, 3 columns for x,y and z respectively
%      nv: nodal norms (vector) calculated from nodesurfnorm.m
%          with dimensions of (size(v,1),3)
%      cn: a integer vector with the length of p, denoting the closest
%          surface nodes (indices of v) for each point in p. this
%          value can be calculated from dist2surf.m
%      radmax: if speicified, the search for elements to project will be
%          limited to those within a bounding box with half-edge-length
%          of radmax centered at the point to be projected
%
%      if nv and cn are not supplied, proj2mesh will project the point
%      cloud onto the surface by the direction pointing to the centroid
%      of the mesh
%
% outputs:
%      newpt: the projected points from p
%      elemid: a vector of length of p, denotes which surface trangle (in elem)
%             contains the projected point
%      weight: the barycentric coordinate for each projected points, these are
%             the weights
%
% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
%
% this function is part of "metch" toobox, see COPYING for license

cent=mean(v);
enum=length(f);
ec=reshape(v(f(:,1:3)',:)', [3 3,enum]);
centroid=squeeze(mean(ec,2));
newpt =zeros(size(pt,1),3);
elemid=zeros(size(pt,1),1);
weight=zeros(size(pt,1),3);
[idoldmesh,loc]=ismember(pt,v,'rows');
idnode=find(idoldmesh);
if(~isempty(idnode))
    [tt,ll]=ismember(loc(idnode),f);
    [p1,p2]=ind2sub(size(f),ll); % p1 is the index in f
    newpt(idnode,:)=pt(idnode,:);
    elemid(idnode)=p1;
    weight(sub2ind(size(weight),idnode,p2))=1;
end
radlimit=-1;

if(nargin>=5)
       % if nv and cn are supplied, use nodal norms to project the points
       direction=nv(cn,:);
       if(nargin>=6)
           radlimit=radmax;
       end
elseif(nargin==3)
       % otherwise, project toward the centroid
       direction=pt-repmat(cent,size(pt,1),1);
end

for t=1:size(pt,1)
    if(idoldmesh(t)~=0) continue; end
    maxdist=sqrt(sum((pt(t,:)-cent).*(pt(t,:)-cent)));

    if(radlimit>0)
        maxdist=radlimit;
    end

    idx=find(centroid(1,:)>pt(t,1)-maxdist & ...
             centroid(1,:)<pt(t,1)+maxdist & ...
             centroid(2,:)>pt(t,2)-maxdist & ...
             centroid(2,:)<pt(t,2)+maxdist & ...
             centroid(3,:)>pt(t,3)-maxdist & ...
             centroid(3,:)<pt(t,3)+maxdist );

    dist=centroid(:,idx);
    dist(1,:)=dist(1,:)-pt(t,1);
    dist(2,:)=dist(2,:)-pt(t,2);
    dist(3,:)=dist(3,:)-pt(t,3);
    c0=sum(dist.*dist);

    % sort the distances to accelate the calculation
    [c1,sorted]=sort(c0);

    for i=1:length(idx)

        % project the point along vector direction and calculate the intersection to a plane

        [inside,p,w]=linextriangle(pt(t,:), pt(t,:)+direction(t,:),v(f(idx(sorted(i)),:),:));

	% the intersection is within the current trianglar surface
        if(inside)
            newpt(t,:)=p;
            weight(t,:)=w;
            elemid(t)=idx(sorted(i));
            break;
        end
    end
end

author: fangq


text:
function [sSubject, iSubject] = db_add_subject( varargin )
% DB_ADD_SUBJECT: Add a subject to data base.
%
% USAGE:                        db_add_subject(sSubject, iSubject)
%        [sSubject, iSubject] = db_add_subject(sSubject)
%        [sSubject, iSubject] = db_add_subject(SubjectName)
%        [sSubject, iSubject] = db_add_subject(SubjectName, iSubject)
%        [sSubject, iSubject] = db_add_subject(SubjectName, iSubject, UseDefaultAnat, UseDefaultChannel)
%
% INPUT:
%     - SubjectName : Name of the subject
%     - sSubject    : subject structure
%     - iSubject    : Indice of the subject in the ProtocolSubjects.Subject array
%                     Set to [] to create a new subject
%     - UseDefaultAnat, UseDefaultChannel : values for new subject
% OUTPUT:
%     - sSubject : Subject structure (set to [] if an error occurs)
%     - iSubject : Subject indice in current protocol

% @=============================================================================
% 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-2011

% Get protocol description
ProtocolInfo     = bst_get('ProtocolInfo');
ProtocolSubjects = bst_get('ProtocolSubjects');
nbSubjects = length(ProtocolSubjects.Subject);

%% ===== PARSE INPUTS =====
% CALL: db_add_subject( sSubject )
if isstruct(varargin{1})
    sSubject = varargin{1};
% CALL: db_add_subject( SubjectName )
elseif ischar(varargin{1})
    % Get subject name
    SubjectName = file_standardize(varargin{1});
    % Create a new empty subject
    sSubject = db_template('Subject');
    sSubject.Name              = SubjectName;
    sSubject.FileName          = bst_fullfile(SubjectName, 'brainstormsubject.mat');
    sSubject.UseDefaultAnat    = ProtocolInfo.UseDefaultAnat;
    sSubject.UseDefaultChannel = ProtocolInfo.UseDefaultChannel;
else
    error('Invalid call to db_add_subject()');
end
% CALL: db_add_subject( ..., iSubject )
if (nargin < 2) || isempty(varargin{2})
    % New indice
    iSubject = nbSubjects + 1;
else
    iSubject = varargin{2};
end
% CALL: db_add_subject( ..., UseDefaultAnat, UseDefaultChannel )
if (nargin >= 3)
    sSubject.UseDefaultAnat = varargin{3};
end
if (nargin >= 4)
    sSubject.UseDefaultChannel = varargin{4};
end

% Is new subject ?
isNewSubject = (iSubject > nbSubjects);

%% ===== CHECK SUBJECT UNICITY =====
% Only if creating a new subject
if isNewSubject
    % Check the subject unicity
    if ~isempty(bst_get('Subject', sSubject.Name, 1))
        % A subject with the same Name is found : display an error box and return to 'Subject editor' window
        bst_error(sprintf('Subject "%s" already exists in protocol.', sSubject.Name), 'Subject editor', 0);
        sSubject = [];
        return
    end
end

%% ===== SAVE SUBJECT FILE =====
SubjectMat = db_template('subjectmat');
SubjectMat.Comments          = sSubject.Comments;
SubjectMat.UseDefaultAnat    = sSubject.UseDefaultAnat;
SubjectMat.UseDefaultChannel = sSubject.UseDefaultChannel;
try
    fullFileName = bst_fullfile(ProtocolInfo.SUBJECTS, sSubject.FileName);
    % Create target directory
    if ~file_exist(bst_fileparts(fullFileName))
        mkdir(bst_fileparts(fullFileName));
    end
    % Save file
    bst_save(fullFileName, SubjectMat, 'v7');
catch
    bst_error(sprintf('Error : cannot save subject file "%s".', sSubject.Name), 'Subject editor');
    sSubject = [];
    return
end

%% ===== UPDATE DATABASE =====
% === Update Subjects ===
% Register subject in ProtocolSubjects
ProtocolSubjects.Subject(iSubject) = sSubject;
% Update database in Matlab preferences
bst_set('ProtocolSubjects', ProtocolSubjects);

% === Create extra system conditions ===
% Only if it is a new subject
if isNewSubject
    % Add conditions: analysis_intra and default_study
    db_add_condition(sSubject.Name, bst_get('DirAnalysisIntra'), 0);
    db_add_condition(sSubject.Name, bst_get('DirDefaultStudy'),  0);
end
% Save database
db_save();

author: brainstorm-tools


text:
function p = fullfile(this)
% Returns full path to the meeg mat file
% FORMAT p = fullfile(this)
% _______________________________________________________________________
% Copyright (C) 2008-2011 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: fullfile.m 5025 2012-10-31 14:44:13Z vladimir $

p = fullfile(path(this), fname(this));

author: spm