text:
function rtk=detslip_tdcpins(rtk,cur_obs,old_obs,nav,rcv)

global glc
dt=timediff(cur_obs(1).time,old_obs(1).time);
if abs(dt)>5
    return;
end

cur_sv=satposs(cur_obs,nav,rtk.opt.sateph);
old_sv=satposs(old_obs,nav,rtk.opt.sateph);

if rcv==1
    cur_pos=rtk.ins.pos;     cur_rr=blh2xyz(cur_pos);
    old_pos=rtk.ins.oldpos;  old_rr=blh2xyz(old_pos);
elseif rcv==2
    cur_rr=rtk.basepos;  cur_pos=xyz2blh(cur_rr);
    old_rr=cur_rr;       old_pos=cur_pos;
end

n_cur=size(cur_obs,1);
n_old=size(old_obs,1);

v=zeros(n_cur,1); H=zeros(n_cur,4); P=zeros(n_cur,n_cur);
nv=0; vsat=zeros(n_cur,1); vv=zeros(n_cur,1);

for i=1:n_cur

    cur_obsi=cur_obs(i); cur_svi=cur_sv(i);
    lam=nav.lam(cur_obsi.sat,:);
    [sys,~]=satsys(cur_obsi.sat);
    if rtk.mask(sys)==0,continue;end

    if lam(1)==0,continue;end

    old_idx=0;
    for j=1:n_old
        if cur_obsi.sat==old_obs(j).sat
            old_idx=j;break;
        end
    end
    if old_idx==0,continue;end

    old_obsi=old_obs(old_idx); old_svi=old_sv(old_idx);

    if cur_svi.svh~=0||old_svi.svh~=0||norm(cur_svi.pos)<=0||norm(old_svi.pos)<=0
        continue;
    end
    if cur_obsi.L(1)==0||old_obsi.L(1)==0,continue;end

    [cur_range,cur_LOS]=geodist(cur_svi.pos,cur_rr);
    cur_azel=satazel(cur_pos,cur_LOS);
    [old_range,old_LOS]=geodist(old_svi.pos,old_rr);
    old_azel=satazel(old_pos,old_LOS); %#ok

    tdcp_obs=cur_obsi.L(1)*lam(1)-old_obsi.L(1)*lam(1);
    tdcp_pre=cur_range-old_range;

    v(nv+1,1)=tdcp_obs-tdcp_pre+glc.CLIGHT*(cur_svi.dts-old_svi.dts);

    H(nv+1,:)=[-cur_LOS,1];

    P(nv+1,nv+1)=(0.05^2/sin(cur_azel(2))^2)^-1;

    vv(i,1)=v(nv+1,1);

    nv=nv+1;

    vsat(i)=1;

end

if nv<4,return;end
if nv<n_cur
    v(nv+1:end,:)=[]; H(nv+1:end,:)=[]; P(nv+1:end,:)=[];P(:,nv+1:end)=[];
end

nv=size(v,1);exc=zeros(nv,1);
for i=1:nv
    vi=abs(v(i));
    v_other=abs(v);v_other(i)=[];
    ave_v=sum(abs(v_other))/(nv-1);
    ave_distance=sum(abs(vi-v_other))/(nv-1);
    if ave_distance>3*ave_v
        exc(i)=1;
    end
end
n=0; v0=zeros(nv,1); H0=zeros(nv,4); P0=zeros(nv,nv);
for i=1:nv
    if exc(i)==1,continue;end
    v0(n+1,1)=v(i);
    H0(n+1,:)=H(i,:);
    P0(n+1,n+1)=P(i,i);
    n=n+1;
end
if n<4,return;end
if n<nv
    v0(n+1:end,:)=[]; H0(n+1:end,:)=[]; P0(n+1:end,:)=[];P0(:,n+1:end)=[];
end

[dx,~]=least_square(v0,H0,P0);
vv=vv-dx(4);

slip=zeros(n_cur,1);
for i=1:n_cur
    sat=cur_obs(i).sat;
    if vsat(i)==0
        rtk.sat(sat).slip(1)=bitor(rtk.sat(sat).slip(1),1);
        continue;
    end

    vi=abs(vv(i));
    v_other=abs(vv);v_other(i)=[];
    ave_v=sum(abs(v_other))/(nv-1);
    ave_distance=sum(abs(vi-v_other))/(nv-1);
    if ave_distance>3*ave_v&&abs(vv(i))>0.15
        slip(i)=1;
    end
end

idx=find(slip==1);
ratio=size(idx,1)/n_cur;
if ratio>=0.5
    return;
end

for i=1:n_cur
    sat=cur_obs(i).sat;
    if slip(i)==1
        rtk.sat(sat).slip(1)=bitor(rtk.sat(sat).slip(1),slip(i));
    end
end

return


text:
function test_suite = testCos
initTestSuite;

function testTooManyInputs
assertExceptionThrown(@() cos(1, 2), 'MATLAB:maxrhs');

author: ricedsp


text:
% demo VB-learning (free learner)
% This demo simulates a VB agent that reacts to outcomes by choosing the
% action that maximizes the expected gain. Here, the agent learns the
% distribution of outcomes associated with each available action. The key
% trick is that this distribution is summmarized in terms of its two first
% moments. The agent a priori belives that these can drift over time, with
% transition variances (i.e. volatilities) exp(theta), where theta is the
% 2x1 evolution parameter vector.
% NB: there are 2 posteriro moments per moment of the outcome distribution,
% which is action-dependent. With 2 available actions, this means there are
% 2x2x2=8 hidden states in this model (to be compared with 2 Q-values for
% RL).
% NB2: the action emission law does not use any utility mapping yet!

close all
clear variables

% simulation parameters
theta = [1;-2]; % volatilities
phi = log(2); % inverse temperature = 4

f_fname = @f_VBfree;
g_fname = @g_ExpUtil;
h_fname = @h_randOutcome;

% allocate feedback struture for simulations
u0 = repmat([ones(1,50),0*ones(1,50)],1,4); % 'correct' answers
fb.inH.u0 = u0;
fb.inH.er = 1; % expected reward when correct answer
fb.inH.vr = .1; % reward variance
fb.h_fname = h_fname;
fb.indy = 1;
fb.indfb = 2;

% choose dummy initial conditions
x0 = repmat([0;1;0;1],2,1);
u = zeros(2,size(fb.inH.u0,2)+1);

n_t = length(u); % number of trials

dim = struct('n',8,'n_theta',2,'n_phi',1);

priors.muPhi = zeros(dim.n_phi,1);
priors.SigmaPhi = 1e0*eye(dim.n_phi);
priors.muTheta = [0;0];
priors.SigmaTheta = 1e0*eye(dim.n_theta);
% priors.SigmaTheta(2,2) = 0;
priors.muX0 = x0;%zeros(dim.n,1);
priors.SigmaX0 = 0e0*eye(dim.n);
priors.a_alpha = Inf;
priors.b_alpha = 0;

options.priors = priors;
options.sources = struct('type',1 ,'out', 1); % one binomial observation;
options.skipf = zeros(1,n_t);
options.skipf(1) = 1; % apply identity mapping from x0 to x1.

[y,x,x0,eta,e,u] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,Inf,Inf,options,x0,fb);

figure
plot(y-e,'r')
hold on
plot(y,'kx')
legend({'p(y=1|theta,phi,m)','binomial data samples'})
figure
ti = {'mu1: E[1st moment of u^{o}]','s1: V[1st moment of u^{o}]','mu2: E[log- 2nd moment of u^{o}]','s2: V[log- 2nd moment of u^{o}]'};
for i=1:4
    subplot(2,2,i),plot(x([i,i+4],:)'),title(ti{i})
end
drawnow
VBA_getSubplots ();

[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

displayResults(posterior,out,y,x,x0,theta,phi,Inf,Inf);

author: MBB-team

repo: VBA-toolbox


text:
function rule_test ( rule_fun, rule_title )

%*****************************************************************************80
%
%% RULE_TEST tests a rule simply by printing out some examples of it.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameter:
%
%    Input, function pointer [ X, W ] = RULE_FUN ( N );
%
%    Input, string RULE_TITLE, a description of the rule.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RULE_TEST\n' );
  fprintf ( 1, '  Compute and display points and weights of a rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rule is %s\n', rule_title );

  for n = 1 : 8
    [ x, w ] = rule_fun ( n );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Order = %d\n', n );
    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %10f  %10f\n', i, x(i), w(i) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [errRadians] = hyperSam(a, b)
% HYPERSAM Computes the spectral angle error (in radians) between two vectors
%
% Usage
%   [errRadians] = hyperSam(a, b)
% Inputs
%   a - Vector 1.
%   b - Vector 2.
% Outputs
%   errRadians - angle between vectors a and b in radians

[p,N] = size(a);
errRadians = zeros(1,N);
for k=1:N
    tmp = a(:,k);
    errRadians(k) = acos(dot(tmp, b)/ (norm(b) * norm(tmp)));
end
return;

author: isaacgerg

repo: matlabHyperspectralToolbox


text:
function g = rbfOutputGradX(model, X)

% RBFOUTPUTGRADX Evaluate derivatives of a RBF model's output with respect to inputs.
% FORMAT
% DESC returns the derivatives of the outputs of an periodic radial basis function model with
% respect to the inputs to the model. Currently a wrapper for rbfjacob.
% ARG model : the model for which the derivatives will be computed.
% ARG X : the locations at which the derivatives will be computed.
% RETURN g : the gradient of the output with respect to the inputs.
%
% SEEALSO : rbfOutputGrad, modelOutputGradX, rbfjacob
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

g = rbfjacob(model, X);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function groupedNodeStruct = findGroups(mesh,nodeList)
% Takes a list of nodes and tries to group it into unconnected unique clusters
%
%   groupedNodeStruct = findGroups(mesh,nodeList)
%
%

numVertices=length(mesh.connectionMatrix);
groupsFound=0;
currentGroup=1;

rowsToSearch=[1:numVertices];
foundAll=0;

searchConnectionMatrix=mesh.connectionMatrix(nodeList,nodeList);

counter=0;

while (~foundAll)
   nodesFoundThisTime=1;
   allNodesFound=[];
   [startNode x]=find(searchConnectionMatrix);
   currentNodes=startNode(1); % First non-zero row is the initial node

   while(nodesFoundThisTime)
      [y connectedNodes]=find(searchConnectionMatrix(currentNodes,:));
      connectedNodes=unique(connectedNodes(:));  % All the things connected to current node set
      nodesFoundThisTime=length(connectedNodes);

      if (nodesFoundThisTime)
         searchConnectionMatrix(:,currentNodes)=0;
         allNodesFound=[allNodesFound;connectedNodes];
         currentNodes=connectedNodes;

      end % if nothing found

   end % Loop while you're finding something
   groupedNodeStruct{currentGroup}.nodeList=unique(nodeList(allNodesFound));
   groupedNodeStruct{currentGroup}.tempList=unique(allNodesFound);
   currentGroup=currentGroup+1;

   foundAll=~(sum(searchConnectionMatrix(:)));

   counter=counter+1;
   disp(counter);

end

fprintf('%d groups found',counter);

% Re-generate this...
searchConnectionMatrix=mesh.connectionMatrix(nodeList,:);
searchConnectionMatrix=searchConnectionMatrix(:,nodeList);
searchConnectionMatrix=searchConnectionMatrix;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function rf=lpcao2rf(ao)
%LPCAO2RF LPC: Convert area ratios to reflection coefficients RF=(AO)

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

rf = (1-ao)./(1+ao);

author: decouples


text:
function [solution, metabolite]=consumeSomething(model,ignoreMets,isNames,minNrFluxes,params,ignoreIntBounds)
% consumeSomething
%   Tries to consume any metabolite using as few reactions as possible.
%   The intended use is when you want to make sure that you model cannot
%   consume anything without producing something. It is intended to be used
%   with no active exchange reactions.
%
%   model           a model structure
%   ignoreMets      either a cell array of metabolite IDs, a logical vector
%                   with the same number of elements as metabolites in the model,
%                   of a vector of indexes for metabolites to exclude from
%                   this analysis (opt, default [])
%   isNames         true if the supplied mets represent metabolite names
%                   (as opposed to IDs). This is a way to delete
%                   metabolites in several compartments at once without
%                   knowing the exact IDs. This only works if ignoreMets
%                   is a cell array (opt, default false)
%   minNrFluxes     solves the MILP problem of minimizing the number of
%                   fluxes instead of the sum. Slower, but can be
%                   used if the sum gives too many fluxes (opt, default
%                   false)
%   params          parameter structure as used by getMILPParams (opt)
%   ignoreIntBounds	true if internal bounds (including reversibility)
%                   should be ignored. Exchange reactions are not affected.
%                   This can be used to find unbalanced solutions which are
%                   not possible using the default constraints (opt,
%                   default false)
%
%   solution        flux vector for the solution
%   metabolite      the index of the metabolite(s) which was consumed. If
%                   possible only one metabolite is reported, but there are
%                   situations where metabolites can only be consumed in
%                   pairs (or more)
%
%   NOTE: This works by forcing at least 1 unit of "any metabolites" to be
%   consumed and then minimize for the sum of fluxes. If more than one
%   metabolite is consumed, it picks one of them to be consumed and then
%   minimizes for the sum of fluxes.
%
%   Usage: [solution, metabolite]=consumeSomething(model,ignoreMets,isNames,...
%               minNrFluxes,params,ignoreIntBounds)

if nargin<2
    ignoreMets=[];
elseif ~islogical(ignoreMets) && ~isnumeric(ignoreMets)
    ignoreMets=convertCharArray(ignoreMets);
end
if nargin<3
    isNames=false;
end
if nargin<4
    minNrFluxes=false;
end
if nargin<5
    params.relGap=0.8;
end
if nargin<6
    ignoreIntBounds=false;
end

if isNames==true && ~isempty(ignoreMets)
    %Check that metsToRemove is a cell array
    if iscellstr(ignoreMets)==false
        EM='Must supply a cell array of strings if isNames=true';
        dispEM(EM);
    end
end

if isNames==false
    indexesToIgnore=getIndexes(model,ignoreMets,'mets');
else
    indexesToIgnore=[];
    for i=1:numel(ignoreMets)
        indexesToIgnore=[indexesToIgnore;find(strcmp(ignoreMets(i),model.metNames))];
    end
end

%Change all internal reactions to be unbounded in both directions
if ignoreIntBounds==true
    [~, I]=getExchangeRxns(model);
    nonExc=true(numel(model.rxns),1);
    nonExc(I)=false;
    model=setParam(model,'lb',nonExc,-1000);
    model=setParam(model,'ub',nonExc,1000);
    model=setParam(model,'rev',nonExc,1);
end

solution=[];
metabolite=[];

nRxns=numel(model.rxns);
nMets=numel(model.mets);

%Add uptake reactions for all metabolites
model.S=[model.S speye(nMets)];

%Add so that they all consume a fake metabolite
model.S=[model.S;[sparse(1,nRxns) ones(1,nMets)*-1]];

%Change so that the ignoreMets have a coefficient 0 in their reactions.
%Does not remove the actual reaction to make mapping easier later
model.S(:,indexesToIgnore+nRxns)=0;

%Add an uptake reaction for this last fake metabolite
model.S(size(model.S,1),size(model.S,2)+1)=1;
model.b=[model.b;zeros(1,size(model.b,2))];
model.lb=[model.lb;zeros(nMets,1);1];
model.ub=[model.ub;inf(nMets+1,1)];
model.rev=[model.rev;zeros(nMets+1,1)];
model.c=zeros(size(model.S,2),1);

%Add padding to the reaction annotation to prevent an error in solveLP
padding=1:numel(model.rev);
padding=num2cell(padding)';
padding=cellfun(@num2str,padding,'uniformoutput',false);
model.rxns=padding;
model.rxnNames=padding;
model.eccodes=padding;
model.rxnMiriams=padding;
model.grRules=padding;
if isfield(model,'genes')
    model.rxnGeneMat=sparse(numel(model.rev),numel(model.genes));
end
model.subSystems=padding;
model.rxnFrom=padding;
model.rxnComps=ones(numel(model.rev),1);
model.rxnNotes=padding;
model.rxnReferences=padding;
model.rxnConfidenceScores=NaN(numel(model.rev),1);

sol=solveLP(model,1);
if any(sol.x)
    %It could be that several metabolites were consumed in order to get the
    %best solution. The setdiff is to avoid including the last fake
    %metabolite
    I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1));

    if any(I) %This should always be true
        %Change the coefficients so that only the first is consumed. This
        %is not always possible, but it is tested for since it it results
        %in more easily interpretable results

        oldS=model.S;
        foundSingle=false;
        %Test if any of the metabolites could be consumed on their own
        for i=1:numel(I)
            model.S=oldS;
            J=nRxns+1:numel(model.lb)-1;
            J(I(i))=[];
            model.S(:,J)=0;
            sol=solveLP(model);
            if any(sol.x)
                foundSingle=true;
                break;
            end
        end
        %This means that there was no metabolite which could be consumed on
        %its own. Then let all the consumeable metabolites be consumed.
        if foundSingle==false
            model.S=oldS;
        end
        if minNrFluxes==true
            %Has to add names for the rxns to prevent an error in
            %minNrFluxes
            model.rxns=cell(numel(model.lb),1);
            model.rxns(:)={'TEMP'};
            model.mets=cell(size(model.b,1),1);
            model.mets(:)={'TEMP'};
            sol=solveLP(model,3,params);
        else
            sol=solveLP(model,1);
        end
        solution=sol.x(1:nRxns);
        metabolite=find(sol.x(nRxns+1:end-1)>0.1);
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
%-----------------------------------------------------------------------
% Job saved on 27-Nov-2015 18:36:22 by cfg_util (rev $Rev: 7703 $)
% spm SPM - SPM12 (12.1)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.tools.beamforming.group.BF = '<UNDEFINED>';
matlabbatch{1}.spm.tools.beamforming.group.prefix = '';
matlabbatch{1}.spm.tools.beamforming.group.plugin.batch.batchfile = {'C:\spm12\toolbox\DAiSS\batch_group_GALA_data_sources.m'};
matlabbatch{2}.spm.tools.beamforming.group.BF(1) = cfg_dep('Group analysis: BF.mat files', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{2}.spm.tools.beamforming.group.prefix = '';
matlabbatch{2}.spm.tools.beamforming.group.plugin.GALA.iter = 3;
matlabbatch{3}.spm.tools.beamforming.group.BF(1) = cfg_dep('Group analysis: BF.mat files', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{3}.spm.tools.beamforming.group.prefix = '';
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.measure = 'lJcov';
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.spread = 4;
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.threshold = 0.01;
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.mincorr = 0.7;
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.maxsize = 50;
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.distratio1 = 4;
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.distratio2 = 2;
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.cluster.maxclust.maxclustsize = 30;
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.linkmeth = 'complete';
matlabbatch{3}.spm.tools.beamforming.group.plugin.functionalROI.similarity = 0;
matlabbatch{4}.spm.tools.beamforming.group.BF(1) = cfg_dep('Group analysis: BF.mat files', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{4}.spm.tools.beamforming.group.prefix = '';
matlabbatch{4}.spm.tools.beamforming.group.plugin.batch.batchfile = {'C:\spm12\toolbox\DAiSS\batch_group_GALA_write.m'};

author: spm


text:
function [D info] = ojw_stereo(images, P, disps, sz, options)
%OJW_STEREO  Global stereo with 2nd-order smoothness prior & occlusion model
%
%   [D info] = ojw_stereo(images, P, disps, sz, options)
%
% Generates a disparity (1/depth) map for the first input image. This
% algorithm implements a "global" stereo algorithm, with asymmetrical
% occlusion modelling, using an alpha-expansion graph cuts style approach,
% but with arbitrary disparity proposals.
%
%IN:
%   images - 1xN cell array of input images, the reference image being
%            images{1}.
%   P - 3x4xN array of projection matrices for the input images, relative
%       to the output image.
%   disps - 1xM list of disparities to sample at.
%   sz - 1x2 vector of output image dimensions: [H W].
%   options - a structure containing the following input parameters:
%       col_thresh - scalar noise parameter for data likelihood.
%       occl_const - scalar occlusion cost.
%       disp_thresh - scalar disparity threshold for smoothness prior.
%       smoothness_kernel - index denoting which smoothness kernel to use.
%                           1: truncated linear; 2: truncated quadratic.
%       lambda_l - scalar smoothness prior weight for cliques crossing
%                  segmentation boundaries.
%       lambda_h - scalar smoothness prior weight for cliques not crossing
%                  segmentation boundaries.
%       seg_params - 1x3 vector of parameters for the mean-shift
%                    over-segmentation of the reference image.
%       visibility - boolean indicating whether to employ the geometrical
%                    visbility contstraint.
%       connect - scalar neighbourhood system of the graph, 4 or 8
%                 connected.
%       max_iters - scalar number of iterations to halt after, if
%                   convergence is not achieved first.
%       converge - scalar percentage decrease in energy per iteration at
%                  which optimization stops.
%       average_over - scalar number of iterations to average over when
%                      checking convergence.
%       contract - scalar number of QPBOP iterations to do.
%       improve - scalar indicating which method to use to label unlabelled
%                 nodes. 0: QPBO-F, 1: QPBOI-F, 2: QPBO-R, 3: QPBO-L,
%                 4: QPBOI-R.
%       independent - boolean indicating whether to use independent, or
%                     merely strongly-connected, regions for improve
%                     methods 2 & 4.
%
%OUT:
%   D - HxW disparity map.
%   info - structure containing other outputs from the algorithm.

% $Id: ojw_stereo.m,v 1.5 2008/11/17 11:27:35 ojw Exp $

% Crude check for a reference image
if max(abs(P([1:6 9])-[1 0 0 0 1 0 1])) > 1e-12
    error('First image must be reference image');
end

% Initialize data arrays
R = images{1}(round(P(8))+(1:sz(1)),round(P(7))+(1:sz(2)),:);
vals.I = images(2:end);
vals.P = permute(P(:,:,2:end), [2 1 3]);
vals.sz = sz;
colors = size(R, 3);
num_in = numel(images);
Rorig = uint8(R);
if colors == 1
    Rorig = repmat(Rorig, [1 1 3]);
end
vals.R = repmat(reshape(single(R), [], colors), [2 1]);
vals.d_min = disps(end);
vals.d_step = disps(1) - vals.d_min;
vals.ndisps = numel(disps);

T = reshape(uint32(1:prod(sz)), sz);
if options.planar
    % Use 2nd order smoothness prior
    SEI = [reshape(T(1:end-2,:), 1, []) reshape(T(:,1:end-2), 1, []); ...
           reshape(T(2:end-1,:), 1, []) reshape(T(:,2:end-1), 1, []); ...
           reshape(T(3:end,:), 1, []) reshape(T(:,3:end), 1, [])];
    if options.connect == 8
        SEI = [SEI [reshape(T(1:end-2,1:end-2), 1, []) reshape(T(3:end,1:end-2), 1, []); ...
                    reshape(T(2:end-1,2:end-1), 1, []) reshape(T(2:end-1,2:end-1), 1, []); ...
                    reshape(T(3:end,3:end), 1, []) reshape(T(1:end-2,3:end), 1, [])]];
    end
else
    % Use 1st order smoothness prior
    SEI = [reshape(T(1:end-1,:), 1, []) reshape(T(:,1:end-1), 1, []); ...
           reshape(T(2:end,:), 1, []) reshape(T(:,2:end), 1, [])];
    if options.connect == 8
        SEI = [SEI [reshape(T(1:end-1,1:end-1), 1, []) reshape(T(2:end,1:end-1), 1, []); ...
                    reshape(T(2:end,2:end), 1, []) reshape(T(1:end-1,2:end), 1, [])]];
    end
end
clear T

% Initialise display
vals.show_output = options.show_output;
if vals.show_output
    vals.show_output = gcf;
    set(0, 'CurrentFigure', vals.show_output);
    subplot('Position', [0 0.5 1/3 0.5]);
    sc(R, [0 255]);
end

% Segment the image using mean shift
info.segment = vgg_segment_ms(Rorig, options.seg_params(1), options.seg_params(2), options.seg_params(3));
% Find smoothness edges which don't cross segmentation boundaries
EW = reshape(~any(diff(int32(info.segment(SEI))), 1), 1, []);
EW = EW * options.lambda_h + ~EW * options.lambda_l;
EW = EW * (num_in / ((options.connect==8) + 1));
EW = reshape(repmat(EW, [4*(1+(options.planar~=0)) 1]), [], 1);

% Set up values for ibr_fuse_depths
vals.visibility = (options.visibility ~= 0) * 1e4;
vals.improve = options.improve;
vals.contract = options.contract;
vals.independent = options.independent;
vals.compress_graph = options.compress_graph;

% Set up our robust kernels
vals.ephoto = @(F) log(2) - log(exp(sum(F .^ 2, 2)*(-1/(options.col_thresh*colors)))+1);
switch options.smoothness_kernel
    case 1
        vals.esmooth = @(F) EW .* min(abs(F), options.disp_thresh);
    case 2
        EW = EW / options.disp_thresh;
        vals.esmooth = @(F) EW .* min(F.^2, options.disp_thresh^2);
    otherwise
        error('Unknown smoothness kernel specified');
end
vals.occl_val = options.occl_const + log(2);
vals.SEI = SEI;
clear T SEI EW Rorig

if nargout > 1
    % Save parameters
    info.params.disp_thresh = options.disp_thresh;
    info.params.col_thresh = options.col_thresh;
    info.params.occl_const = options.occl_const;
    info.params.lambda_l = options.lambda_l;
    info.params.lambda_h = options.lambda_h;
end

if isnumeric(options.proposal_method) && size(options.proposal_method, 1) == 1
    % Use the proposal methods:
    for a = options.proposal_method
        switch a
            case 0
                % Ordered fronto-parallel
                [D info.samedisc_optim] = ojw_stereo_optim(vals, @(n) 3, options);
                info.samedisc_optim.D = D;
            case 1
                % SameUni (random fronto-parallel)
                [D info.sameuni_optim] = ojw_stereo_optim(vals, @(n) 1, options);
                info.sameuni_optim.D = D;
            case 2
                % SegPln (prototypical segment-based stereo proposals)
                [Dproposals info.segpln_gen] = ojw_segpln(images, P, disps, R, options);
                clear R
                Dproposals = @(n) Dproposals(:,:,mod(n-1, size(Dproposals, 3))+1);
                [D info.segpln_optim] = ojw_stereo_optim(vals, Dproposals, options);
                clear Dproposals
                info.segpln_optim.D = D;
            case 3
                % Smooth*
                Dproposals = {info.segpln_optim.D, info.sameuni_optim.D, 2, 2, 2, 2};
                Dproposals = @(n) Dproposals{mod(n-1, 6)+1};
                [D info.smooth_optim] = ojw_stereo_optim(vals, Dproposals, options);
                clear Dproposals
                info.smooth_optim.D = D;
            case 4
                % Smooth
                Dproposals = {D, 2};
                Dproposals = @(n) Dproposals{(n>1)+1};
                [D info.smooth2_optim] = ojw_stereo_optim(vals, Dproposals, options);
                clear Dproposals
                info.smooth2_optim.D = D;
        end
    end
else
    % Input fixed set of proposals
    if isnumeric(options.proposal_method)
        Dproposals = @(n) options.proposal_method(:,:,mod(n-1, size(options.proposal_method, 3))+1);
    else
        Dproposals = options.proposal_method;
    end
    [D info.udprop_optim] = ojw_stereo_optim(vals, Dproposals, options);
    clear Dproposals
    info.udprop_optim.D = D;
end
return

repo: SelfExSR


text:
function d = dcircle ( p, xc, yc, r )

%*****************************************************************************80
%
%% DCIRCLE returns the signed distance of one or more points to a circle.
%
%  Discussion:
%
%    The corresponding routine in 3D is DSPHERE.
%
%  Copyright:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,2), the point coordinates.
%
%    Input, real XC, YC, the coordinates of the center of the circle.
%
%    Input, real R, the radius of the circle.
%
%    Output, real D(NP), the signed distance of each point to the
%    circle.  The point is inside, on, or outside the circle depending
%    on whether D is negative, zero, or positive.
%
  d = sqrt ( ( p(:,1) - xc ).^2 + ( p(:,2) - yc ).^2 ) - r;

  return
end


text:
function out = putBinA(A,B,r,c)
% Function that places a 2D matrix B into a larger 2D matrix A. The
% upper-left corner of B is placed at coordinates (r,c) in A. If (r,c) are
% not given, then B is put in the upper-left corner of A, a default of
% (1,1). Some error checking is done first to make sure that A and B are 2D
% matrices, A is big enough to contain B, and the given (r,c) will not
% spill B outside of A.

%% Setup variables
[ar ac] = size(A);
[br bc] = size(B);
if nargin < 4
    r = 1; c = 1;
else
    r = round(r); c = round(c);
end

%% Error checking
% check that A and B are 2D matricies
if (numel(size(A))~=2) || (numel(size(B))~=2)
    error('The input matrices must be 2D arrays');
end
% check that A is >= B on both dims
if ar<br || ac<bc
    error('Matrix "A" must be big enough to contain matrix "B"');
end
% check that the (r,c) placement lands B entirely inside of A
if (ar < (r+br-1)) || (ac < (c+bc-1))
    error('Matrix "B" will fall outside matrix "A" with these coordinates');
end

%% Place B in A at (r,c) position
out = A;
rr = r + br - 1;
cc = c + bc - 1;
out(r:rr,c:cc) = B;

author: Sable

repo: mcbench-benchmarks


text:
function [in3] = mi32in3(mi3)
% Convert volume from cubic miles to cubic inches.
% Chad Greene 2012
in3 = mi3*254358061050000;


text:
function test10 ( dim_num, n, z, ns, sample_routine, seed_init )

%*****************************************************************************80
%
%% TEST10 tests TAU_MEASURE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST10\n' );
  fprintf ( 1, '  TAU_MEASURE computes the TAU measure of quality.\n' );
  fprintf ( 1, '  2nd moment trace measure       Tau = %14f\n', ...
    tau_measure ( dim_num, n, z, ns, sample_routine, seed_init ) );

  return
end


text:
%% Copyright (C) 2003 Willem J. Atsma
%% Copyright (C) 2014-2016, 2019 Colin B. Macdonald
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public
%% License as published by the Free Software Foundation;
%% either version 3, or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.  If not,
%% see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod  @@sym {@var{doublec} =} sym2poly (@var{p})
%% @deftypemethodx @@sym {@var{c} =} sym2poly (@var{p}, @var{x})
%% Return vector of coefficients of a symbolic polynomial.
%%
%% In the two-input form, the second argument @var{x} specifies the free
%% variable; in this case this function returns a row vector @var{c} of
%% symbolic expressions. The coefficients correspond to decreasing exponent
%% of the free variable.  Example:
%% @example
%% @group
%% syms x y
%% sym2poly(2*x^2 + 3*x - pi, x)
%%    @result{} (sym) [2  3  -π]  (1×3 matrix)
%% sym2poly(x^2 + y*x, x)
%%    @result{} (sym) [1  y  0]  (1×3 matrix)
%% @end group
%% @end example
%%
%% @strong{Warning}: Using the single-argument form, the coefficient vector
%% @var{c} is a plain numeric vector (double).  This is for compatibility
%% with the Matlab Symbolic Math Toolbox.
%% We suggest making this clear in your code by explicitly casting to @code{double},
%% as in:
%% @example
%% @group
%% syms x
%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% double(sym2poly(pi*x^3 + 3*x/2 + exp(sym(1))))
%%    @result{}     3.1416      0   1.5000   2.7183
%% @end group
%% @end example
%% You may prefer specifying @var{X} or using @code{coeffs}:
%% @example
%% @group
%% coeffs(pi*x^3 + 3*x/2 + exp(sym(1)), 'all')
%%    @result{} (sym) [π  0  3/2  ℯ]  (1×4 matrix)
%% @end group
%% @end example
%%
%% If @var{p} is not a polynomial the result has no warranty.  SymPy can
%% certainly deal with more general concepts of polynomial but we do not
%% yet expose all of that here.
%%
%% @seealso{poly2sym, @@sym/coeffs, polyval, roots}
%% @end deftypemethod

%% Created: 18 April 2003
%% Changed: 25 April 2003
%%    Removed the use of differentiate to get to coefficients - round-off
%%     errors cause problems. Now using newly created sumterms().
%% Changed: 6 May 2003
%%    Removed the attempt to use ldegree(), degree() and coeff() - results
%%     with these are inconsistent.
%% Changed: 16 April 2014
%%    Used the comment header and tests in OctSymPy, but rewrote
%%    the body (by Colin Macdonald).

function c = sym2poly(p, x)

  if ~(isscalar(p))
    error ('sym2poly: works for scalar input only');
  end

  if (nargin == 1)
    ss = findsymbols(p);
    if (length (ss) >= 2)
      error ('sym2poly: input has more than one symbol: not clear what you want me to do')
    elseif (length (ss) == 1)
      x = ss{1};
    else
      x = sym('x');
    end
    convert_to_double = true;
  elseif (nargin == 2)
    convert_to_double = false;
  else
    print_usage ();
  end

  cmd = { 'f = _ins[0]'
          'x = _ins[1]'
          'p = Poly.from_expr(f,x)'
          'return p.all_coeffs(),' };

  c2 = pycall_sympy__ (cmd, sym(p), sym(x));
  if (isempty(c2))
    error ('sym2poly: empty python output, can this happen?  A bug.')
  end

  % FIXME: should be able to convert c2 to array faster than array
  % expansion!  Particularly in the case where we just convert to
  % double anyway!
  c = sym([]);
  for j = 1:numel(c2)
    % Bug #17
    %c(j) = c2{j};
    idx.type = '()'; idx.subs = {j};
    c = subsasgn(c, idx, c2{j});
  end

  if (convert_to_double)
    c = double(c);
  end

end

%!shared x,y,a,b,c
%! syms x y a b c
%!assert (isequal (sym2poly (x^2 + 3*x - 4), [1 3 -4]))
%!assert (isequal (sym2poly (x^6 - x^3), [1 0 0 -1 0 0 0]))
%!assert (isequal (sym2poly (x^2 + 3*x - 4, x), [1 3 -4]))
%!assert (norm (sym2poly (pi*x^2 + exp(sym(1))) - [pi 0 exp(1)]) < 10*eps)
%% types
%!assert (isa (sym2poly (x^2 + 3*x - 4), 'double'))
%!assert (isa (sym2poly (x^2 + 3*x - 4, x), 'sym'))
%% tests with other vars
%!assert (isequal (sym2poly (x^2+y*x, x), [sym(1) y sym(0)]))
%!assert (isequal (sym2poly (x^2+y*x, y), [x x^2]))
%% inverse relationship
%!assert (isequal (sym2poly (poly2sym ([a b c], x), x), [a b c]))
%!assert (isequal (poly2sym (sym2poly(a*x^2 + c, x), x), a*x^2 + c))
%!assert (isequal (sym2poly (poly2sym ([1 2 3])), [1 2 3]))

%!error <more than one symbol>
%! % too many symbols for single-input
%! p = a*x^2 + 2;
%! c = sym2poly (p);

%!assert (isequal (sym2poly (sym(5)), sym(5)))

repo: octsympy


text:
function r8vec_sorted_undex_test ( )

%*****************************************************************************80
%
%% R8VEC_SORTED_UNDEX_TEST tests R8VEC_SORTED_UNDEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_SORTED_UNDEX_TEST\n' );
  fprintf ( 1, '  R8VEC_SORTED_UNDEX produces index vectors which create a sorted\n' );
  fprintf ( 1, '  list of the unique elements of a sorted R8VEC, \n' );
  fprintf ( 1, '  and a map from the original vector to the (implicit) \n' );
  fprintf ( 1, '  vector of sorted unique elements.\n' );

  x_num = 9;
  x_val = [ 11.0, 11.0, 11.0, 22.0, 22.0, 33.0, 33.0, 55.0, 55.0 ];

  r8vec_print ( x_num, x_val, '  The vector X:' );

  tol = r8_epsilon ( );
  x_unique_num = r8vec_sorted_unique_count ( x_num, x_val, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tolerance for equality is %e\n', tol );
  fprintf ( 1, '  Number of unique entries in X is %d\n', x_unique_num );

  [ undx, xdnu ] = r8vec_sorted_undex ( x_num, x_val, x_unique_num, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  UNDX can be used to list the unique elements of X\n' );
  fprintf ( 1, '  in sorted order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I  UNDX   X(UNDX)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : x_unique_num
    fprintf ( 1, '  %4d  %4d  %8f\n', i, undx(i), x_val(undx(i)) );
  end

  xu_val(1:x_unique_num) = x_val(undx(1:x_unique_num));

  fprintf ( 1, '\n' );
  fprintf ( 1, '  UNDX can be used to created XU, a copy of X\n' );
  fprintf ( 1, '  containing only the unique elements, in sorted order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I  UNDX     XU(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : x_unique_num
    fprintf ( 1, '  %4d  %4d  %8f\n', i, undx(i), xu_val(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XDNU can be used to match each element of X with one of the\n' );
  fprintf ( 1, '  unique elements\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I  XDNU    X(I)       XU(XDNU(I))\n' );
  fprintf ( 1, '\n' );

  for i = 1 : x_num
    fprintf ( 1, '  %4d  %4d  %8f  %12f\n', i, xdnu(i), x_val(i), xu_val(xdnu(i)) );
  end

  return
end


text:
%% Simulation Data
simu = simulationClass();               % Initialize Simulation Class
simu.simMechanicsFile = 'RM3.slx';      % Specify Simulink Model File
simu.rampTime = 100;                    % Wave Ramp Time [s]
simu.endTime=200;                       % Simulation End Time [s]
simu.dt = 0.1;                          % Simulation Time-Step [s]
simu.explorer = 'off';

%% Wave Information
% Regular Waves
waves = waveClass('regular');           % Initialize Wave Class and Specify Type
waves.height = 2.5;                          % Wave Height [m]
waves.period = 8;                            % Wave Period [s]

%% Body Data
% Float
body(1) = bodyClass('../../../../examples/RM3/hydroData/rm3.h5');
body(1).geometryFile = '../../../../examples/RM3/geometry/float.stl';
body(1).mass = 'equilibrium';
body(1).inertia = [20907301 21306090.66 37085481.11];

% Spar/Plate
body(2) = bodyClass('../../../../examples/RM3/hydroData/rm3.h5');
body(2).geometryFile = '../../../../examples/RM3/geometry/plate.stl';
body(2).mass = 'equilibrium';
body(2).inertia = [94419614.57 94407091.24 28542224.82];

%% PTO and Constraint Parameters
% Floating (3DOF) Joint
constraint(1) = constraintClass('Constraint1'); % Initialize Constraint Class for Constraint1
constraint(1).location = [0 0 0];                    % Constraint Location [m]

% Translational PTO
pto(1) = ptoClass('PTO1');                      % Initialize PTO Class for PTO1
pto(1).stiffness = 0;                                   % PTO Stiffness [N/m]
pto(1).damping = 1200000;                             % PTO Damping [N/(m/s)]
pto(1).location = [0 0 0];                           % PTO Location [m]

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text: RTW.TargetRegistry.getInstance('reset');

author: Firmament-Autopilot

repo: FMT-Model


text:
function [ err errFrame errTot SAbsoluteT anim ] = computeError(anim, ...
  varargin)
% Compute several SFM errors
%
% It can compute:
%  - the reprojection error
%  - the 3D error between several recovered 3D shapes and a ground truth 3D
%  shape
%
% USAGE
%  [ err errFrame errTot ] = computeError( 'reproj', isProj, varargin )
%  [ err errFrame errTot SBest ] = computeError( '3D', isProj, varargin)
%
% INPUTS
%  anim       - Animation object
%  varargin   - list of paramaters in quotes alternating with their values
%       - 'doFillP', false, will replace/create anim.P with the projection
%       matrices
%       - 'animGT'   [] the ground truth anim to compare to
%       - 'doCheckAmbiguity' [true] if true, check for ambiguities due to
%       the camera model (scale for isProj and necker reversal for ~isProj)
%       - 'checkTransform' ['camera']. Transform that can be applied to get
%       a better match. Can be 'rigid', 'rigid+scale', 'homography' or
%       'camera' (cf the alignTo function)
%
% OUTPUTS 'reproj'
%  err       - [ 1 x 3 ] array of errors:
%              if no animGT is given:
%              err = mean( errFrame, 2 );
%              if an animGT is givne:
%              err(1) = mean( errFrame{1}(3,:), 2 )
%              err(2) = mean( errFrame{2}(3,:), 2 )
%              err(3) = mean( errFrame{3} )
%  errFrame  - if no animGT is given:
%              errFrame(1,:) = sum of squared reprojected distance
%                              differences per frame
%              errFrame(2,:) = sum of reprojected distances per frame
%                              divided by the span of the object
%              errFrame(3,:) = same as first one divided by the sum of
%                              squared distances
%  errFrame  - if an animGT is given, cell object:
%              {1} [ 3 x nFrame ] Bregler, CVPR00 error (SSD)
%              The rows contain the sums over the features (and for each
%              frame) of the squared errors in X-Y, in Z, and in 3D.
%              I.e. the rows contain for each frame:
%                sum_{i=1}^nPoint  (xReconstructed_i-xOriginal_i)^2
%                                         +(yReconstructed_i-yOriginal_i)^2
%                sum_{i=1}^nPoint  (zReconstructed_i-zOriginal_i)^2
%                row1+row2
%              {2} [ 3 x nFrame ] Torresani, Bregler, CVPR01 error
%              Same as above but distances are considered (not squared
%              distances), average and not sum, and the results are
%              normalized by the span of the object
%              {3} [ 1 x nFrame ] Torresani, PAMI08 error
%              At each frame:
%              \Vert SReconstructed-SOriginal \Vert_F/
%                               \Vert SOriginal \Vert_F
%  errTot   - [3 x nPoint x nFrame] error for each point at each frame
%             as a distance (simply abs of diff, not squared)
%  SAbsoluteT - the best transformed S to match SGTAbsolute (SGT
%               transformed by its R and t)
%  anim     - the original Animation object but modified is aksed to be
%
% EXAMPLE
%
% See also Animation
%
% Vincent's Structure From Motion Toolbox      Version 3.1.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

global IS_EXIST_PDIST;
if isempty(IS_EXIST_PDIST); IS_EXIST_PDIST=exist('pdist','file'); end

[ doFillP animGT doCheckAmbiguity checkTransform ] = ...
  getPrmDflt( varargin, ...
  { 'doFillP' false 'animGT' [] ...
  'doCheckAmbiguity' true 'checkTransform' 'camera'}, 1 );

nFrame = anim.nFrame; nPoint = anim.nPoint;
method='reproj';
if ~isempty(animGT); method = '3D'; end

switch method
  case 'reproj'
    % compute the projected points
    [ W disc anim ]=generateW(anim,'doFillP', doFillP);

    % precompute stuff for error number 2
    spanW = zeros(1,nFrame);
    WIsnan=isnan(W);
    hasAnyNan=any(WIsnan(:));
    if hasAnyNan
      if IS_EXIST_PDIST
        for i = 1 : nFrame
          spanW(i) = max(pdist( anim.W(:,~WIsnan(1,:,i),i)));
        end
      else
        for i = 1 : nFrame
          tmp=anim.W(:,~WIsnan(1,:,i),i);
          spanW(i) = max(max(pdist2( tmp, tmp)));
        end
      end
    else
      if IS_EXIST_PDIST
        for i = 1 : nFrame
          spanW(i) = max(pdist( anim.W(:,:,i)));
        end
      else
        for i = 1 : nFrame
          spanW(i) = max(max(pdist2( anim.W(:,:,i), anim.W(:,:,i))));
        end
      end
    end

    % compute the different errors
    errFrame = zeros( 3, nFrame );
    tmp = W - anim.W;

    % if there is a mask, only keep certain points on certain frames
    WNorm = W;
    if ~isempty(anim.mask)
      % set to 0 the ones that are 0 in the mask
      mask = repmat(reshape(~anim.mask,1,nPoint,nFrame), [2, 1, 1]);
      tmp(mask) = 0;
      WNorm(mask(:)) = 0;
    end
    WNorm = sum(reshape(WNorm,[],nFrame).^2, 1);

    errTot = abs(tmp);
    errFrame(1,:) = sum(reshape(tmp,[],nFrame).^2,1);
    errFrame(2,:) = sum(reshape(errTot,[],nFrame),1)./spanW;
    errFrame(3,:) = sqrt(errFrame(1,:)./WNorm);

    err = mean( errFrame, 2 );
  case '3D'
    [ disc disc disc disc SAbsoluteT SGTAbsolute ] = alignTo(anim, ...
      animGT, checkTransform);

    SDiff=SAbsoluteT-SGTAbsolute;

    % Update the errors
    errTot=SDiff;

    errFrame=cell(1,3);
    % compute the first type of error
    %            {1} [ 3 x nFrame ] Bregler, CVPR00 error (SSD)
    %            The rows contain the sums over the features (and for each
    %            frame) of the squared errors in X-Y, in Z, and in 3D.
    %            I.e. the rows contain for each frame:
    %              sum_{i=1}^nPoint  (xReconstructed_i-xOriginal_i)^2
    %                                     +(yReconstructed_i-yOriginal_i)^2
    %              sum_{i=1}^nPoint  (zReconstructed_i-zOriginal_i)^2
    %                row1+row2
    errFrame{1} = zeros( 3, nFrame );
    errFrame{1}(1,:) = reshape(sum(sum( SDiff(1:2,:,:).^2, 2 ),1),1,[]);
    errFrame{1}(2,:) = reshape(sum(SDiff(3,:,:).^2, 2 ),1,[]);
    errFrame{1}(3,:) = errFrame{1}(1,:)+errFrame{1}(2,:);

    %       {2} [ 3 x nFrame ] Torresani, Bregler, CVPR01 error
    %       Same as above but distances are considered (not squared
    %       distances), average and not sum, and the results are normalized
    %       by the span of the object
    normS = sqrt(reshape(sum(sum(animGT.S.^2,1),2),1,[]));
    errFrame{2} = zeros( 3, nFrame );
    errFrame{2}(1,:) = mean(sqrt(reshape(sum(sum( SDiff(1:2,:,:).^2, ...
      2),1),1,[])));
    errFrame{2}(2,:) = mean(reshape(SDiff(3,:,:),1,[]));
    errFrame{2}(3,:) = mean(sqrt(reshape(sum(sum( SDiff(:,:,:).^2, ...
      2),1),1,[])));
    % precompute stuff for error number 2
    spanS = zeros(1,size(animGT.S,3));
    if exist('pdist','file')
      for i = 1 : size(animGT.S,3)
        spanS(i) = max(pdist( animGT.S(:,:,i)));
      end
    else
      for i = 1 : size(animGT.S,3)
        spanS(i) = max(max(pdist2( animGT.S(:,:,i), animGT.S(:,:,i))));
      end
    end
    if size(spanS,2)==1; spanS = repmat(spanS,1,nFrame); end
    errFrame{2} = bsxfun(@rdivide, errFrame{2}, spanS);

    %              {3} [ 1 x nFrame ] Torresani, PAMI08 error
    %              At each frame:
    %              \Vert SReconstructed-SOriginal \Vert_F/
    %                               \Vert SOriginal \Vert_F
    errFrame{3} = sqrt(errFrame{1}(3,i))./normS;

    % Update other errors
    err=[ mean(errFrame{1}(3,:)) mean(errFrame{2}(3,:)) ...
      mean(errFrame{3}) ];
end

author: vrabaud


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  thresholds calculation from ideal and mismatched    %
%         (both random and linear gradient)            %
%               resistive binary string                %
% code by Fabrizio Conso, university of pavia, student %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
  function [v,v_err]=trs(nbit,vref)
  vref_err=vref*0.1*randn(1);
  for i=1:nbit      %calculation of resistance values
      R(1)=1;
      R1(1)=1+0.1*0.5*randn(1)+0.01*(i); % ideal DAC
      R(i+1) =2^(i-1);
      R1(i+1)=2^(i-1)+0.1*0.5*randn(1)+0.01*(i); % mismatched DAC
  end

  for j=1:nbit        %calculation of tension values
      v(j)    =vref*sum(R(1:j))/sum(R);
      v_err(j)=vref_err*sum(R1(1:j))/sum(R1);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Output             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                %
 v=v;
 if nargout > 1
	v_err=v_err;
 end
%                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function [data_rec] = gtmexpimpute(net, data)
%GTMEXPIMPUTE Impute missing data using GTM expected values
%
%	Description
%	 DATA_REC = GTMEXPIMPUTE(NET, DATA) takes a GTM structure NET, and
%	imputes the missing values in DATA according to expectation of hidden
% variables given observed data.
%
%	See also
%	GTM, GTMEM, GTMLMEAN, GMLMODE, GMMPROB

% Copyright (c) Tommi Vatanen (2012)

% Check for consistency
errstring = consist(net, 'gtm', data);
if ~isempty(errstring)
  error(errstring);
end

data_rec = data;
missing = isnan(data_rec);

net.gmmnet.centres = rbffwd(net.rbfnet, net.X);
R = gmmpost(net.gmmnet, data);

rec = R*net.gmmnet.centres;
data_rec(missing) = rec(missing);

author: ilarinieminen

repo: SOM-Toolbox


text:
function geometry_test207 ( )

%*****************************************************************************80
%
%% TEST207 tests TRIANGLE_CONTAINS_POINT_2D_1, TRIANGLE_CONTAINS_POINT_2D_2, TRIANGLE_CONTAINS_POINT_2D_3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 June 2006
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  test_num = 7;

  p_test = [ ...
     0.25,   0.25; ...
     0.75,   0.25; ...
     1.00,   1.00; ...
    11.00,   0.50; ...
     0.00,   1.00; ...
     0.50, -10.00; ...
     0.60,   0.60 ]';
  t = [ ...
    0.0, 1.0; ...
    0.0, 0.0; ...
    1.0, 0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST207\n' );
  fprintf ( 1, '  For a triangle in 2D,\n' );
  fprintf ( 1, '  TRIANGLE_CONTAINS_POINT_2D_1 reports if a point \n' );
  fprintf ( 1, '    is inside a triangle (and doesn''t care about\n' );
  fprintf ( 1, '    the ordering of the vertices)\n' );
  fprintf ( 1, '  TRIANGLE_CONTAINS_POINT_2D_2 reports if a point \n' );
  fprintf ( 1, '    is inside a triangle (and DOES care about\n' );
  fprintf ( 1, '    the ordering of the vertices)\n' );
  fprintf ( 1, '  TRIANGLE_CONTAINS_POINT_2D_3 reports if a point \n' );
  fprintf ( 1, '    is inside a triangle (and doesn''t care about\n' );
  fprintf ( 1, '    the ordering of the vertices)\n' );

  r8mat_transpose_print ( dim_num, 3, t, '  Triangle vertices:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     X       Y     In1  In2  In3\n' );
  fprintf ( 1, '\n' );

  for j = 1 : test_num

    p(1:dim_num) = p_test(1:dim_num,j);

    inside1 = triangle_contains_point_2d_1 ( t, p );
    inside2 = triangle_contains_point_2d_2 ( t, p );
    inside3 = triangle_contains_point_2d_3 ( t, p );

    fprintf ( 1, '  %10f  %10f  %1d  %1d  %1d\n', ...
      p(1:dim_num), inside1, inside2, inside3 );

  end
%
%  Make a copy of the triangle with vertices in reverse order.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeat the test, but reverse the triangle vertex\n' );
  fprintf ( 1, '  ordering.\n' );

  for j = 1 : 3
    t2(1:2,j) = t(1:2,4-j);
  end

  r8mat_transpose_print ( dim_num, 3, t2, '  Triangle vertices (reversed):' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     X       Y     In1  In2  In3\n' );
  fprintf ( 1, '\n' );

  for j = 1 : test_num

    p(1:dim_num) = p_test(1:dim_num,j);

    inside1 = triangle_contains_point_2d_1 ( t2, p );
    inside2 = triangle_contains_point_2d_2 ( t2, p );
    inside3 = triangle_contains_point_2d_3 ( t2, p );

    fprintf ( 1, '  %10f  %10f  %1d  %1d  %1d\n', ...
      p(1:dim_num), inside1, inside2, inside3 );

  end

  return
end


text:
% Build and initialize the computational graph for regression based speech
% enhancement/dereverberation
%
function [layer, para] = Build_EnhanceNet_Regression(Data_tr, para)
para.output = 'tmp';

layer = genNetworkDereverb_Regression(para.topology);     % generate the network graph
para.preprocessing{1} = {};                     % optional preprocessing for each data stream
para.preprocessing{2} = {};
para.cost_func.layer_idx = length(layer);       % specify which layers are cost function layers
para.cost_func.layer_weight = [1];              % set the weights of each cost function layer
para.IO.inputFeature = [1 1];
para.IO.isTensor = [1 1];
para = ParseOptions2(para);

% generating the scaling factor for the input, as we will need to use a
% small constant in the logarithm. We need to make sure that the power of
% speech are larger than this constant most of the time.
scale = 1e4;        % we hard code the scale to be a constant so that all network will use the same number
scale = scale/2^16; % note that we are using int16 to store waveform samples, so need to scale down
if para.topology.useWav
    layer = InitWavScaleLayer(layer, scale);
end

if strcmpi(para.topology.RegressionNetType, 'DNN')    % if use DNN, splice the frames
    idx = ReturnLayerIdxByName(layer, 'splice');
else
    idx = ReturnLayerIdxByName(layer, 'delta'); % if use LSTM, use dynamic features
end
fft_net_length = idx(1);
fft_net = layer(1:fft_net_length);
paraTmp = para;
paraTmp.out_layer_idx = fft_net_length;
paraTmp.IO = RemoveIOStream(paraTmp.IO, [2]);
paraTmp.IO.nStream = 1;
paraTmp.IO.inputFeature = paraTmp.IO.DynamicDistortion.inputFeature(1);
if isfield(paraTmp.IO, 'fileReader'); paraTmp.IO = rmfield(paraTmp.IO, 'fileReader'); end
paraTmp.IO.fileReader(1)= paraTmp.IO.DynamicDistortion.fileReader(1);
fprintf('Generate global MVN weights for mask subnet - %s\n', datestr(now));
[layer{fft_net_length+1}.W, layer{fft_net_length+1}.b] = computeGlobalCMVN(Data_tr(1), 100, paraTmp, fft_net);
VerifyPreprocessingTree(layer(1:fft_net_length+1), Data_tr(1), paraTmp, 100);

% set weight of static, velocity, and accelration features in the MSE cost
% function.
delta_idx = ReturnLayerIdxByName(layer, 'delta');
weight_idx = delta_idx+1;
for i=1:length(weight_idx)
    if isfield(layer{weight_idx(i)}, 'W') && numel(layer{weight_idx(i)}.W) == prod(layer{weight_idx(i)}.dim)
        continue;
    end
    layer{weight_idx(i)}.W = diag([ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(1); ...
        ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(2); ...
        ones(para.topology.nFreqBin,1)*para.topology.MSECostWeightSDA(3)]);
    layer{weight_idx(i)}.b = zeros(para.topology.nFreqBin*3,1);
end

end

author: singaxiong

repo: SignalGraph


text:
function apply_recursivly(in_dir,cmd,pattern)
% apply cmd recursively to all file in a directory
%
%% Input
%
%% Output

sub_dirs = getSubDirs(in_dir);

for id = 1:length(sub_dirs)

  %disp(sub_dirs{id});
  files = dir(fullfile(sub_dirs{id},pattern));

  for i = 1:length(files)
    cmd(fullfile(sub_dirs{id},files(i).name));
  end

end

author: mtex-toolbox

repo: mtex

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


text:
function display(t)
%DISPLAY Command window display of a sptenmat.
%
%   DISPLAY(T) displays the tensor with its name.
%
%   See also SPTENMAT/DISP.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

disp(t,inputname(1));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function out=resampleImageStack(in, zoom, interpolation_method)
% Uses resample.m function on the stack on images. Resamples each of the
% image in the scatck by zoom value.
% out=resampleImageStack(in, zoom, interpolation_method)
% default interpolation_method = '4-cubic';
% out(:,:,ii)=resample(in(:,:,ii),zoom, interpolation_method)
if ~exist('interpolation_method','var')
    interpolation_method = '4-cubic';
    fprintf('Default interpolation method: %s\n',interpolation_method)
else
    fprintf('Interpolation method used: %s\n',interpolation_method)
end
nz=size(in,3);

out=zeros(zoom*size(in,1),zoom*size(in,2), size(in,3));
for ii=1:nz
    out(:,:,ii)=resample(in(:,:,ii), zoom,0,interpolation_method);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ ntri, tri ] = sphere_imp_gridfaces_3d ( maxtri, nlat, nlong )

%*****************************************************************************80
%
%% SPHERE_IMP_GRIDFACES_3D produces a grid of triangles on an implicit sphere in 3D.
%
%  Discussion:
%
%    The point numbering system is the same used in SPHERE_IMP_GRIDPOINTS_3D,
%    and that routine may be used to compute the coordinates of the points.
%
%    An implicit sphere in 3D satisfies the equation:
%
%      sum ( ( P(1:DIM_NUM) - CENTER(1:DIM_NUM) )**2 ) = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer MAXTRI, the maximum number of triangles.
%
%    Input, integer NLAT, NLONG, the number of latitude and longitude
%    lines to draw.  The latitudes do not include the North and South
%    poles, which will be included automatically, so NLAT = 5, for instance,
%    will result in points along 7 lines of latitude.
%
%    Output, integer NTRI, the number of triangles.
%
%    Output, integer TRI(3,MAXTRI), the triangle vertices.
%
  dim_num = 3;
  ntri = 0;
%
%  The first row.
%
  n = 1;

  sw = 2;
  se = sw + 1;

  s_min = 2;
  s_max = nlong + 1;

  for j = 0 : nlong - 1

    if ( ntri < maxtri )
      ntri = ntri + 1;
      tri(1:dim_num,ntri) = [ sw, se, n ]';
    end

    sw = se;

    if ( se == s_max )
      se = s_min;
    else
      se = se + 1;
    end

  end
%
%  The intermediate rows.
%
  for i = 1 : nlat

    n_max = s_max;
    n_min = s_min;

    s_max = s_max + nlong;
    s_min = s_min + nlong;

    nw = n_min;
    ne = nw + 1;
    sw = s_min;
    se = sw + 1;

    for j = 0 : nlong - 1

      if ( ntri < maxtri )
        ntri = ntri + 1;
        tri(1:dim_num,ntri) = [ sw, se, nw ]';
      end

      if ( ntri < maxtri )
        ntri = ntri + 1;
        tri(1:dim_num,ntri) = [ ne, nw, se ]';
      end

      sw = se;
      nw = ne;

      if ( se == s_max )
        se = s_min;
      else
        se = se + 1;
      end

      if ( ne == n_max )
        ne = n_min;
      else
        ne = ne + 1;
      end

    end

  end
%
%  The last row.
%
  n_max = s_max;
  n_min = s_min;

  s = n_max + 1;

  nw = n_min;
  ne = nw + 1;

  for j = 0 : nlong - 1

    if ( ntri < maxtri )
      ntri = ntri + 1;
      tri(1:dim_num,ntri) = [ ne, nw, s ]';
    end

    nw = ne;

    if ( ne == n_max )
      ne = n_min;
    else
      ne = ne + 1;
    end

  end

  return
end


text:
function [out1,out2,out3] = myginput(arg1,strpointertype)
%MYGINPUT Graphical input from mouse with custum cursor pointer.
%   [X,Y] = MYGINPUT(N) gets N points from the current axes and returns
%   the X- and Y-coordinates in length N vectors X and Y.
%
%   [X,Y] = MYGINPUT(N, POINTER) also specifies the cursor pointer, e.g.
%   'crosshair', 'arrow', 'circle' etc. See "Specifying the Figure Pointer"
%   in Matlab's documentation  to see the list of available pointers.
%
%   MYGINPUT is strictly equivalent to Matlab's original GINPUT, except
%   that a second argument specifies the cursor pointer instead of the
%   default 'fullcrosshair' pointer.
%
%   Example:
%     plot(1:2,1:2,'s');
%     hold on
%     [x,y] = myginput(1,'crosshair');
%     plot(x,y,'o');
%     hold off
%
%   MYGINPUT is copied from Matlab's GINPUT rev. 5.32.4.4.
%
%   See also GINPUT.

%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 1.02,  Date: 2006/10/24

% History:
% 2005/10/31: v1.00, first version, from GINPUT rev. 5.32.4.4.
% 2005/11/25: v1.01, line 'uisuspend' modified (for compatibility with
%                    ML7.00)
% 2006/10/24: v1.02, help text improved

out1 = []; out2 = []; out3 = []; y = [];

if nargin<1     % modified MYGINPUT
    strpointertype='crosshair'; % default GINPUT pointer
end

c = computer;
if ~strcmp(c(1:2),'PC')
   tp = get(0,'TerminalProtocol');
else
   tp = 'micro';
end

if ~strcmp(tp,'none') && ~strcmp(tp,'x') && ~strcmp(tp,'micro'),
   if nargout == 1,
      if nargin == 1,
         out1 = trmginput(arg1);
      else
         out1 = trmginput;
      end
   elseif nargout == 2 || nargout == 0,
      if nargin == 1,
         [out1,out2] = trmginput(arg1);
      else
         [out1,out2] = trmginput;
      end
      if  nargout == 0
         out1 = [ out1 out2 ];
      end
   elseif nargout == 3,
      if nargin == 1,
         [out1,out2,out3] = trmginput(arg1);
      else
         [out1,out2,out3] = trmginput;
      end
   end
else

   fig = gcf;
   figure(gcf);

   if nargin == 0
      how_many = -1;
      b = [];
   else
      how_many = arg1;
      b = [];
      if  isstr(how_many) ...
            || size(how_many,1) ~= 1 || size(how_many,2) ~= 1 ...
            || ~(fix(how_many) == how_many) ...
            || how_many < 0
         error('Requires a positive integer.')
      end
      if how_many == 0
         ptr_fig = 0;
         while(ptr_fig ~= fig)
            ptr_fig = get(0,'PointerWindow');
         end
         scrn_pt = get(0,'PointerLocation');
         loc = get(fig,'Position');
         pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
         out1 = pt(1); y = pt(2);
      elseif how_many < 0
         error('Argument must be a positive integer.')
      end
   end

   % Suspend axes functions
       %haxes = findobj(fig,'type','axes');
        state = uisuspend(fig);
   %haxes = findobj(fig,'type','axes');
   %state = uisuspend(haxes);
   pointer = get(gcf,'pointer');
   set(gcf,'pointer',strpointertype);  % modified MYGINPUT
   fig_units = get(fig,'units');
   char = 0;

   % We need to pump the event queue on unix
   % before calling WAITFORBUTTONPRESS
   drawnow

   while how_many ~= 0
      % Use no-side effect WAITFORBUTTONPRESS
      waserr = 0;
      try
	keydown = wfbp;
      catch
	waserr = 1;
      end
      if(waserr == 1)
         if(ishandle(fig))
            set(fig,'units',fig_units);
	    uirestore(state);
            error('Interrupted');
         else
            error('Interrupted by figure deletion');
         end
      end

      ptr_fig = get(0,'CurrentFigure');
      if(ptr_fig == fig)
         if keydown
            char = get(fig, 'CurrentCharacter');
            button = abs(get(fig, 'CurrentCharacter'));
            scrn_pt = get(0, 'PointerLocation');
            set(fig,'units','pixels')
            loc = get(fig, 'Position');
            pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
            set(fig,'CurrentPoint',pt);
         else
            button = get(fig, 'SelectionType');
            if strcmp(button,'open')
               button = 1;
            elseif strcmp(button,'normal')
               button = 1;
            elseif strcmp(button,'extend')
               button = 2;
            elseif strcmp(button,'alt')
               button = 3;
            else
               error('Invalid mouse selection.')
            end
         end
         pt = get(gca, 'CurrentPoint');

         how_many = how_many - 1;

         if(char == 13) % & how_many ~= 0)
            % if the return key was pressed, char will == 13,
            % and that's our signal to break out of here whether
            % or not we have collected all the requested data
            % points.
            % If this was an early breakout, don't include
            % the <Return> key info in the return arrays.
            % We will no longer count it if it's the last input.
            break;
         end

         out1 = [out1;pt(1,1)];
         y = [y;pt(1,2)];
         b = [b;button];
      end
   end

   uirestore(state);
   set(fig,'units',fig_units);

   if nargout > 1
      out2 = y;
      if nargout > 2
         out3 = b;
      end
   else
      out1 = [out1 y];
   end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key = wfbp
%WFBP   Replacement for WAITFORBUTTONPRESS that has no side effects.

fig = gcf;
current_char = [];

% Now wait for that buttonpress, and check for error conditions
waserr = 0;
try
  h=findall(fig,'type','uimenu','accel','C');   % Disabling ^C for edit menu so the only ^C is for
  set(h,'accel','');                            % interrupting the function.
  keydown = waitforbuttonpress;
  current_char = double(get(fig,'CurrentCharacter')); % Capturing the character.
  if~isempty(current_char) & (keydown == 1)           % If the character was generated by the
	  if(current_char == 3)                       % current keypress AND is ^C, set 'waserr'to 1
		  waserr = 1;                             % so that it errors out.
	  end
  end

  set(h,'accel','C');                                 % Set back the accelerator for edit menu.
catch
  waserr = 1;
end
drawnow;
if(waserr == 1)
   set(h,'accel','C');                                % Set back the accelerator if it errored out.
   error('Interrupted');
end

if nargout>0, key = keydown; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: dipum

repo: dipum-toolbox

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


text:
%Question No:5
%IDEAL LOW-PASS FILTER

function idealfilter(X,P)
f=imread(X);
[M,N]=size(f);
F=fft2(double(f));
u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H=double(D<=P);
G=H.*F;
g=real(ifft2(double(G)));
imshow(f),figure,imshow(g,[ ]);
end


text:
%% positive integer dumping
if dumpmsgpack(int8(0)) ~= uint8(0)
    error('Dumping 0 failed')
end
if any(dumpmsgpack(uint8(1)) ~= uint8(1))
    error('Dumping positive fixnum failed')
end
if any(dumpmsgpack(uint8(128)) ~= uint8([204, 128]))
    error('Dumping uint8  failed')
end
if any(dumpmsgpack(uint16(256)) ~= uint8([205, 1, 0]))
    error('Dumping uint16 failed')
end
if any(dumpmsgpack(uint32(2^16)) ~= uint8([206, 0, 1, 0, 0]))
    error('Dumping uint32 failed')
end
if any(dumpmsgpack(uint64(2^32)) ~= uint8([207, 0, 0, 0, 1, 0, 0, 0, 0]))
    error('Dumping uint64 failed')
end

%% negative integer dumping
if dumpmsgpack(int8(-1)) ~= uint8(255)
    error('Dumping negative fixnum failed')
end
if any(dumpmsgpack(int8(-128)) ~= uint8([208, 128]))
    error('Dumping int8  failed')
end
if any(dumpmsgpack(int16(-256)) ~= uint8([209, 255, 0]))
    error('Dumping int16 failed')
end
if any(dumpmsgpack(int32(-2^16)) ~= uint8([210, 255, 255, 0, 0]))
    error('Dumping int32 failed')
end
if any(dumpmsgpack(int64(-2^32)) ~= uint8([211, 255, 255, 255, 255, 0, 0, 0, 0]))
    error('Dumping int64 failed')
end

%% float dumping
if any(dumpmsgpack(single(1.5)) ~= uint8([202, 63, 192, 0, 0]))
    error('Dumping float32 failed')
end
if any(dumpmsgpack(double(1.5)) ~= uint8([203, 63, 248, 0, 0, 0, 0, 0, 0]))
    error('Dumping float64 failed')
end

%% string dumping
if any(dumpmsgpack('foo') ~= uint8([163, 102, 111, 111]))
    error('Dumping fixstr failed')
end
if any(dumpmsgpack(repmat('a', [1, 32])) ~= uint8([217, 32, ones(1, 32)*'a']))
    error('Dumping str8 failed')
end
if any(dumpmsgpack(repmat('a', [1, 2^8])) ~= uint8([218, 1, 0, ones(1, 2^8)*'a']))
    error('Dumping str16 failed')
end
if any(dumpmsgpack(repmat('a', [1, 2^16])) ~= uint8([219, 0, 1, 0, 0, ones(1, 2^16)*'a']))
    error('Dumping str16 failed')
end

%% bin dumping
if any(dumpmsgpack(repmat(uint8(42), [1, 32])) ~= uint8([196, 32, ones(1, 32)*42]))
    error('Dumping str8 failed')
end
if any(dumpmsgpack(repmat(uint8(42), [1, 2^8])) ~= uint8([197, 1, 0, ones(1, 2^8)*42]))
    error('Dumping str16 failed')
end
if any(dumpmsgpack(repmat(uint8(42), [1, 2^16])) ~= uint8([198, 0, 1, 0, 0, ones(1, 2^16)*42]))
    error('Dumping str16 failed')
end

%% array dumping
if any(dumpmsgpack({uint8(1), uint8(2)}) ~= uint8([146, 1, 2]))
    error('Dumping fixarray failed')
end
if any(dumpmsgpack(num2cell(repmat(uint8(42), [1, 16]))) ~= uint8([220, 0, 16, repmat(42, [1, 16])]))
    error('Dumping array16 failed')
end
% takes too long:
% if any(dumpmsgpack(num2cell(repmat(uint8(42), [1, 2^16]))) ~= uint8([221, 0, 1, 0, 0 repmat(42, [1, 2^16])]))
%     error('Dumping array32 failed')
% end

%% map dumping
if any(dumpmsgpack(struct('one', uint8(1), 'two', uint8(2))) ~= uint8([130, dumpmsgpack('one'), 1, dumpmsgpack('two'), 2]))
    error('Dumping fixmap failed')
end
data = struct();
msgpack = uint8([]);
for n=[1 10 11 12 13 14 15 16 2 3 4 5 6 7 8 9] % default struct field order
    data.(['x' num2str(n)]) = uint8(n);
    msgpack = [msgpack dumpmsgpack(['x' num2str(n)]) uint8(n)];
end
if any(dumpmsgpack(data) ~= uint8([222, 0, 16, msgpack]))
    error('Dumping map16 failed')
end
% map32 takes too long

author: bastibe

repo: transplant

save_path: github-repos/MATLAB/bastibe-transplant


text:
function y = powercone(z,x,y,alpha)
%POWERCONE Defines a power cone x^alpha y ^(1-alpha) > |z|
%
% Input
%    z,y,x   : sclar SDPVAR objects.
%    alpha   : scalar double 0<=alpha<=1
%
% Example
%    F = powercone(z,x,y,alpha)

if numel(z)>1
    error('x must be a scalar')
end
if numel(x)>1
    error('x must be a scalar')
end
if numel(y)>1
    error('y must be a scalar')
end

try
    y = [x;y;z;alpha];
    y.typeflag = 20;
    y = lmi(y);
catch
    rethrow(lasterror)
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [R,SS] = fit_rotations(S,varargin)
  % FIT_ROTATIONS Given an input mesh and new positions find rotations for
  % every vertex that best maps its one ring to the new one ring
  %
  % R = fit_rotations(S,'ParamName',ParamValue)
  %
  % Inputs:
  %   S  dim by dim by #rotations list of covariance matrices to fit rotations
  %     to
  %   Optional parameters
  %     'AllowFlips'  optionally followed by true or false, find best fitting
  %       rotation OR reflection
  %     'SinglePrecision'  use single precision if available.
  % Outputs:
  %   R  dim by dim by #F list of rotations
  %   SS  dim by dim by #rotations list of svd diagonals
  %

  dim = size(S,1);
  assert(dim == size(S,2));
  nr = size(S,3);
  allow_flips = false;
  use_mex = true;

  single_precision = true;
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'AllowFlips','SinglePrecision','Mex'}, ...
    {'allow_flips','single_precision','use_mex'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  % Even faster way to check if mex exists
  if use_mex && fit_rotations_mex && ~allow_flips
    nr = size(S,3);
    dim = size(S,1);
    SS = reshape(permute(S,[3 1 2]),[nr*dim dim]);
    R = fit_rotations_mex(SS,varargin{:});
    R = reshape(R,[dim dim nr]);
    return;
  end

  %R = cellfun(@fit_rotation,S,'UniformOutput',false);
  %R = reshape(cell2mat(R'),[dim dim nr]);
  % For loop is faster
  R = zeros([dim dim nr]);
  SS = zeros([dim dim nr]);
  for ii = 1:nr
    % svd
    [su,ss,sv]=svd(S(:,:,ii));
    Ri = sv*su';
    SS(:,:,ii) = ss;
    % if reflection then flip last column
    if(~allow_flips && det(Ri) < 0 )
      su(:,end) = -su(:,end);
      Ri = sv*su';
    end
    % should definitely be rotation now
    %assert( det(Ri) >= 0 );
    R(:,:,ii) = Ri;
  end

  function D = det3(M)
    % DET3 compute the determinant of a 3x3 matrix
    % Input:
    %   M  3 by 3 matrix
    % Output:
    %   D  determinant
    %
    % http://en.wikipedia.org/wiki/Determinant#3-by-3_matrices
    D = ...
      M(1,1) * M(2,2) * M(3,3) + ...
      M(1,2) * M(2,3) * M(3,1) + ...
      M(1,3) * M(2,1) * M(3,2) - ...
      M(1,3) * M(2,2) * M(3,1) - ...
      M(1,2) * M(2,1) * M(3,3) - ...
      M(1,1) * M(2,3) * M(3,2);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function [ a, R, T, T3D, params, error, shapeOrtho ] = fit_PDM_ortho_proj_to_2D_no_reg( M, E, V, shape2D)
%FITPDMTO2DSHAPE Summary of this function goes here
%   Detailed explanation goes here

    hidden = false;

    % if some of the points are unavailable modify M, V, and shape2D (can
    % later infer the actual shape from this)
    if(sum(shape2D(:)==0) > 0)

        hidden = true;
        % which indices to remove
        inds_to_rem = shape2D(:,1) == 0 | shape2D(:,2) == 0;

        shape2D = shape2D(~inds_to_rem,:);

        inds_to_rem = repmat(inds_to_rem, 3, 1);

        M_old = M;
        V_old = V;

        M = M(~inds_to_rem);
        V = V(~inds_to_rem,:);

    end

    num_points = numel(M) / 3;

    m = reshape(M, num_points, 3)';
    width_model = max(m(1,:)) - min(m(1,:));
    height_model = max(m(2,:)) - min(m(2,:));

    bounding_box = [min(shape2D(:,1)), min(shape2D(:,2)),...
                    max(shape2D(:,1)), max(shape2D(:,2))];

    a = (((bounding_box(3) - bounding_box(1)) / width_model) + ((bounding_box(4) - bounding_box(2))/ height_model)) / 2;

    tx = (bounding_box(3) + bounding_box(1))/2;
    ty = (bounding_box(4) + bounding_box(2))/2;

    % correct it so that the bounding box is just around the minimum
    % and maximum point in the initialised face
    tx = tx - a*(min(m(1,:)) + max(m(1,:)))/2;
    ty = ty - a*(min(m(2,:)) + max(m(2,:)))/2;

    R = eye(3);
    T = [tx; ty];

    params = zeros(size(E));

    currShape = getShapeOrtho(M, V, params, R, T, a);

    currError = getRMSerror(currShape, shape2D);

    reg_rigid = zeros(6,1);
    regFactor = 0.25;
    regularisations = [reg_rigid; regFactor ./ E]; % the above version, however, does not perform as well
    regularisations = diag(regularisations)*diag(regularisations);

    red_in_a_row = 0;

    for i=1:1000

        shape3D = M + V * params;
        shape3D = reshape(shape3D, numel(shape3D) / 3, 3);

        % Now find the current residual error
        currShape = a * R(1:2,:)*shape3D' + repmat(T, 1, numel(M)/3);
        currShape = currShape';

        error_res = shape2D - currShape;

        eul = Rot2Euler(R);

        p_global = [a; eul'; T];

        % get the Jacobians
        J = CalcJacobian(M, V, params, p_global);

        % RLMS style update
        p_delta = (J'*J + regularisations) \ (J'*error_res(:) - regularisations*[p_global;params]);

        % not to overshoot
        p_delta = 0.5 * p_delta;

        [params, p_global] = CalcReferenceUpdate(p_delta, params, p_global);

        a = p_global(1);
        R = Euler2Rot(p_global(2:4));
        T = p_global(5:6);

        shape3D = M + V * params;
        shape3D = reshape(shape3D, numel(shape3D) / 3, 3);
        currShape = a * R(1:2,:)*shape3D' + repmat(T, 1, numel(M)/3);
        currShape = currShape';

        error = getRMSerror(currShape, shape2D);

        if(0.999 * currError < error)
            red_in_a_row = red_in_a_row + 1;
            if(red_in_a_row == 5)
                break;
            end
        end

        currError = error;

    end

    if(hidden)
        shapeOrtho = getShapeOrtho(M_old, V_old, params, R, T, a);
    else
        shapeOrtho = currShape;
    end

    T3D = [0;0;0];

end

function [shape2D] = getShapeOrtho(M, V, p, R, T, a)

    % M - mean shape vector
    % V - eigenvectors
    % p - parameters of non-rigid shape
    % R - rotation matrix
    % T - translation vector (tx, ty)
    shape3D = getShape3D(M, V, p);
    shape2D = a * R(1:2,:)*shape3D' + repmat(T, 1, numel(M)/3);
    shape2D = shape2D';
end

function [shape2D] = getShapeOrthoFull(M, V, p, R, T, a)

    % M - mean shape vector
    % V - eigenvectors
    % p - parameters of non-rigid shape
    % R - rotation matrix
    % T - translation vector (tx, ty)
    T = [T; 0];
    shape3D = getShape3D(M, V, p);
    shape2D = a * R*shape3D' + repmat(T, 1, numel(M)/3);
    shape2D = shape2D';
end

function [shape3D] = getShape3D(M, V, params)

    shape3D = M + V * params;
    shape3D = reshape(shape3D, numel(shape3D) / 3, 3);

end

function [error] = getRMSerror(shape2Dv1, shape2Dv2)

    error = sqrt(mean(reshape(shape2Dv1 - shape2Dv2, numel(shape2Dv1), 1).^2));

end

% This calculates the combined rigid with non-rigid Jacobian
function J = CalcJacobian(M, V, p, p_global)

    n = size(M, 1)/3;

    non_rigid_modes = size(V,2);

    J = zeros(n*2, 6 + non_rigid_modes);

    % now the layour is
    % ---------- Rigid part -------------------|----Non rigid part--------|
    % dx_1/ds, dx_1/dr1, ... dx_1/dtx, dx_1/dty dx_1/dp_1 ... dx_1/dp_m
    % dx_2/ds, dx_2/dr1, ... dx_2/dtx, dx_2/dty dx_2/dp_1 ... dx_2/dp_m
    % ...
    % dx_n/ds, dx_n/dr1, ... dx_n/dtx, dx_n/dty dx_n/dp_1 ... dx_n/dp_m
    % dy_1/ds, dy_1/dr1, ... dy_1/dtx, dy_1/dty dy_1/dp_1 ... dy_1/dp_m
    % ...
    % dy_n/ds, dy_n/dr1, ... dy_n/dtx, dy_n/dty dy_n/dp_1 ... dy_n/dp_m

    % getting the rigid part
    J(:,1:6) = CalcRigidJacobian(M, V, p, p_global);

    % constructing the non-rigid part
    R = Euler2Rot(p_global(2:4));
    s = p_global(1);

    % 'rotate' and 'scale' the principal components

    % First reshape to 3D
    V_X = V(1:n,:);
    V_Y = V(n+1:2*n,:);
    V_Z = V(2*n+1:end,:);

    J_x_non_rigid = s*(R(1,1)*V_X + R(1,2)*V_Y + R(1,3)*V_Z);
    J_y_non_rigid = s*(R(2,1)*V_X + R(2,2)*V_Y + R(2,3)*V_Z);

    J(1:n, 7:end) = J_x_non_rigid;
    J(n+1:end, 7:end) = J_y_non_rigid;

end

function J = CalcRigidJacobian(M, V, p, p_global)

 	n = size(M, 1)/3;

	% Get the current 3D shape (not affected by global transform, as this
	% is how the Jacobian was derived (for derivation please see
	% ../derivations/orthoJacobian
	shape3D = GetShape3D(M, V, p);

	% Get the rotation matrix corresponding to current global orientation
	R = Euler2Rot(p_global(2:4));
	s = p_global(1);

    % Rigid Jacobian is laid out as follows
    % dx_1/ds, dx_1/dr1, dx_1/dr2, dx_1/dr3, dx_1/dtx, dx_1/dty
    % dx_2/ds, dx_2/dr1, dx_2/dr2, dx_2/dr3, dx_2/dtx, dx_2/dty
    % ...
    % dx_n/ds, dx_n/dr1, dx_n/dr2, dx_n/dr3, dx_n/dtx, dx_n/dty
    % dy_1/ds, dy_1/dr1, dy_1/dr2, dy_1/dr3, dy_1/dtx, dy_1/dty
    % ...
    % dy_n/ds, dy_n/dr1, dy_n/dr2, dy_n/dr3, dy_n/dtx, dy_n/dty

    J = zeros(n*2, 6);

    % dx/ds = X * r11  + Y * r12 + Z * r13
    % dx/dr1 =  s*(r13 * Y - r12 * Z)
    % dx/dr2 = -s*(r13 * X - r11 * Z)
    % dx/dr3 =  s*(r12 * X - r11 * Y)
    % dx/dtx = 1
    % dx/dty = 0

    % dy/ds = X * r21  + Y * r22 + Z * r23
    % dy/dr1 =  s * (r23 * Y - r22 * Z)
    % dy/dr2 = -s * (r23 * X - r21 * Z)
    % dy/dr3 =  s * (r22 * X - r21 * Y)
    % dy/dtx = 0
    % dy/dty = 1

    % set the Jacobian for x's

    % with respect to scaling factor
    J(1:n,1) = shape3D * R(1,:)';

    % with respect to angular rotation around x, y, and z axes

    % Change of x with respect to change in axis angle rotation
	dxdR = [      0,  R(1,3), -R(1,2);
            -R(1,3),       0,  R(1,1);
             R(1,2), -R(1,1),      0];

    J(1:n,2:4) = s*(dxdR * shape3D')';

    % with respect to translation
    J(1:n,5) = 1;
    J(1:n,6) = 0;

    % set the Jacobian for y's

    % with respect to scaling factor
    J(n+1:end,1) = shape3D * R(2,:)';

     % with respect to angular rotation around x, y, and z axes

    % Change of y with respect to change in axis angle rotation
    dydR = [      0,  R(2,3), -R(2,2);
            -R(2,3),       0,  R(2,1);
             R(2,2), -R(2,1),      0];

     J(n+1:end,2:4) = s*(dydR * shape3D')';

     % with respect to translation
     J(n+1:end,5) = 0;
     J(n+1:end,6) = 1;

end

% This updates the parameters based on the updates from the RLMS
function [non_rigid, rigid] = CalcReferenceUpdate(params_delta, current_non_rigid, current_global)

    rigid = zeros(6, 1);
    % Same goes for scaling and translation parameters
    rigid(1) = current_global(1) + params_delta(1);
    rigid(5) = current_global(5) + params_delta(5);
    rigid(6) = current_global(6) + params_delta(6);

    % for rotation however, we want to make sure that the rotation matrix
    % approximation we have
	% R' = [1, -wz, wy
	%       wz, 1, -wx
	%       -wy, wx, 1]
    % is a legal rotation matrix, and then we combine it with current
    % rotation (through matrix multiplication) to acquire the new rotation

	R = Euler2Rot(current_global(2:4));

    wx = params_delta(2);
    wy = params_delta(3);
    wz = params_delta(4);

    R_delta = [1, -wz, wy;
               wz, 1, -wx;
               -wy, wx, 1];

	% Make sure R_delta is orthonormal
	R_delta = OrthonormaliseRotation(R_delta);

    % Combine rotations
	R_final = R * R_delta;

	% Extract euler angle
	euler = Rot2Euler(R_final);

	rigid(2:4) = euler;

    if(length(params_delta) > 6)
        % non-rigid parameters can just be added together
        non_rigid = params_delta(7:end) +  current_non_rigid;
    else
        non_rigid = current_non_rigid;
    end

end

function R_ortho = OrthonormaliseRotation(R)

    % U * V' is basically what we want, as it's guaranteed to be
    % orthonormal
    [U, ~, V] = svd(R);

    % We also want to make sure no reflection happened

    % get the orthogonal matrix from the initial rotation matrix
    X = U*V';

    % This makes sure that the handedness is preserved and no reflection happened
    % by making sure the determinant is 1 and not -1
    W = eye(3);
    W(3,3) = det(X);
    R_ortho = U*W*V';
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function grafica2000frame
%%%%%%%%%%%%%%%%%%           SAP2000 v14           %%%%%%%%%%%%%%%%%%%%%
%programa registrado por MECAVIL ingenieria-ciencia
%bresler_lin@hotmail.com
%numera los frames
%clf
graficaTRIDIMENSIONAL

[frini frfin elmprp ielmprp conecc inud nud]=frame3dsap2000longitud;
%%%%%%%%%%%%%%%%%%%%%%
%enumeracion frames%%%
%%%%%%%%%%%%%%%%%%%%%%
M=(frfin(:,[2 3 4])+frini(:,[2 3 4]))./2;
axis('equal')
figure(1),hold on
for i=1:1:ielmprp
 text(M(i,1),M(i,2),M(i,3),[{'\color{blue}'} [,int2str( conecc(i,1) ),]],'FontSize',7)
end
hold off


text:
function varargout = chunking(varargin)

% function f = chunking(v,num)
%
% <v> is a vector
% <num> is desired length of a chunk
%
% return a cell vector of chunks.  the last vector
% may have fewer than <num> elements.
%
% example:
% isequal(chunking(1:5,3),{1:3 4:5})
%
% OR
%
% function [f,xbegin,xend] = chunking(v,num,n)
%
% <v> is a vector
% <num> is length of a chunk
% <n> is chunk number desired
%
% return the desired chunk in <f>.
% also return the beginning and ending indices associated with
% this chunk in <xbegin> and <xend>.
%
% example:
% isequal(chunking([4 2 3],2,2),3)

switch length(varargin)

case 2
  v = varargin{1};
  num = varargin{2};

  f = {};
  for p=1:ceil(length(v)/num)
    f{p} = v((p-1)*num+1 : min(length(v),p*num));
  end

  varargout = {f};

case 3
  v = varargin{1};
  num = varargin{2};
  n = varargin{3};

  xbegin = (n-1)*num+1;
  xend = min(length(v),n*num);
  f = v(xbegin:xend);

  varargout = {f xbegin xend};

end

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function varargout = sv_result(varargin)

if (nargin == 1) & ~ischar(varargin{1})  % LAUNCH GUI

  fig = openfig(mfilename,'reuse');

  % Generate a structure of handles to pass to callbacks, and store it.
  handles = guihandles(fig);
  % Create the figure-window
  handles.hPlotFigure = figure_w_normalized_uicontrolunits('Name', 'Result plot', 'NumberTitle', 'off');
  handles.hMainAxes = newplot;
  % Store the datasets internally
  handles.vResult = varargin{1}(1);
  handles.vResults = varargin{1};
  % Set up the listbox with descriptions of the stored datasets
  sList = [];
  for nCnt = 1:length(handles.vResults)
    result = handles.vResults(nCnt);
    try
      sList = [sList; cellstr(result.sComment)];
    catch
      sList = [sList; cellstr(num2str(nCnt))];
    end
  end
  set(handles.lstSets, 'String', sList);
  set(handles.lstSets, 'Value', 1);

  guidata(fig, handles);
  % Plot the very first dataset and map
  plotit(handles)

  try
    if length(varargin{1}) > 1
      set(handles.lstMaps, 'String', varargin{1}(1).vcsGridNames);
    else
      set(handles.lstMaps, 'String', varargin{1}.vcsGridNames);
    end
  catch
  end

  if nargout > 0
    varargout{1} = fig;
  end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

  try
    [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
  catch
    disp(lasterr);
  end

end

% --------------------------------------------------------------------
function varargout = btnClose_Callback(h, eventdata, handles, varargin)

% Closes plot-figure and this figure
try
  delete(handles.hPlotFigure);
catch
end
delete(handles.dlgResult);

% --------------------------------------------------------------------
function plotit(handles)

try
  % Get selected dataset
  handles.vResult = handles.vResults(get(handles.lstSets, 'Value'));

  % Create new values for plotting
  nIdx = get(handles.lstMaps, 'Value');
  vTmp = ones(length(handles.vResult.vX) * length(handles.vResult.vY), 1) * nan;
  vTmp(handles.vResult.vUsedNodes) = handles.vResult.mValueGrid(:,nIdx);
  handles.mPlotValues = reshape(vTmp, length(handles.vResult.vY), length(handles.vResult.vX));

  % Renew plot
  figure_w_normalized_uicontrolunits(handles.hPlotFigure);
  delete(handles.hMainAxes);
  handles.hMainAxes = newplot;
  set(handles.hPlotFigure, 'renderer', 'zbuffer');

  % Plot data
  handles.hSurface = surface(handles.vResult.vX, handles.vResult.vY, zeros(size(handles.mPlotValues)), handles.mPlotValues, 'Parent', handles.hMainAxes);
  axis(handles.hMainAxes, 'image');
  if get(handles.chkInterpolated, 'Value') == 1
    shading(handles.hMainAxes, 'interp');
  else
      set(handles.hSurface,'Edgecolor','none');
  end
  set(handles.hMainAxes, 'box', 'on');
  hold on;

  % Plotting coastline & faults
  if handles.vResult.bMap
    if get(handles.chkCoastlines, 'Value') == 1
      try
        handles.hCoastline = plot(handles.vResult.vCoastline(:,1), handles.vResult.vCoastline(:,2),'LineWidth',1.0,'Color',[0  0  0 ]);
      catch
      end
    end
    if get(handles.chkFaults, 'Value') == 1
      try
        handles.hFaults = plot(handles.vResult.vFaults(:,1), handles.vResult.vFaults(:,2),'k','LineWidth',0.2);
      catch
      end
    end
  end

  if get(handles.chkSeismicity, 'Value') == 1
    if ~isempty(handles.vResult.mCatalog)
      if handles.vResult.bMap
        handles.hEQPlot = plot(handles.vResult.mCatalog(:,1),handles.vResult.mCatalog(:,2),'.k','MarkerSize',6,'Marker','.','Color','k','Visible','on');
      else
        handles.hEQPlot = plot(handles.vResult.mCatalog(:,length(handles.vResult.mCatalog(1,:))),-handles.vResult.mCatalog(:,7),'.k','MarkerSize',6,'Marker','.','Color','k','Visible','on');
      end
    end
  end
  % Set information values
  set(handles.lblMin, 'String', num2str(min(min(handles.mPlotValues)), 3));
  set(handles.lblMax, 'String', num2str(max(max(handles.mPlotValues)), 3));
  set(handles.lblMean, 'String', num2str(mean(handles.vResult.mValueGrid(:,nIdx),'omitnan'), 3));
  if get(handles.chkKeepColormap, 'Value') == 0
    set(handles.txtMin, 'String', num2str(min(min(handles.mPlotValues)), 3));
    set(handles.txtMax, 'String', num2str(max(max(handles.mPlotValues)), 3));
  end

  % Set limits
  fLimMin = str2double(get(handles.txtMin, 'String'));
  if (isnan(fLimMin) | isinf(fLimMin))
    fLimMin = min(min(handles.mPlotValues));
    if isnan(fLimMin)
      fLimMin = 0;
    elseif isinf(fLimMin) % Set to the lowest non-inf value
      mSelection = ~isinf(handles.mPlotValues);
      mNoInfValues = handles.mPlotValues(mSelection);
      fLimMin = min(min(mNoInfValues));
    end
    set(handles.txtMin, 'String', num2str(fLimMin, 3));
  end
  fLimMax = str2double(get(handles.txtMax, 'String'));
  if (isnan(fLimMax) | isinf(fLimMax))
    fLimMax = max(max(handles.mPlotValues));
    if isnan(fLimMax)
      fLimMax = 1;
    elseif isinf(fLimMax) % Set to the highest non-inf value
      mSelection = ~isinf(handles.mPlotValues);
      mNoInfValues = handles.mPlotValues(mSelection);
      fLimMax = max(max(mNoInfValues));
    end
    set(handles.txtMax, 'String', num2str(fLimMax, 3));
  end
  % If both of the values have the same value, colorbar is causing trouble
  if fLimMax == fLimMin
    fLimMax = fLimMax + 0.001;
  end
  set(handles.hMainAxes, 'CLim', [fLimMin fLimMax]);

  % Create a new colorbar
  nColormap = get(handles.lstColormap, 'Value');
  switch nColormap
  case  1, vColormap = autumn(256);
  case  2, vColormap = bone(256);
  case  3, vColormap = colorcube(256);
  case  4, vColormap = cool(256);
  case  5, vColormap = copper(256);
  case  6, vColormap = flag(256);
  case  7, vColormap = gray(256);
  case  8, vColormap = hot(256);
  case  9, vColormap = hsv(256);
  case 10, vColormap = jet(256);
  case 11, vColormap = lines(256);
  case 12, vColormap = pink(256);
  case 13, vColormap = prism(256);
  case 14, vColormap = spring(256);
  case 15, vColormap = summer(256);
  case 16, vColormap = white(256);
  case 17, vColormap = winter(256);
  case 18, vColormap = gui_Colormap_HotCut(256);
  case 19, vColormap = gui_Colormap_HSVCut(256);
  case 20, vColormap = gui_Colormap_RedGrayGreen(256);
  case 21, vColormap = gui_Colormap_RedGrayBlue(256);
  case 22, vColormap = gui_Colormap_BlueYellowRed(256);
  case 23, vColormap = gui_Colormap_Rastafari(256);
  case 24, vColormap = gui_Colormap_4Colors(256);
  case 25, vColormap = gui_Colormap_Skyline(256);
  case 26, vColormap = gui_Colormap_6Colors(256);
  otherwise, vColormap = jet(256);
  end
  % Flip colormap
  if get(handles.chkFlipColormap, 'Value') == 1
    vColormap = flipud(vColormap);
  end
  if get(handles.chkColorscale, 'Value') == 0
    nLimit = str2double(get(handles.txtScaleLimit, 'String'));
    if (isnan(nLimit)) | (nLimit ~= round(nLimit))
      nLimit = 20;
      set(handles.txtScaleLimit, 'String', num2str(nLimit));
    end
    vColormap = ex_dcolor(vColormap, nLimit);
  end
  colormap(handles.hMainAxes, vColormap);
  handles.hColorBar = colorbar('horiz', 'peer', handles.hMainAxes);

  % Save changes
  guidata(handles.dlgResult, handles);
catch
end

% --------------------------------------------------------------------
function varargout = mnuFileLoad_Callback(h, eventdata, handles, varargin)

% Get the filename of new result file
[bOK, sFilePath] = ex_GetFile('*.mat', 'Select a result file');
if bOK
  % Load the new result records
  vTmpResults = load(sFilePath, 'vResults');
  % Add the new result records
  for nCnt = 1:length(vTmpResults.vResults)
    result=vTmpResults.vResults(nCnt);
    handles.vResults = [handles.vResults; result];
  end
  % Create a new list of result records
  sList = [];
  for nCnt = 1:length(handles.vResults)
    result = handles.vResults(nCnt);
    try
      sList = [sList; cellstr(result.sComment)];
    catch
      sList = [sList; cellstr(num2str(nCnt))];
    end
  end
  set(handles.lstSets, 'String', sList);
  set(handles.lstSets, 'Value', 1);
  % Store the new result records
  guidata(handles.dlgResult, handles);
  % Plot the very first dataset and map
  plotit(handles)
  try
    set(handles.lstMaps, 'String', result.vcsGridNames);
  catch
  end
end

% --------------------------------------------------------------------
function varargout = mnuFileSave_Callback(h, eventdata, handles, varargin)

% Save the data
[newfile, newpath] = uiputfile('result.mat', 'Save resutls');
newfile = [newpath newfile];
if length(newfile) > 1
  vResults = handles.vResults;
  save(newfile);
end

% --------------------------------------------------------------------
function varargout = mnuFMD_Callback(h, eventdata, handles, varargin)

aux_McEMR(handles.vResult, handles.dlgResult);

% --------------------------------------------------------------------
function varargout = update_Callback(h, eventdata, handles, varargin)

plotit(handles);

% Buttons
% --------------------------------------------------------------------
function varargout = btnAnalyze_Callback(h, eventdata, handles, varargin)

sAuxiliaryFunction = get(handles.txtAuxiliaryFunction, 'String');
if ~isempty(sAuxiliaryFunction)
  eval([sAuxiliaryFunction '(handles.vResult, handles.dlgResult)']);
end

% --------------------------------------------------------------------
function varargout = btnPlotCDF_Callback(h, eventdata, handles, varargin)

try
  nIdx = get(handles.lstMaps, 'Value');
  vDistribution = handles.vResult.mValueGrid(:,nIdx);
  plot_CDF(vDistribution);
  assignin('base', 'vDistribution', vDistribution);
catch
end

% --------------------------------------------------------------------
function varargout = btnPlotPDF_Callback(h, eventdata, handles, varargin)

try
  nIdx = get(handles.lstMaps, 'Value');
  vDistribution = handles.vResult.mValueGrid(:,nIdx);
  vSel = ~isnan(vDistribution);
  vDistribution = vDistribution(vSel,:);
  nNumberBins = str2double(get(handles.txtNumberBins, 'String'));
  figure;
  histogram(vDistribution, nNumberBins);
catch
end

%% Getting handles
% --------------------------------------------------------------------
function varargout = GetAxesHandle(h, eventdata, handles, varargin)

varargout{1} = handles.hMainAxes;

% --------------------------------------------------------------------
function varargout = GetColorbarHandle(h, eventdata, handles, varargin)

varargout{1} = handles.hColorBar;

% --------------------------------------------------------------------
function varargout = GetFigureHandle(h, eventdata, handles, varargin)

varargout{1} = handles.hPlotFigure;

% --------------------------------------------------------------------
function varargout = GetSurfaceHandle(h, eventdata, handles, varargin)

varargout{1} = handles.hSurface;

% --------------------------------------------------------------------
function varargout = GetEQPlotHandle(h, eventdata, handles, varargin)

varargout{1} = handles.hEQPlot;

author: CelsoReyes


text:
function res=rot2eulZYX(R)
% calculates ZYX euler rotation angles from
% the rotation matrix R

% Matlab implementation of the method descirped in:
% "Computing Euler angles from a rotation matrix"
% Author: Gregory G. Slabaugh

% Copyright, Mohammd SAFEEA, 9th August 2018

% test:
% R=rotz(45)*roty(22)*rotx(30);
% res=rot2eulZYX(R)*180/pi

% res=rot2eulZYX(eye(3))

    [res]=eulerAngles(R) ;
end

function y=closeEnough(a,b)
    epsilon=0.000001;
    if (epsilon>abs(a-b))
        y=1;
    else
        y=0;
    end
end

function [res]=eulerAngles(R)
    %% check for gimbal lock
    R=R';
    if (closeEnough(R(1,3), -1.0))
        b = 0; %% gimbal lock, value of x doesn't matter
        c = pi / 2;
        a = b + atan2(R(2,1), R(3,1));
    elseif (closeEnough(R(1,3), 1.0))
        b = 0;
        c = -pi / 2;
        a = -b + atan2(-R(2,1), -R(3,1));
    else  %% two solutions exist
        b1 = -asin(R(1,3));
        b2 = pi - b1;

        c1 = atan2(R(2,3) / cos(b1), R(3,3) / cos(b1));
        c2 = atan2(R(2,3) / cos(b2), R(3,3) / cos(b2));

        a1 = atan2(R(1,2) / cos(b1), R(1,1)/ cos(b1));
        a2 = atan2(R(1,2) / cos(b2), R(1,1) / cos(b2));

        %% choose one solution to return
        %% for example the "shortest" rotation
        condition=((abs(b1) + abs(c1) + abs(a1)) <= (abs(b2) + abs(c2) + abs(a2)));
        if condition
            a=a1;
            b=b1;
            c=c1;
        else
            a=a2;
            b=b2;
            c=c2;
        end
    end
    res=[a,b,c];
end

repo: KST-Kuka-Sunrise-Toolbox