text:
function hypograph( varargin )

%HYPOGRAPH Declares a hypograph variable.
%   HYPOGRAPH VARIABLE x
%   where x is a valid MATLAB variable name, declares a scalar
%   variable for the current cvx problem, and specifies it as
%   the objective of a maximization. There is no need to declare
%   a separate MAXIMIZE command.
%
%   This keyword should be used only when attempting to
%   create a new function for the CVX library. Suppose you have a
%   function F(X) whose hypograph {(X,Y)}{F(X)>=Y} can be
%   described as a CVX feasibiliy problem. Then declaring Y
%   as an EPIGRAPH turns the model into one that computes the
%   value of the function itself, maximizing over Y.
%
%   Other uses of this keyword are NOT supported and can lead
%   to numerical errors.
%
%   HYPOGRAPH VARIABLE x(n1,n2,...,nk)
%   declares a vector, matrix, or array hypograph variable with
%   dimensions n1, n2, ..., nk, each of which must be positive
%   integers.
%
%   Structure modifiers such as "symmetric", "toeplitz", etc.
%   are permitted with hypograph variables.
%
%   Examples:
%      hypograph variable x
%      hypograph variable x(100)
%
%   See also VARIABLE, EPIGRAPH.

if nargin < 2,
    error( 'Incorrect syntax for HYPOGRAPH VARIABLE. Type HELP HYPOGRAPH for details.' );
elseif ~iscellstr( varargin ),
    error( 'All arguments must be strings.' );
elseif ~strcmp( varargin{1}, 'variable' ),
    error( 'Incorrect syntax for HYPOGRAPH VARIABLE. Type HELP HYPOGRAPH for details.' );
end
evalin( 'caller', sprintf( '%s ', 'variable', varargin{2:end}, ' hypograph_' ) );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% testing split_sentence()

inputtext = 'Hi There. How Are you? Jane doesnt have a strange cup. Ms. Brown doesn''t own a little camera. She doesn''t have a light pencil.';

n = split_sentence(inputtext);

for i= 1:size(n,2)
    text = sprintf('%s\n', n{i});
    disp(text)

end

% Reading corpora text
addpath(genpath('corpora\'));
fid1 = fopen('declaration of independence.txt','r'); %# open csv file for reading

while ~feof(fid1)
    line = fgets(fid1); %# read line by line
    n = split_sentence(line);

    for i= 1:size(n,2)
        text = sprintf('%s\n', n{i});
        disp(text)

    end
end
fclose(fid1);

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function y=lagrange(x,pointx,pointy)
%
%LAGRANGE   approx a point-defined function using the Lagrange polynomial interpolation
%
%      LAGRANGE(X,POINTX,POINTY) approx the function definited by the points:
%      P1=(POINTX(1),POINTY(1)), P2=(POINTX(2),POINTY(2)), ..., PN(POINTX(N),POINTY(N))
%      and calculate it in each elements of X
%
%      If POINTX and POINTY have different number of elements the function will return the NaN value
%
%      function wrote by: Calzino
%      7-oct-2001
%
n=size(pointx,2);
L=ones(n,size(x,2));
if (size(pointx,2)~=size(pointy,2))
   fprintf(1,'\nERROR!\nPOINTX and POINTY must have the same number of elements\n');
   y=NaN;
else
   for i=1:n
      for j=1:n
         if (i~=j)
            L(i,:)=L(i,:).*(x-pointx(j))/(pointx(i)-pointx(j));
         end
      end
   end
   y=0;
   for i=1:n
      y=y+pointy(i)*L(i,:);
   end
end

author: Sable

repo: mcbench-benchmarks


text:
function [tri, bc] = gridDataFast3D(x, y, z, xi, yi, zi)
% Copyright (c) 2012, Chao Huang
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% enforce x, y, z, xi, yi, zi to be column vectors
x = x(:);
y = y(:);
z = z(:);
xi = xi(:);
yi = yi(:);
zi = zi(:);

% check if the MATLAB version is new enough to have DelaunayTri
if exist('DelaunayTri') %#ok<EXIST>
    % triangulize the data
    tri = DelaunayTri(x, y, z);
else
    % give an error
    error('''CartInterp'' set to ''linear'' for 3D simulations is not supported in your version of MATLAB. Try running the simulation with the optional input ''CartInterp'', ''nearest'' or using a different data type.');
end

% catch trinagulation error
if isempty(tri)
    error('Data cannot be triangulated.');
end

% find the nearest triangle and the corresponding Barycentric coordinates
[t, bc] = pointLocation(tri,[xi yi zi]);

% check points are valid
if any(isnan(t))
    error('Cartesian points must lie within the k-space grid defined by kgrid');
end

tri = tri(t,:);

repo: Aero-matlab


text:
close all;
clear;
clc;

%% load dataset
%load('../dataset/ORL_Face_img.mat');
%load('../dataset/Brodatz_texture_img_small_set.mat');
%load('../dataset/AR_Face_img_27x20.mat');
%load('../dataset/AR_Face_img_60x43.mat');
load('../dataset/USPS.mat');
%load('../dataset/MNIST.mat');
%load('../dataset/COIL20.mat');
%load('../dataset/COIL100.mat');

%% reduce dataset for a quick test if necessary
max_class_num = 10;
max_train_samples = 2;
max_test_samples = 15;
[TrainSet, TestSet, train_num, test_num, class_num] = reduce_dataset(TrainSet, TestSet, max_class_num, max_train_samples, max_test_samples);

%% set paramters
eigenface_flag = true;
eigenface_dim = floor(train_num/2); % example
dim = size(TrainSet.X, 1);
if eigenface_dim > dim
    eigenface_dim = dim;
end

%% normalize dataset
[TrainSet_normalized.X, TrainSet_normalized.y] = data_normalization(TrainSet.X, TrainSet.y, 'std');
[TestSet_normalized.X, TestSet_normalized.y] = data_normalization(TestSet.X, TestSet.y, 'std');

%% SVM
options.verbose = true;
options.eigenface = eigenface_flag;
options.eigenface_dim = eigenface_dim;
accuracy_svm = svm_classifier(TrainSet, TestSet, train_num, test_num, class_num, options);
fprintf('# SVM: Accuracy = %5.5f\n', accuracy_svm);

%% LSR
lambda = 0.001;
options.verbose = true;
[accuracy_lsr, ~, ~] = lsr(TrainSet, TestSet, train_num, test_num, class_num, lambda, options);
fprintf('# LSR: Accuracy = %5.5f\n', accuracy_lsr);

%% LRC
clear options;
options.verbose = true;
accuracy_lrc = lrc(TrainSet_normalized, TestSet_normalized, test_num, class_num, options);
fprintf('# LRC: Accuracy = %5.5f\n', accuracy_lrc);

%% LDRC
clear options;
options.verbose = true;
accuracy_ldrc = ldrc(TrainSet_normalized, TestSet_normalized, train_num, test_num, class_num, eigenface_dim, options);
fprintf('# LDRC: Accuracy = %5.5f\n', accuracy_ldrc);

%% LCDRC
clear options;
options.verbose = true;
accuracy_lcdrc = lcdrc(TrainSet_normalized, TestSet_normalized, train_num, test_num, class_num, eigenface_dim, options);
fprintf('# LCDRC: Accuracy = %5.5f\n', accuracy_lcdrc);

%% display accuracy
fprintf('\n\n## Summary of results\n\n')
fprintf('# SVM: Accuracy = %5.5f\n', accuracy_svm);
fprintf('# LSR: Accuracy = %5.5f\n', accuracy_lsr);
fprintf('# LRC: Accuracy = %5.5f\n', accuracy_lrc);
fprintf('# LDRC: Accuracy = %5.5f\n', accuracy_ldrc);
fprintf('# LCDRC: Accuracy = %5.5f\n', accuracy_lcdrc);

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [intT,T]=ChebyshevPolyInt(tau,n,tauStart,tauEnd)
%%CHEBYSHEVPOLYDERIV Compute the indefinite integral (with zero additive
%                    constant) of Chebyshev polynomials of the first kind
%                    from order 0 to order n evaluated at the points given
%                    in tau.
%
%INPUTS: tau An NX1 or 1XN vector of values from tauStart to tauEnd, or
%            from -1 to 1 if tauStart and tauEnd are omitted, where one
%            wishes to evaluate the indefinite integrals (antiderivatives)
%            of the Chebyshev polynomials.
%          n The non-negative integer maximum order of the Chebyshev
%            polynomials evaluated.
% tauStart,tauEnd The possible range of the inputs. If omitted, a range of
%            -1 to 1 is assumed --the normal range for Chebyshev
%            polynomials. The option for mapping to a wider range is useful
%            when using Chebyshev polynomials for interpolation. Note that
%            such polynomials are generally not useful for interpolating
%            far outside of the valid range.
%
%OUTPUTS: dT An (n+1)XN matrix of the antiderivatives of the Chebyshev
%            polynomials from order 0 to n evaluated at each of the values
%            of tau. T(i,j) is the antiderivative of the (i-1)th order
%            Chebyshev polynomial evaluated at tau(j), taking into account
%            that the function has been mapped to the range
%            (tauStart,tauEnd), if necessary.
%          T Since the Chebyshev functions have to be computed to find the
%            indefinite integrals, an (n+2)XN matrix of the Chebyshev
%            functions can also be returned, if desired
%
%The recursion for finding the indefinite integral values can be obtained
%by examining how the function ChebyshevPolyIntCoeffs modifies the
%coefficients of a weighted series of Chebyshev polynomials when only a
%single 1 is present for each given order.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=length(tau);
tau=tau(:)';%Make tau into a column vector.
%Map the tau values to the -1 to 1 range, if necessary.
if(nargin>2)
    tau=(tau-0.5*(tauStart+tauEnd))/(0.5*(tauEnd-tauStart));
end

%Get the Chebyshev polynomial function values.
T=ChebyshevPoly(tau,n+1);
intT=zeros(n+1,N);

%Deal with the n=0 case.
curN=0;
intT(curN+1,:)=T(curN+1+1,:);

if(n>0)
    %Deal with the n=1 case.
    curN=1;
    intT(curN+1,:)=0.5*(T(curN+1+1,:)/(curN+1)+0.5);
end

%Perform recusion for the rest of the values.
for curN=2:n
    intT(curN+1,:)=0.5*(T(curN+1+1,:)/(curN+1)-T(curN-1+1,:)/(curN-1));
end

%The constant to handle a mapping from -1->1 to some other range of
%parameters.
if(nargin>2)
    intT=intT*((tauEnd-tauStart)/2);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function model = fgplvmAddConstraint(model,options,varargin)

% FGPLVMADDCONSTRAINT Add latent constraints to FGPLVM model
% FORMAT
% DESC Adds constraint structure to FGPLVM model
% ARG model : fgplvm model
% ARG options : options strucure as returnded from
% constraintOptions
% RETURN model : the GP-LVM model.
%
% COPYRIGHT : Carl Henrik Ek, 2009
%
% SEEALSO : constraintOptions

% FGPLVM

%
if(isfield(model,'constraints')&&~isempty(model.constraints))
  % constraints init
  model = constraintCreate(model,[],options,varargin{:});
else
  % no constraints
  model.constraints = [];
  model.constraints.q = model.q;
  model.constraints.N = model.N;
  model.constraints.id = [];
  model.constraints.numConstraints = 0;
  model.constraints.comp = {};
  model = constraintCreate(model,[],options,varargin{:});
end

model.constraints.numConstraints = model.constraints.numConstraints + 1;
model.constraints.id = [model.constraints.id; false*ones(1,model.q)];
if(~isempty(options.dim))
  model.constraints.id(end,options.dim) = true;
else
  warning('No Dimensions set for constraints, applying to full latent space');
  model.constraints.id(end,:) = true;
end

% update constraints
for(i = 1:1:model.constraints.numConstraints)
  model.constraints.comp{i} = constraintExpandParam(model.constraints.comp{i},model.X);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [x,w,P,Z] = fnnlsb(XtX,Xty,P_old,Z_old,tol)
%FNNLSb	Non-negative least-squares.
%
%     	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	%							 %
%	% Please note this version of NNLS is for advanced users %
%	% Please refer to the m-file FNNLS for a simpler fast    %
%	% version of NNLS				         %
%	%						         %
%     	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% 	Adapted from NNLS of Mathworks, Inc.
%
%	x = fnnlsb(A,b) returns the vector X that solves x = pinv(A)*b
%	in a least squares sense, subject to x >= 0.
%
%	A default tolerance of TOL = MAX(SIZE(A)) * NORM(A,1) * EPS
%	is used for deciding when elements of x are less than zero.
%	This can be overridden with x = fnnlsb(A,b,P_old,Z_old,TOL).
%
%	[x,w] = fnnlsb(A,b) also returns dual vector w where
%	w(i) < 0 where x(i) = 0 and w(i) = 0 where x(i) > 0.
%       [X,W,P,Z] = fnnlsb(A,b) also returns the index set P and Z
%	that explicitly tells which constraints (variables) are
% 	active, which are not. These can be used as subsequent
%	input to fnnlsb if only minor changes are expected in a new
%	problem. P designates the passive (unconstrained) set, while
%	designates the active constraints. For initialisation use
% 	P=zeros(size(Xty));
%	Z=[1:length(Xty)]';
%
%	See also FNNLS and NNLS

%	L. Shure 5-8-87
%	Revised, 12-15-88,8-31-89 LS.
%	(Partly) Copyright (c) 1984-94 by The MathWorks, Inc.

%	Modified by R. Bro 5-7-96 according to
%       Bro R., de Jong S., Journal of Chemometrics, 1997, 11, 393-401
% 	Corresponds to the FNNLSb algorithm in the paper
%
%
%	Rasmus bro
%	Chemometrics Group, Food Technology
%	Dept. Dairy and Food Science
%	Royal Vet. & Agricultural
%	DK-1958 Frederiksberg C
%	Denmark
%	rb@kvl.dk
%	http://newton.foodsci.kvl.dk/rasmus.html

%  Reference:
%  Lawson and Hanson, "Solving Least Squares Problems", Prentice-Hall, 1974.

% initialize variables
if nargin < 5
    tol = 10*eps*norm(XtX,1)*max(size(XtX));
end
[m,n] = size(XtX);
P = P_old(:)';
PP = find(P);
Z = Z_old(:)';
ZZ=find(Z);
iter = 0;
itmax = 30*n;
z=zeros(n,1);
x=z;
z(PP)=(Xty(PP)'/XtX(PP,PP)');

    while any((z(PP) <= tol)) & iter < itmax
        iter = iter + 1;
        QQ = find((z <= tol) & P');
        alpha = min(x(QQ)./(x(QQ) - z(QQ)));
        x = x + alpha*(z - x);
        ij = find(abs(x) < tol & P' ~= 0);
        Z(ij)=ij';
        P(ij)=zeros(1,max(size(ij)));
        PP = find(P);
        ZZ = find(Z);
        nzz = size(ZZ);
        z(PP)=(Xty(PP)'/XtX(PP,PP)');
        z(ZZ) = zeros(nzz(2),nzz(1));
        z=z(:);
    end
    x = z;
    w = Xty-XtX*x;

% set up iteration criterion
iter = 0;

% outer loop to put variables into set to hold positive coefficients
while any(Z) & any(w(ZZ) > tol)
    [wt,t] = max(w(ZZ));
    t = ZZ(t);
    P(1,t) = t;
    Z(t) = 0;
    PP = find(P);
    ZZ = find(Z);
    nzz = size(ZZ);
    z(PP)=(Xty(PP)'/XtX(PP,PP)');
    z(ZZ) = zeros(nzz(2),nzz(1))';
    z=z(:);
    while any((z(PP) <= tol)) & iter < itmax
        iter = iter + 1;
        QQ = find((z <= tol) & P');
        alpha = min(x(QQ)./(x(QQ) - z(QQ)));
        x = x + alpha*(z - x);
        ij = find(abs(x) < tol & P' ~= 0);
        Z(ij)=ij';
        P(ij)=zeros(1,max(size(ij)));
        PP = find(P);
        ZZ = find(Z);
        nzz = size(ZZ);
        z(PP)=(Xty(PP)'/XtX(PP,PP)');
        z(ZZ) = zeros(nzz(2),nzz(1));
        z=z(:);
    end
    x = z;
    w = Xty-XtX*x;
end

x=x(:);


text:
function checkNNGradients(lambda)
%CHECKNNGRADIENTS Creates a small neural network to check the
%backpropagation gradients
%   CHECKNNGRADIENTS(lambda) Creates a small neural network to check the
%   backpropagation gradients, it will output the analytical gradients
%   produced by your backprop code and the numerical gradients (computed
%   using computeNumericalGradient). These two gradient computations should
%   result in very similar values.
%

if ~exist('lambda', 'var') || isempty(lambda)
    lambda = 0;
end

input_layer_size = 3;
hidden_layer_size = 5;
num_labels = 3;
m = 5;

% We generate some 'random' test data
Theta1 = debugInitializeWeights(hidden_layer_size, input_layer_size);
Theta2 = debugInitializeWeights(num_labels, hidden_layer_size);
% Reusing debugInitializeWeights to generate X
X  = debugInitializeWeights(m, input_layer_size - 1);
y  = 1 + mod(1:m, num_labels)';

% Unroll parameters
nn_params = [Theta1(:) ; Theta2(:)];

% Short hand for cost function
costFunc = @(p) nnCostFunction(X,y,p,input_layer_size, hidden_layer_size, ...
                               num_labels,lambda);

[cost, grad] = costFunc(nn_params);
numgrad = computeNumericalGradient(costFunc, nn_params);

% Visually examine the two gradient computations.  The two columns
% you get should be very similar.
disp([numgrad grad]);
fprintf(['The above two columns you get should be very similar.\n' ...
         '(Left-Your Numerical Gradient, Right-Analytical Gradient)\n\n']);

% Evaluate the norm of the difference between two solutions.
% If you have a correct implementation, and assuming you used EPSILON = 0.0001
% in computeNumericalGradient.m, then diff below should be less than 1e-9
diff = norm(numgrad-grad)/norm(numgrad+grad);

fprintf(['If your backpropagation implementation is correct, then \n' ...
         'the relative difference will be small (less than 1e-9). \n' ...
         '\nRelative Difference: %g\n'], diff);

end

author: ShiMengjie

repo: Machine-Learning-Andrew-Ng


text:
function Bc = brownian_bridge(P,K,a,b,sigma)

t = linspace(0,1,P)';
% generate a bunch of bridges
B = (randn(P-1,K)+1i*randn(P-1,K))/sqrt(P);
B = [zeros(1,K); cumsum(B)];
B = B - t * B(end,:);
B = B*sigma;
Bc = B + repmat(1-t,[1 K])*a + repmat(t,[1 K])*b;

end

author: optimaltransport


text:
function [ftlb] = Btu2ftlb(Btu)
% Convert energy or work from British thermal units to foot-pounds.
% Chad A. Greene 2012
ftlb = Btu*778.16932495 ;


text:
function [QDu,area,center,QDlambda,stbElem] = graduVEM(node,elem,u)
%% GRADUVEM the projected gradient of a linear virtual element function.
%
% QDu = GRADUVEM(node,elem,u) compute the L2 projection of the gradient of
% a virtual element function u on a polygonal mesh representing by (node,elem).
%
% [QDu,area,center,QDlambda,stbElem] = GRADUVEM(node,elem,u) also outputs areas, centroids,
% elementwise stabilization, and the L2 projection of Dlambda which is the gradient of P1
% conforming VEM basis. QDu{t}(i,1) is the x-component of the i-th vertex's basis in t-th element.
% stbElem(t,:) is the \ell^2 norm of (I-P)u_h on each element.
%
% Remark: this routine is mostly following Long's vectorization in PoissonVEM,
% however is ONLY implemented for Matlab version > R2016b
% see: https://scaomath.github.io/blog/MATLAB-update-on-array-compatibility/
%
% See also gradu, gradbasis, PoissonVEM
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.
%%
if ~iscell(elem); elem = num2cell(elem,2); end

%%
NT = size(elem,1);
QDu = zeros(NT,2);
area = zeros(NT,1);
center = zeros(NT,2);
QDlambda = cell(NT,1);
stbElem = zeros(NT,1); % || (I- P) u ||_{l^2,K}
%%
elemVertexNumber = cellfun('length',elem);% the number of vertices per element
minNv = min(elemVertexNumber);
maxNv = max(elemVertexNumber);

%%
for nV = minNv:maxNv
    isNv = (elemVertexNumber == nV); % index of elements with Nv vertices
    if ~any(isNv); continue; end
    elemNv = cell2mat(elem(isNv));
    NelemNv = sum(isNv);
    x1 = reshape(node(elemNv,1),[NelemNv,nV]);
    y1 = reshape(node(elemNv,2),[NelemNv,nV]);
    x2 = circshift(x1,[0,-1]);
    y2 = circshift(y1,[0,-1]);
    bdIntegral = x1.*y2 - y1.*x2;
    areaNv = sum(bdIntegral,2)/2;
    xc = sum((x1+x2).*bdIntegral,2)./abs(areaNv)/6;
    yc = sum((y1+y2).*bdIntegral,2)./abs(areaNv)/6;
    center(isNv,:) = [xc, yc];
    normVecx = y2 - y1; % normal vector is a rotation of edge vectors
    normVecy = x1 - x2;
    Bx = (normVecx + circshift(normVecx,[0,1]))./(2*areaNv); % average of normal vectors
    By = (normVecy + circshift(normVecy,[0,1]))./(2*areaNv); % in adjacent edges
    QDlambdaNv = reshape([Bx, By]',[nV,2,NelemNv]);
    QDlambda(isNv) = squeeze(num2cell(QDlambdaNv, [1, 2]));
    uh2elemNv = u(elemNv);
    if NelemNv == 1; uh2elemNv = uh2elemNv'; end % in this case: size(uh2elemNv,2) == 1
    QDudx = sum(uh2elemNv.*Bx, 2);
    QDudy = sum(uh2elemNv.*By, 2);
    QDlambda(isNv) = squeeze(num2cell(QDlambdaNv, [1, 2]));
    %%
    area(isNv,:) = abs(areaNv);
    QDu(isNv,:) = [QDudx, QDudy];

    %% compute stablization
    if nargout > 4
        h = sign(areaNv).*sqrt(abs(areaNv)); % h = sqrt(area) not the diameter
        cx = mean(reshape(node(elemNv(:),1),[NelemNv,nV]),2);
        cy = mean(reshape(node(elemNv(:),2),[NelemNv,nV]),2);
        Dx = (x1 - cx)./h; Dy = (y1 - cy)./h; %  m = (x - cx)/h
        IminusP = zeros(NelemNv,nV,nV);
        for i = 1:nV
            for j = 1:nV
                IminusP(:,i,j) = - 1/nV - Dx(:,i).*Bx(:,j).*abs(h) - Dy(:,i).*By(:,j).*abs(h);
            end
            IminusP(:,i,i) = ones(NelemNv,1) + IminusP(:,i,i);
        end
        if NelemNv == 1
            IminusP = squeeze(IminusP);
            utIminusPu = uh2elemNv*IminusP*uh2elemNv';
        else
            IminusPu = sum(bsxfun(@times, IminusP, uh2elemNv), 2);
            utIminusPu = sum(uh2elemNv.*squeeze(IminusPu),2);
        end

        stbElem(isNv) = sqrt(abs(utIminusPu));
    end
end

repo: ifem


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Jellyfish Example Courtesy of Alexander P. Hoover, PhD
%
% Converted from IBAMR: 1/16/2018 by NAB.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Make_Jelly_Geometry()

close all;
clear all;
L = 8;                              % height of computational domain (m) for keeping desired resolution
Lh = 10;                            % actual height of computational domain (m) (MATCHES INPUT2D)
Lw = 3;                             % width of computational domain (m) (MATCHES INPUT2D)
N = 96;                            % number of Cartesian grid meshwidths at the finest level of the AMR grid
dx = L/N;                           % Cartesian mesh width (m)
ds = dx/2;

a=.5;                               % bell radius (semi-minor axis, horizontal axis, note width=2a)
b=.75;                              % bell semi-major axis
d=-0.25;
factor_a=.8;

F=1e5; %5e0

theta=zeros(1000,1);
theta_lim=asin(d/b);
theta_test=pi/2;

x_points=zeros(1000,1);
z_points=zeros(1000,1);
id_points=zeros(1000,1);
offset = 0;

kappa_spring = 1e7; %1e5               % spring constant (Newton)
kappa_beam = 2.5e5; %1e5    %5e3              % beam stiffness constant (Newton m^2)
%kappa_beam_flexible = kappa_beam/5;   % beam stiffness constant (Newton m^2)
kappa_target = kappa_spring;           % target point penalty spring constant (Newton)

c=0;
while(theta_test<(pi-theta_lim))
    c=c+1;
    theta(c)=theta_test;

    x_points(c)=a*cos(theta(c));
    z_points(c)=b*sin(theta(c));
    id_points(c)=c-1;

    theta_test=ds/((a*sin(theta(c)))^(2)+(b*cos(theta(c)))^(2))^(.5)+theta(c);

end

c_stiff=c;

npts=2*c-1;
npts_wing=floor(npts/2);
npts_musc=floor(npts_wing/4);

for j=(c+1):(npts)
    x_points(j)=-1*x_points(j-c+1);
    z_points(j)=z_points(j-c+1);
    id_points(j)=j-1;
end

mesh_name = 'jelly';
xShift = 1.5;
yShift = 2;

x_points=x_points(1:npts)+xShift;
z_points=z_points(1:npts)+yShift;
it_points=id_points(1:npts);

plot(x_points(:),z_points(:),'*'); hold on;
axis([0 8 0 8])

% Lag Pts to Mess up Flow At Edge
xBlock = ds:4*ds:Lw-ds;
yBlock = (Lh-5*ds)*ones(1,length(xBlock))+ds;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Print .vertex information
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vertex_fid = fopen([mesh_name num2str(N) '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', npts + npts_musc*2 + length(xBlock));
    lag_ct = 0;

    %
    % bell
    %
    for j=1:npts
        fprintf(vertex_fid, '%1.16e %1.16e\n', x_points(j), z_points(j));
        lag_ct = lag_ct + 1;
    end

    %
    % muscles
    %
    for s = 1:npts_musc
        fprintf(vertex_fid, '%1.16e %1.16e\n', x_points(npts_wing+1-npts_musc+s), z_points(npts_wing+1-npts_musc+s));
        plot(x_points(npts_wing+1-npts_musc+s),z_points(npts_wing+1-npts_musc+s),'r*'); hold on;
        lag_ct = lag_ct + 1;
    end
    for s = 1:npts_musc
        fprintf(vertex_fid, '%1.16e %1.16e\n', x_points(npts-npts_musc+s), z_points(npts-npts_musc+s));
        plot(x_points(npts-npts_musc+s),z_points(npts-npts_musc+s),'r*'); hold on;
        lag_ct = lag_ct + 1;
    end

    for ii=1:length(xBlock)
        fprintf(vertex_fid, '%1.16e %1.16e\n', xBlock(ii), yBlock(ii));
    end

fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Print .spring information
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

spring_fid = fopen([mesh_name num2str(N) '.spring'], 'w');

    npts_spring_type1=npts-1;

    fprintf(spring_fid, '%d\n', npts-1 + npts_musc);

    fprintf('\nNumber of springs before muscles: %d \n\n',npts-1)

    factor = 1;%ds^2/ds;

    %
    % bell
    %
    for s = 1:c-1
        resting=sqrt((x_points(s)-x_points(s+1))^(2)+(z_points(s)-z_points(s+1))^(2));
        fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n', id_points(s)+1, id_points(s+1)+1, kappa_spring*ds/(ds^2)*factor, resting, 1);
    end
    for s = c+1:npts-1
        resting=sqrt((x_points(s)-x_points(s+1))^(2)+(z_points(s)-z_points(s+1))^(2));
        fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n', id_points(s)+1, id_points(s+1)+1, kappa_spring*ds/(ds^2)*factor, resting, 1);
    end
    resting=sqrt((x_points(1)-x_points(c+1))^(2)+(z_points(1)-z_points(c+1))^(2));
    fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n', id_points(1)+1, id_points(c+1)+1, kappa_spring*ds/(ds^2)*factor, resting, 1);

    %
    % muscles
    %
    for s = 1:npts_musc
        fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n',npts+s-1+1, npts+s+npts_musc-1+1, F, 0, 1);
    end

    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Print .nonInv_beam information
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

beam_fid = fopen([mesh_name num2str(N) '.nonInv_beam'], 'w');

    fprintf(beam_fid, '%d\n', npts-2);

    factor=1;% = (ds^4)/ds;

    for s = 2:c-1
        C1 = x_points(s-1)+x_points(s+1)-2*x_points(s);
        C2 = z_points(s-1)+z_points(s+1)-2*z_points(s);
        fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(s-1)+1, id_points(s)+1, id_points(s+1)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);
    end
    for s = c+2:npts-1
        C1 = x_points(s-1)+x_points(s+1)-2*x_points(s);
        C2 = z_points(s-1)+z_points(s+1)-2*z_points(s);
        fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(s-1)+1, id_points(s)+1, id_points(s+1)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);
    end

    C1 = x_points(c+2)+x_points(1)-2*x_points(c+1);
    C2 = z_points(c+2)+z_points(1)-2*z_points(c+1);
    fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(c+2)+1, id_points(c+1)+1, id_points(1)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);

    C1 = x_points(c+1)+x_points(2)-2*x_points(1);
    C2 = z_points(c+1)+z_points(2)-2*z_points(1);
    fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(c+1)+1, id_points(1)+1, id_points(2)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);

    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PRINT TARGET POINTS!!!
%
% print target points (flow blocker along edge)
k_Target = 2.5e6;
nBefore = lag_ct; % Counts pts in jellyfish for bookkeeping for .target file
struct_name = ['jelly' num2str(N)];
print_Lagrangian_Target_Pts(xBlock,k_Target,struct_name,nBefore)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called 'struct_name'.target
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name,nBefore)

    N = length(xLag);
    Nstart = nBefore+1;
    Nend = nBefore+N;

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = Nstart:Nend
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

author: nickabattista


text:
function n = p17_n ( )

%*****************************************************************************80
%
%% P17_N returns the number of variables for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 4;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function net_pingIIWA(ip)
%% About
% Ping kuka iiwa through the network

%% Syntax
% net_pingIIWA(ip)

%% Arreguments
% ip: is the IP of the kuka controller

% Copyright: Mohammad SAFEEA 16th-April-2018

command=['ping ',ip];
dos(command);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function lfmVisualise(model, visualiseFunction, visualiseModify, varargin)

% LFMVISUALISE Visualise the outputs in a latent force model
% FORMAT
% DESC visualises a latent force model with two latent forces as inputs
% ARG model :  the model to visualise.
% ARG visualiseFunction : the function that draws the visualisation (in
% data space) when the graphs are first drawn.
% ARG visualiseModify : the function that modifies the visualisation as
% you create the latent function.
% ARG arg1, arg2, arg3, ... : various additional arguments to be passed to the
% visualisation commands.
%
%
% COPYRIGHT : Mauricio Alvarez and Neil D. Lawrence, 2008

% MLTOOLS

global visualiseInfo
kernType = model.kernType(1:3);
if ~strcmpi(kernType, 'lfm')
  error('This function is only implemented for "LFM" kernels');
else
  if length(varargin)~=4,
    error('Include the skeleton, the initial position and two posteriors for the outputs')
  end
end
range = 1.1*max(max(abs([varargin{3}{1} varargin{3}{2}])));
figure(1)
set(gcf, 'Position', [19 202 582 527]);
clf
% Create a black panel as background
f1 = linspace(-range, range, 200);
f2 = linspace(-range, range, 200);
ax = axes('position', [0.15 0.15 0.75 0.75]);
f1Lim = [min(f1) max(f1)];
f2Lim = [min(f2) max(f2)];
set(ax, 'xLim', f1Lim);
set(ax, 'yLim', f2Lim);
set(ax, 'fontname', 'arial');
set(ax, 'fontsize', 15);
set(ax, 'TickLength', [0 0]);

plot(varargin{3}{1}, varargin{3}{2}, 'b.', 'markersize', 20, 'LineWidth', 2)
hold on
plot(varargin{4}{1}, varargin{4}{2}, 'r.', 'markersize', 20, 'LineWidth', 2)
for i = 2:length(varargin{3}{1})
  arrow([varargin{3}{1}(i-1) varargin{3}{1}(i)], ...
        [varargin{3}{2}(i-1) varargin{3}{2}(i)], [], 'b', 2);
end
for i = 2:length(varargin{4}{1})
  arrow([varargin{4}{1}(i-1) varargin{4}{1}(i)], ...
        [varargin{4}{2}(i-1) varargin{4}{2}(i)], [], 'r', 2);
end
a = text(varargin{3}{1}(1), varargin{3}{2}(1), ' Motion 18');
set(a, 'horizontalalignment', 'left', 'color', [0 0 1], 'fontsize', 16)
a = text(varargin{4}{1}(1), varargin{4}{2}(1), 'Motion 19 ');
set(a, 'horizontalalignment', 'right', 'color', [1 0 0], 'fontsize', 16)
xlabel('f_1(t)')
ylabel('f_2(t)')
title('Latent Forces')
grid on
hold off
%
visualiseInfo.plotAxes = ax;

if verLessThan('matlab', 'R2014a')
  visualiseInfo.latentHandle = line(varargin{3}{1}(1), varargin{3}{2}(1), ...
                                    'markersize', 20, ...
                                    'color', [0.5 0.5 0.5], 'marker', ...
                                    '.', 'visible', ...
                                    'on', 'erasemode', 'xor');
else
  visualiseInfo.latentHandle = line(varargin{3}{1}(1), varargin{3}{2}(1), ...
                                    'markersize', 20, ...
                                    'color', [0.5 0.5 0.5], 'marker', ...
                                    '.', 'visible', ...
                                    'on');
end
visualiseInfo.latentLine = [];
visualiseInfo.clicked = 0;
visualiseInfo.digitAxes = [];
visualiseInfo.digitIndex = [];

% visualiseInfo.dynamicsSlider = ...
%     uicontrol('Style', 'slider', ...
%               'String', 'Time', ...
%               'sliderStep', [0.01, 0.1], ...
%               'units', 'normalized', ...
%               'position', [0 0.95 1 0.05], ...
%               'callback', 'lfmClassVisualise(''dynamicsSliderChange'')');

% set(visualiseInfo.dynamicsSlider, 'visible', 'off');
set(gcf, 'WindowButtonMotionFcn', 'lfmClassVisualise(''move'')')
set(gcf, 'WindowButtonDownFcn', 'lfmClassVisualise(''click'')')

visualiseInfo.timer.series = 0;
visualiseInfo.timer.stepTime = 0.03;

%figure(2)
%set(gcf,'Position',[418 332 445 395]);
%clf
%subplot(2,1,1);
%visualiseInfo.f1.handle = plot(0,1.5, 'LineWidth', 2);
%title('f_1(t)','FontSize', 15, 'FontName', 'arial')
%set(gca, 'yLim', [-range range]);
%set(gca, 'fontname', 'arial');
%set(gca, 'fontsize', 15);
visualiseInfo.f1.series = 0;
%subplot(2,1,2);
%visualiseInfo.f2.handle = plot(0,1.5,'LineWidth', 2);
%title('f_2(t)','FontSize', 15, 'FontName', 'arial')
%set(gca, 'yLim', [-range range]);
%set(gca, 'fontname', 'arial');
%set(gca, 'fontsize', 15);
visualiseInfo.f2.series = 0;

figure(2)
set(gcf,'Position',[635 201 582 527]);
clf
visualiseInfo.visualiseFunction = str2func(visualiseFunction);
visHandle = visualiseInfo.visualiseFunction(varargin{1:2});
%set(gca, 'xlim', [-8  18], ...
%         'ylim', [-2 15], ...
%         'zlim', [0 35]);
set(gca, 'xlim', [-20 10], ...
         'ylim', [-20 20], ...
         'zlim', [0 32]);
title('Synthetic Motion', 'FontSize', 15);

% Pass the data to visualiseInfo
visualiseInfo.model = model;
visualiseInfo.varargin = varargin;
visualiseInfo.visualiseModify = str2func(visualiseModify);
visualiseInfo.visHandle = visHandle;

% figure(4)
% clf
% set(gcf,'Position',[418 30 445 395]);
% visualiseInfo.fPos = plot(0,1.5,'LineWidth', 2);
% %set(gca, 'xLim', [-2*range 2*range]);
% %set(gca, 'yLim', [-2*range 2*range]);
% xlabel('f_1(t)')
% ylabel('f_2(t)')
% set(gca, 'fontname', 'arial');
% set(gca, 'fontsize', 15);


text:
function bnd = prepare_mesh_manual(cfg, mri)

% PREPARE_MESH_MANUAL is called by PREPARE_MESH and opens a GUI to manually
% select points/polygons in an mri dataset.
%
% It allows:
%   Visualization of 3d data in 3 different projections
%   Adjustment of brightness for every slice
%   Storage of the data points in an external .mat file
%   Retrieval of previously saved data points
%   Slice fast scrolling with keyboard arrows
%   Polygons or points selection/deselection
%
% See also PREPARE_MESH_SEGMENTATION, PREPARE_MESH_HEADSHAPE

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

% FIXME: control slice's cmap referred to abs values
% FIXME: clean structure slicedata
% FIXME: check function assign3dpoints

global obj

mri = ft_checkdata(mri, 'datatype', {'volume', 'segmentation'}, 'hasunit', 'yes');

hasheadshape = isfield(cfg, 'headshape');
hasbnd       = isfield(cfg, 'bnd');  % FIXME why is this in cfg?
hasmri       = nargin>1;

% check the consistency of the input arguments
if hasheadshape && hasbnd
  ft_error('you should not specify cfg.headshape and cfg.bnd simultaneously');
end

% check the consistency of the input arguments
if ~hasmri
  % FIXME give a warning or so?
  mri.anatomy = [];
  mri.transform = eye(4);
else
  % ensure that it is double precision
  mri.anatomy = double(mri.anatomy);
end

% start with an empty boundary
bnd.pnt = [];
bnd.tri = [];

if hasheadshape
  if ~isempty(cfg.headshape)
    % start with the headshape
    [bnd.pos, bnd.tri] = headsurface([], [], 'headshape', cfg.headshape);
  end
elseif hasbnd
  if ~isempty(cfg.bnd)
    % start with the prespecified boundaries
    bnd = cfg.bnd;
  end
end

% creating the GUI
fig = figure;

% initialize some values
slicedata = []; points = bnd.pnt;
axs = floor(size(mri.anatomy,1)/2);

% initialize properties
set(fig, 'KeyPressFcn',@keypress);
set(fig, 'CloseRequestFcn', @cb_close);
setappdata(fig,'data',mri.anatomy);
setappdata(fig,'prop',{1,axs,0,[]});
setappdata(fig,'slicedata',slicedata);
setappdata(fig,'points',points);
setappdata(fig,'box',[]);

% add GUI elements
cb_creategui(gca);
cb_redraw(gca);
waitfor(fig);

% close sequence
global pnt
try
  tmp = pnt;
  clear global pnt
  pnt = tmp;
  clear tmp
  [tri] = projecttri(pnt);
  bnd.pnt = pnt;
  bnd.tri = tri;
catch
  bnd.pnt = [];
  bnd.tri = [];
end

function cb_redraw(hObject, eventdata, handles)
fig  = get(hObject, 'parent');
prop = getappdata(fig,'prop');
data = getappdata(fig,'data');
slicedata = getappdata(fig,'slicedata');
points = getappdata(fig,'points');

% draw image
try
  cla
  % i,j,k axes
  if (prop{1}==1) %jk
    imh=imagesc(squeeze(data(prop{2},:,:)));
    xlabel('k');
    ylabel('j');
    colormap bone
  elseif (prop{1}==2) %ik
    imh=imagesc(squeeze(data(:,prop{2},:)));
    xlabel('k');
    ylabel('i');
    colormap bone
  elseif (prop{1}==3) %ij
    imh=imagesc(squeeze(data(:,:,prop{2})));
    xlabel('j');
    ylabel('i');
    colormap bone
  end
  axis equal; axis tight
  brighten(prop{3});
  title([ 'slice ' num2str(prop{2}) ])
catch
end
if ~ishold,hold on,end
% draw points and lines
for zz = 1:length(slicedata)
  if(slicedata(zz).proj==prop{1})
    if(slicedata(zz).slice==prop{2})
      polygon = slicedata(zz).polygon;
      for kk=1:length(polygon)
        if ~isempty(polygon{kk})
          [xs,ys]=deal(polygon{kk}(:,1),polygon{kk}(:,2));
          plot(xs, ys, 'g.-');
        end
      end
    end
  end
end % end for

% draw other slices points
points = getappdata(fig,'points');
pnts   = round(points);
if ~isempty(pnts)
  % exclude polygons inside the same slice
  indx  = find(points(:,prop{1})==prop{2});
  tmp   = ones(1,size(points,1));
  tmp(indx)  = 0;
  pnts  = points(find(tmp),:);
  % calculate indexes of other proj falling in current slice
  indx2 = find(round(pnts(:,prop{1}))==prop{2});
  rest  = setdiff([1 2 3],prop{1});
  % plot the points
  for jj=1:length(indx2)
    [x,y] = deal(pnts(indx2(jj),rest(1)),pnts(indx2(jj),rest(2)));
    plot(y,x,'g*')
  end
end

% add blue cross markers in case of box selection
box        = getappdata(fig,'box');
point2mark = getappdata(fig,'point2mark');

if ~isempty(point2mark)
  plot(point2mark.x,point2mark.y,'marker','+')
end

if ~isempty(box)
  for kk=1:size(box,1)
    proj   = box(kk,1);
    slice  = box(kk,2);
    rowmin = box(kk,3);
    rowmax = box(kk,4);
    colmin = box(kk,5);
    colmax = box(kk,6);
    aaa = get(findobj(fig, 'color', 'g'));
    for ii=1:length(aaa)
      if ispolygon(aaa(ii))
        L = length(aaa(ii).YData);
        for jj=1:L
          cond = lt(aaa(ii).YData(jj),rowmax).*gt(aaa(ii).YData(jj),rowmin).* ...
            lt(aaa(ii).XData(jj),colmax).*gt(aaa(ii).XData(jj),colmin);
          if cond
            plot(aaa(ii).XData(jj),aaa(ii).YData(jj),'marker','+')
          end
        end
      elseif ispoint(aaa(ii))
        cond = lt(aaa(ii).YData,rowmax).*gt(aaa(ii).YData,rowmin).* ...
          lt(aaa(ii).XData,colmax).*gt(aaa(ii).XData,colmin);
        % the test after cond takes care the box doesnt propagate in 3D
        if (cond && (proj == prop{1}) && (slice == prop{2}))
          plot(aaa(ii).XData,aaa(ii).YData,'marker','+')
        end
      end
    end
  end
end

if get(findobj(fig, 'tag', 'toggle axes'), 'value')
  axis on
else
  axis off
end

function cb_creategui(hObject, eventdata, handles)

fig = get(hObject, 'parent');
% define the position of each GUI element
% constants
CONTROL_WIDTH  = 0.10;
CONTROL_HEIGHT = 0.075;
CONTROL_HOFFSET = 0.8;
CONTROL_VOFFSET = 0.6;

% control buttons
uicontrol('tag', 'view', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'View', 'value', 1, 'callback', @which_task1);
ft_uilayout(fig, 'tag', 'view',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+4*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);

uicontrol('tag', 'ins', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'Ins', 'value', 0, 'callback', @which_task2);
ft_uilayout(fig, 'tag', 'ins',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+2.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'del', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'Del', 'value', 0, 'callback', @which_task3);
ft_uilayout(fig, 'tag', 'del',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+2.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'slice', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'slice', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'slice',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+1.25*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'min', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '<', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'min',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET+1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'max', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '>', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'max',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET+1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'brightness', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'brightness', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'brightness',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET-0.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'plus', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'plus',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'minus', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'minus',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'toggle axes', 'parent', fig, 'units', 'normalized', 'style', 'checkbox', 'string', 'axes', 'value', 0, 'callback', @cb_redraw);
ft_uilayout(fig, 'tag', 'toggle axes',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET-2.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);

uicontrol('tag', 'plane', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'plane', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'plane',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-3.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'one', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'jk', 'value', 0, 'callback', @cb_btnp1);
ft_uilayout(fig, 'tag', 'one',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-3.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);
uicontrol('tag', 'two', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'ik', 'value', 0, 'callback', @cb_btnp2);
ft_uilayout(fig, 'tag', 'two',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-4*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'three', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'ij', 'value', 0, 'callback', @cb_btnp3);
ft_uilayout(fig, 'tag', 'three',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-4*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'mesh', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'mesh', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'mesh',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET-5.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);

uicontrol('tag', 'smoothm', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'smooth', 'value', [], 'callback', @smooth_mesh);
ft_uilayout(fig, 'tag', 'smoothm',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'viewm', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'view', 'value', [], 'callback', @view_mesh);
ft_uilayout(fig, 'tag', 'viewm',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'cancel', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'cancel', 'value', [], 'callback', @cancel_mesh);
ft_uilayout(fig, 'tag', 'cancel',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'ok', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'ok', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'ok',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

function cb_btn(hObject, eventdata, handles)
fig  = get(gca, 'parent');
prop = getappdata(fig,'prop');
slice = prop{2};
br    = prop{3};
datdim = size(getappdata(fig,'data'));
% p1 = [];p2 = [];p3 = [];
p1  = get(findobj('string','jk'),'value');
p2  = get(findobj('string','ik'),'value');
p3  = get(findobj('string','ij'),'value');
set(findobj('string','-'),'Value',prop{3});
set(findobj('string','+'),'Value',prop{3});
beta = get(findobj('string','-'),'Value');

if (p1)     %jk
  setappdata(fig,'prop',{1,round(datdim(1)/2),br,prop{4}});
  cb_redraw(gca);
elseif (p2) %ik
  setappdata(fig,'prop',{2,round(datdim(2)/2),br,prop{4}});
  cb_redraw(gca);
elseif (p3) %ij
  setappdata(fig,'prop',{3,round(datdim(3)/2),br,prop{4}});
  cb_redraw(gca);
end
if strcmp(get(hObject,'string'),'-')
  beta=beta-0.05;
  set(findobj('string','-'),'Value',beta);
  set(findobj('string','+'),'Value',beta);
  setappdata(fig,'prop',{prop{1},prop{2},beta,prop{4}});
  cb_redraw(gca);
elseif strcmp(get(hObject,'string'),'+')
  beta=beta+0.05;
  set(findobj('string','+'),'Value',beta);
  setappdata(fig,'prop',{prop{1},prop{2},beta,prop{4}});
  cb_redraw(gca);
end

if strcmp(get(hObject,'string'),'<')
  setappdata(fig,'prop',{prop{1},slice-1,0,prop{4}});
  cb_redraw(gca);
end
if strcmp(get(hObject,'string'),'>')
  setappdata(fig,'prop',{prop{1},slice+1,0,prop{4}});
  cb_redraw(gca);
end
if strcmp(get(hObject,'string'),'OK')
  close(fig)
end

function keypress(h, eventdata, handles, varargin)
fig   = get(gca, 'parent');
prop  = getappdata(fig,'prop');
dat   = guidata(gcbf);
key   = get(gcbf, 'CurrentCharacter');
slice = prop{2};
if  key
  switch key
    case 28
      setappdata(fig,'prop',{prop{1},slice-1,0,prop{4}});
      cb_redraw(gca);
    case 29
      setappdata(fig,'prop',{prop{1},slice+1,0,prop{4}});
      cb_redraw(gca);
    otherwise
  end
end
uiresume(h)

function build_polygon
fig   = get(gca, 'parent');
prop  = getappdata(fig,'prop');
data  = getappdata(fig,'data');
ss = size(data);
proj  = prop{1};
slice = prop{2};
thispolygon =1;
polygon{thispolygon} = zeros(0,2);
maskhelp = [ ...
  '------------------------------------------------------------------------\n' ...
  'specify polygons for masking the topographic interpolation\n' ...
  'press the right mouse button to add another point to the current polygon\n' ...
  'press backspace on the keyboard to remove the last point\n' ...
  'press "c" on the keyboard to close this polygon and start with another\n' ...
  'press "q" or ESC on the keyboard to continue\n' ...
  ];
again = 1;
if ~ishold,hold,end
fprintf(maskhelp);
fprintf('\n');

while again
  for i=1:length(polygon)
    fprintf('polygon %d has %d points\n', i, size(polygon{i},1));
  end

  [x, y, k] = ginput(1);
  okflag = 0;

  % check points do not fall out of image boundaries
  if get(findobj('string','Ins'),'value')
    if (proj == 1 && y<ss(2) && x<ss(3) && x>1 && y>1)
      okflag = 1;
    elseif (proj == 2 && y<ss(1) && x<ss(3) && x>1 && y>1)
      okflag = 1;
    elseif (proj == 3 && y<ss(1) && x<ss(2) && x>1 && y>1)
      okflag = 1;
    else
      okflag = 0;
    end
  end

  if okflag
    switch lower(k)
      case 1
        polygon{thispolygon} = cat(1, polygon{thispolygon}, [x y]);
        % add the last line segment to the figure
        if size(polygon{thispolygon},1)>1
          x = polygon{i}([end-1 end],1);
          y = polygon{i}([end-1 end],2);
        end
        plot(x, y, 'g.-');

      case 8 % backspace
        if size(polygon{thispolygon},1)>0
          % redraw existing polygons
          cb_redraw(gca);
          % remove the last point of current polygon
          polygon{thispolygon} = polygon{thispolygon}(1:end-1,:);
          for i=1:length(polygon)
            x = polygon{i}(:,1);
            y = polygon{i}(:,2);
            if i~=thispolygon
              % close the polygon in the figure
              x(end) = x(1);
              y(end) = y(1);
            end
            set(gca,'nextplot','new')
            plot(x, y, 'g.-');
          end
        end

      case 'c'
        if size(polygon{thispolygon},1)>0
          % close the polygon
          polygon{thispolygon}(end+1,:) = polygon{thispolygon}(1,:);
          % close the polygon in the figure
          x = polygon{i}([end-1 end],1);
          y = polygon{i}([end-1 end],2);
          plot(x, y, 'g.-');
          % switch to the next polygon
          thispolygon = thispolygon + 1;
          polygon{thispolygon} = zeros(0,2);
          slicedata = getappdata(fig,'slicedata');
          m = length(slicedata);
          slicedata(m+1).proj    = prop{1};
          slicedata(m+1).slice   = prop{2};
          slicedata(m+1).polygon = polygon;
          setappdata(fig,'slicedata',slicedata);
        end

      case {'q', 27}
        if size(polygon{thispolygon},1)>0
          % close the polygon
          polygon{thispolygon}(end+1,:) = polygon{thispolygon}(1,:);
          % close the polygon in the figure
          x = polygon{i}([end-1 end],1);
          y = polygon{i}([end-1 end],2);
          plot(x, y, 'g.-');
        end
        again = 0;
        %set(gcf,'WindowButtonDownFcn','');
        slicedata = getappdata(fig,'slicedata');
        m = length(slicedata);
        if ~isempty(polygon{thispolygon})
          slicedata(m+1).proj    = prop{1};
          slicedata(m+1).slice   = prop{2};
          slicedata(m+1).polygon = polygon;
          setappdata(fig,'slicedata',slicedata);
        end
      otherwise
        ft_warning('invalid button (%d)', k);
    end

  end
  assign_points3d
end

function cb_btn_dwn(hObject, eventdata, handles)
build_polygon
uiresume

function cb_btn_dwn2(hObject, eventdata, handles)
global obj
erase_points2d(obj)
uiresume

function erase_points2d (obj)
fig       = get(gca, 'parent');
prop      = getappdata(fig,'prop');
slicedata = getappdata(fig,'slicedata');
pntsslice = [];

% 2d slice selected box boundaries
if strcmp(get(obj,'string'),'box')
  [x, y] = ft_select_box;
  rowmin = min(y);    rowmax = max(y);
  colmin = min(x);    colmax = max(x);
  box = getappdata(fig,'box');
  box = [box; prop{1} prop{2} rowmin rowmax colmin colmax];
  setappdata(fig,'box',box);
else
  % converts lines to points
  pos = slicedata2pnts(prop{1},prop{2});
  tmp = ft_select_point(pos);
  point2mark.x = tmp(1); point2mark.y = tmp(2);
  setappdata(fig,'point2mark',point2mark);
end

maskhelp = [ ...
  '------------------------------------------------------------------------\n' ...
  'Delete points:\n' ...
  'press "c" on the keyboard to undo all selections\n' ...
  'press "d" or DEL on the keyboard to delete all selections\n' ...
  ];
fprintf(maskhelp);
fprintf('\n');
again = 1;
if ~ishold,hold,end
cb_redraw(gca);

% waits for a key press
while again
  [junk1, junk2, k] = ginput(1);

  switch lower(k)
    case 'c'
      if strcmp(get(obj,'string'),'box')
        % undoes all the selections, but only in the current slice
        ind = [];
        for ii=1:size(box,1)
          if box(ii,1)==prop{1} && box(ii,2)==prop{2}
            ind = [ind,ii];
          end
        end
        box(ind,:) = [];
        setappdata(fig,'box',box);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      else
        setappdata(fig,'point2mark',[]);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      end
    case {'d', 127}
      if strcmp(get(obj,'string'),'box')
        update_slicedata
        % deletes only the points selected in the current slice
        ind = [];
        for ii=1:size(box,1)
          if box(ii,1)==prop{1} && box(ii,2)==prop{2}
            ind = [ind,ii];
          end
        end
        box(ind,:) = [];
        setappdata(fig,'box',box);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      else
        update_slicedata
        setappdata(fig,'point2mark',[]);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      end
    otherwise
  end
end
assign_points3d

function assign_points3d
fig       = get(gca, 'parent');
slicedata = getappdata(fig,'slicedata');

points = [];
for zz = 1:length(slicedata)
  slice   = slicedata(zz).slice;
  proj    = slicedata(zz).proj;
  polygon = slicedata(zz).polygon;
  for kk=1:length(polygon)
    if ~isempty(polygon{kk})
      [xs,ys]=deal(polygon{kk}(:,1),polygon{kk}(:,2));
      if (proj==1)     %jk
        tmp    = [ slice*ones(length(xs),1),ys(:),xs(:) ];
        points = [points; unique(tmp,'rows')];
      elseif (proj==2) %ik
        tmp    = [ ys(:),slice*ones(length(xs),1),xs(:) ];
        points = [points; unique(tmp,'rows')];
      elseif (proj==3) %ij
        tmp    = [ ys(:),xs(:),slice*ones(length(xs),1) ];
        points = [points; unique(tmp,'rows')];
      end
    end
  end
end
setappdata(fig,'points',points);

function update_slicedata
fig       = get(gca, 'parent');
prop      = getappdata(fig,'prop');
slicedata = getappdata(fig,'slicedata');
box       = getappdata(fig,'box');
point2mark = getappdata(fig,'point2mark');
% case of box selection
if ~isempty(box) && ~isempty(slicedata)
  for zz = 1:length(slicedata)
    slice   = slicedata(zz).slice;
    proj    = slicedata(zz).proj;
    polygon = slicedata(zz).polygon;
    for uu=1:size(box,1)
      if (box(uu,1)==proj) && (box(uu,2)==slice)
        for kk=1:length(polygon)
          if ~isempty(polygon{kk})
            [xs,ys] = deal(polygon{kk}(:,1),polygon{kk}(:,2));
            tmpind = lt(ys,ones(length(xs),1)*box(uu,4)).*gt(ys,ones(length(xs),1)*box(uu,3)).* ...
              lt(xs,ones(length(xs),1)*box(uu,6)).*gt(xs,ones(length(xs),1)*box(uu,5));
            slicedata(zz).polygon{kk} = [ xs(find(~tmpind)),ys(find(~tmpind)) ];
          end
        end
      end
    end
  end
  % case of point selection
else
  if ~isempty(slicedata)
    for zz = 1:length(slicedata)
      slice   = slicedata(zz).slice;
      proj    = slicedata(zz).proj;
      polygon = slicedata(zz).polygon;
      for kk=1:length(polygon)
        if ~isempty(polygon{kk})
          [xs,ys] = deal(polygon{kk}(:,1),polygon{kk}(:,2));
          tmpind = eq(xs,point2mark.x).*eq(ys,point2mark.y);
          slicedata(zz).polygon{kk} = [ xs(find(~tmpind)),ys(find(~tmpind)) ];
        end
      end
    end
  end
end
setappdata(fig,'slicedata',slicedata)

% FIXME: obsolete: replaced by headshape at the beginning
function smooth_mesh(hObject, eventdata, handles)
fig = get(gca,'parent');
disp('not yet implemented')

% FIXME: make it work for pre-loaded meshes
function view_mesh(hObject, eventdata, handles)
fig  = get(gca, 'parent');
assign_points3d
pnt  = getappdata(fig,'points');
pnt_ = pnt;
if ~isempty(pnt)
  tri = projecttri(pnt);
  bnd.pnt = pnt_;
  bnd.tri = tri;
  slicedata = getappdata(fig,'slicedata');
  figure
  ft_plot_mesh(bnd,'vertexcolor','k');
end

function cancel_mesh(hObject, eventdata, handles)
fig  = get(gca, 'parent');
close(fig)

function cb_close(hObject, eventdata, handles)
% get the points from the figure
fig    = hObject;
global pnt
pnt   = getappdata(fig, 'points');
set(fig, 'CloseRequestFcn', @delete);
delete(fig);

function cb_btnp1(hObject, eventdata, handles)
set(findobj('string','jk'),'value',1);
set(findobj('string','ik'),'value',0);
set(findobj('string','ij'),'value',0);
cb_btn(hObject);

function cb_btnp2(hObject, eventdata, handles)
set(findobj('string','jk'),'value',0);
set(findobj('string','ik'),'value',1);
set(findobj('string','ij'),'value',0);
cb_btn(hObject);

function cb_btnp3(hObject, eventdata, handles)
set(findobj('string','jk'),'value',0);
set(findobj('string','ik'),'value',0);
set(findobj('string','ij'),'value',1);
cb_btn(hObject);

function which_task1(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');
set(findobj('string','View'),'value',1);
set(findobj('string','Ins'), 'value',0);
set(findobj('string','Del'), 'value',0);
set(gcf,'WindowButtonDownFcn','');
change_panel

function which_task2(hObject, eventdata, handles)
set(findobj('string','View'),'value',0);
set(findobj('string','Ins'), 'value',1);
set(findobj('string','Del'), 'value',0);
set(gcf,'WindowButtonDownFcn',@cb_btn_dwn);
change_panel

function which_task3(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');
set(findobj('string','View'),'value',0);
set(findobj('string','Ins'), 'value',0);
set(findobj('string','Del'), 'value',1);
change_panel

function change_panel
% affect panel visualization
w1 = get(findobj('string','View'),'value');
w2 = get(findobj('string','Ins'), 'value');
w3 = get(findobj('string','Del'), 'value');

if w1
  set(findobj('tag','slice'),'string','slice');
  set(findobj('tag','min'),'string','<');
  set(findobj('tag','max'),'string','>');
  set(findobj('tag','min'),'style','pushbutton');
  set(findobj('tag','max'),'style','pushbutton');
  set(findobj('tag','min'),'callback',@cb_btn);
  set(findobj('tag','max'),'callback',@cb_btn);
  set(findobj('tag','mesh'),'visible','on');
  set(findobj('tag','openm'),'visible','on');
  set(findobj('tag','viewm'),'visible','on');
  set(findobj('tag','brightness'),'visible','on');
  set(findobj('tag','plus'),'visible','on');
  set(findobj('tag','minus'),'visible','on');
  set(findobj('tag','toggle axes'),'visible','on');
  set(findobj('tag','plane'),'visible','on');
  set(findobj('tag','one'),'visible','on');
  set(findobj('tag','two'),'visible','on');
  set(findobj('tag','three'),'visible','on');
elseif w2
  set(findobj('tag','slice'),'string','select points');
  set(findobj('tag','min'),'string','close');
  set(findobj('tag','max'),'string','next');
  set(findobj('tag','min'),'style','pushbutton');
  set(findobj('tag','max'),'style','pushbutton');
  set(findobj('tag','min'),'callback',@tins_close);
  set(findobj('tag','max'),'callback',@tins_next);
  set(findobj('tag','mesh'),'visible','off');
  set(findobj('tag','openm'),'visible','off');
  set(findobj('tag','viewm'),'visible','off');
  set(findobj('tag','brightness'),'visible','off');
  set(findobj('tag','plus'),'visible','off');
  set(findobj('tag','minus'),'visible','off');
  set(findobj('tag','toggle axes'),'visible','off');
  set(findobj('tag','plane'),'visible','off');
  set(findobj('tag','one'),'visible','off');
  set(findobj('tag','two'),'visible','off');
  set(findobj('tag','three'),'visible','off');
elseif w3
  set(findobj('tag','slice'),'string','select points');
  set(findobj('tag','min'),'string','box');
  set(findobj('tag','max'),'string','point');
  set(findobj('tag','min'),'style','pushbutton');
  set(findobj('tag','max'),'style','pushbutton');
  set(findobj('tag','min'),'callback',@tdel_box);
  set(findobj('tag','max'),'callback',@tdel_single);
  set(findobj('tag','mesh'),'visible','off');
  set(findobj('tag','openm'),'visible','off');
  set(findobj('tag','viewm'),'visible','off');
  set(findobj('tag','brightness'),'visible','off');
  set(findobj('tag','plus'),'visible','off');
  set(findobj('tag','minus'),'visible','off');
  set(findobj('tag','toggle axes'),'visible','off');
  set(findobj('tag','plane'),'visible','off');
  set(findobj('tag','one'),'visible','off');
  set(findobj('tag','two'),'visible','off');
  set(findobj('tag','three'),'visible','off');
end

function tins_close(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');

function tins_next(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');

function tdel_box(hObject, eventdata, handles)
global obj
obj = gco;
set(gcf,'WindowButtonDownFcn',@cb_btn_dwn2);

function tdel_single(hObject, eventdata, handles)
global obj
obj = gco;
set(gcf,'WindowButtonDownFcn',@cb_btn_dwn2);

% accessory functions:
function [pnts] = slicedata2pnts(proj,slice)
fig       = get(gca, 'parent');
slicedata = getappdata(fig,'slicedata');
pnts = [];
for i=1:length(slicedata)
  if slicedata(i).slice==slice && slicedata(i).proj == proj
    for j=1:length(slicedata(i).polygon)
      if ~isempty(slicedata(i).polygon{j})
        tmp = slicedata(i).polygon{j};
        xs = tmp(:,1); ys = tmp(:,2);
        pnts = [pnts;[unique([xs,ys],'rows')]];
      end
    end
  end
end

function h = ispolygon(x)
h = length(x.XData)>1;

function h = ispoint(x)
h = length(x.XData)==1;

function [T,P] = points2param(pnt)
x = pnt(:,1); y = pnt(:,2); z = pnt(:,3);
[phi,theta,R] = cart2sph(x,y,z);
theta = theta + pi/2;
phi   = phi + pi;
T     = theta(:);
P     = phi(:);

function [bnd2,a,Or] = spherical_harmonic_mesh(bnd, nl)

% SPHERICAL_HARMONIC_MESH realizes the smoothed version of a mesh contained in the first
% argument bnd. The boundary argument (bnd) contains typically 2 fields
% called .pnt and .tri referring to vertices and triangulation of a mesh.
% The degree of smoothing is given by the order in the second input: nl
%
% Use as
%   bnd2 = spherical_harmonic_mesh(bnd, nl);

% nl    = 12; % from van 't Ent paper

bnd2 = [];

% calculate midpoint
Or = mean(bnd.pnt);
% rescale all points
x = bnd.pnt(:,1) - Or(1);
y = bnd.pnt(:,2) - Or(2);
z = bnd.pnt(:,3) - Or(3);
X = [x(:), y(:), z(:)];

% convert points to parameters
[T,P] = points2param(X);

% basis function
B     = shlib_B(nl, T, P);
Y     = B'*X;

% solve the linear system
a     = pinv(B'*B)*Y;

% build the surface
Xs = zeros(size(X));
for l = 0:nl-1
  for m = -l:l
    if m<0
      Yml = shlib_Yml(l, abs(m), T, P);
      Yml = (-1)^m*conj(Yml);
    else
      Yml = shlib_Yml(l, m, T, P);
    end
    indx = l^2 + l + m+1;
    Xs = Xs + Yml*a(indx,:);
  end
  fprintf('%d of %d\n', l, nl);
end
% Just take the real part
Xs = real(Xs);
% reconstruct the matrices for plotting.
xs = reshape(Xs(:,1), size(x,1), size(x,2));
ys = reshape(Xs(:,2), size(x,1), size(x,2));
zs = reshape(Xs(:,3), size(x,1), size(x,2));

bnd2.pnt = [xs(:)+Or(1) ys(:)+Or(2) zs(:)+Or(3)];
[bnd2.tri] = projecttri(bnd.pnt);

function B = shlib_B(nl, theta, phi)
% function B = shlib_B(nl, theta, phi)
%
% Constructs the matrix of basis functions
% where b(i,l^2 + l + m) = Yml(theta, phi)
%
% See also: shlib_Yml.m, shlib_decomp.m, shlib_gen_shfnc.m
%
% Dr. A. I. Hanna (2006).
B = zeros(length(theta), nl^2+2*nl+1);
for l = 0:nl-1
  for m = -l:l
    if m<0
      Yml = shlib_Yml(l, abs(m), theta, phi);
      Yml = (-1)^m*conj(Yml);
    else
      Yml = shlib_Yml(l, m, theta, phi);
    end
    indx = l^2 + l + m+1;
    B(:, indx) = Yml;
  end
end
return;

function Yml = shlib_Yml(l, m, theta, phi)
% function Yml = shlib_Yml(l, m, theta, phi)
%
% MATLAB function that takes a given order and degree, and the matrix of
% theta and phi and constructs a spherical harmonic from these. The
% analogue in the 1D case would be to give a particular frequency.
%
% Inputs:
%  m - order of the spherical harmonic
%  l - degree of the spherical harmonic
%  theta - matrix of polar coordinates \theta \in [0, \pi]
%  phi - matrix of azimuthal cooridinates \phi \in [0, 2\pi)
%
% Example:
%
% [x, y] = meshgrid(-1:.1:1);
% z = x.^2 + y.^2;
% [phi,theta,R] = cart2sph(x,y,z);
% Yml = spharm_aih(2,2, theta(:), phi(:));
%
% See also: shlib_B.m, shlib_decomp.m, shlib_gen_shfnc.m
%
% Dr. A. I. Hanna (2006)
Pml=legendre(l,cos(theta));
if l~=0
  Pml=squeeze(Pml(m+1,:,:));
end
Pml = Pml(:);
% Yml = sqrt(((2*l+1)/4).*(factorial(l-m)/factorial(l+m))).*Pml.*exp(sqrt(-1).*m.*phi);
% new:
Yml = sqrt(((2*l+1)/(4*pi)).*(factorial(l-m)/factorial(l+m))).*Pml.*exp(sqrt(-1).*m.*phi);
return;

author: spm


text:
function y = stepLTEMIMOChan(in, prmLTEPDSCH)
%#codegen

% Get simulation params
numTx = prmLTEPDSCH.numTx;
numRx = prmLTEPDSCH.numRx;
chanSRate = prmLTEPDSCH.chanSRate;
SampleTime = 1/chanSRate;

EPAPathDelays = [0 30 70 90 110 190 410]*1e-9;
EPAPathGains  = [0 -1 -2 -3 -8 -17.2 -20.8];
EVAPathDelays = [0 30 150 310 370 710 1090 1730 2510]*1e-9;
EVAPathGains  = [0 -1.5 -1.4 -3.6 -0.6 -9.1 -7 -12 -16.9];
ETUPathDelays = [0 50 120 200 230 500 1600 2300 5000]*1e-9;
ETUPathGains  = [-1 -1 -1 0 0 0 -3 -5 -7];

switch prmLTEPDSCH.chanMdl
    case 'Static MIMO'
        PathDelays          = 0;
        AveragePathGains    = 0;
        MaximumDopplerShift = 0;
    case 'EPA 5Hz'
        PathDelays          = EPAPathDelays;
        AveragePathGains    = EPAPathGains;
        MaximumDopplerShift = 5;
    case 'EVA 5Hz'
        PathDelays          = EVAPathDelays;
        AveragePathGains    = EVAPathGains;
        MaximumDopplerShift = 5;
    case 'EVA 70Hz'
        PathDelays          = EVAPathDelays;
        AveragePathGains    = EVAPathGains;
        MaximumDopplerShift = 70;
    case 'ETU 70Hz'
        PathDelays          = ETUPathDelays;
        AveragePathGains    = ETUPathGains;
        MaximumDopplerShift = 70;
    case 'ETU 300Hz'
        PathDelays          = ETUPathDelays;
        AveragePathGains    = ETUPathGains;
        MaximumDopplerShift = 300;
end

% Create the MIMO Channel object
persistent chanObj;
if isempty(chanObj)
    chanObj = mimochan(numTx, numRx, SampleTime, ...
                       MaximumDopplerShift, PathDelays, AveragePathGains);
end

% Process data through MIMO Channel
y  = filter(chanObj, in);


text:
function [wts,cfreqs] = gammatone_matrix(nfft, sr, nfilts, width, minfreq, maxfreq, maxlen)
% [wts,cfreqa] = fft2gammatonemx(nfft, sr, nfilts, width, minfreq, maxfreq, maxlen)
%      Generate a matrix of weights to combine FFT bins into
%      Gammatone bins.  nfft defines the source FFT size at
%      sampling rate sr.  Optional nfilts specifies the number of
%      output bands required (default 64), and width is the
%      constant width of each band in Bark (default 1).
%      minfreq, maxfreq specify range covered in Hz (100, sr/2).
%      While wts has nfft columns, the second half are all zero.
%      Hence, aud spectrum is
%      fft2gammatonemx(nfft,sr)*abs(fft(xincols,nfft));
%      maxlen truncates the rows to this many bins.
%      cfreqs returns the actual center frequencies of each
%      gammatone band in Hz.
%
% 2004-09-05  Dan Ellis dpwe@ee.columbia.edu  based on rastamat/audspec.m
% Last updated: $Date: 2009/02/22 02:29:25 $

if nargin < 2;    sr = 16000; end
if nargin < 3;    nfilts = 64; end
if nargin < 4;    width = 0.5; end
if nargin < 5;    minfreq = 50; end
if nargin < 6;    maxfreq = sr/2; end
if nargin < 7;    maxlen = nfft/2; end

wts = zeros(nfilts, nfft);

% after Slaney's MakeERBFilters
EarQ = 9.26449;
minBW = 24.7;
order = 1;

cfreqs = -(EarQ*minBW) + exp((1:nfilts)'*(-log(maxfreq + EarQ*minBW) + ...
                log(minfreq + EarQ*minBW))/nfilts) * (maxfreq + EarQ*minBW);
cfreqs = flipud(cfreqs);

GTord = 4;

ucirc = exp(j*2*pi*[0:(nfft/2)]/nfft);

justpoles = 0;

for i = 1:nfilts
  cf = cfreqs(i);
  ERB = width*((cf/EarQ).^order + minBW^order).^(1/order);
  B = 1.019*2*pi*ERB;
  r = exp(-B/sr);
  theta = 2*pi*cf/sr;
  pole = r*exp(j*theta);

  if justpoles == 1
    % point on unit circle of maximum gain, from differentiating magnitude
    cosomegamax = (1+r*r)/(2*r)*cos(theta);
    if abs(cosomegamax) > 1
      if theta < pi/2;  omegamax = 0;
      else              omegamax = pi;   end
    else
      omegamax = acos(cosomegamax);
    end
    center = exp(j*omegamax);
    gain = abs((pole-center).*(pole'-center)).^GTord;
    wts(i,1:(nfft/2+1)) = gain * (abs((pole-ucirc).*(pole'- ...
                                                     ucirc)).^-GTord);
  else
    % poles and zeros, following Malcolm's MakeERBFilter
    T = 1/sr;
    A11 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3+2^1.5)*T*sin(2* ...
                                                      cf*pi*T)./exp(B*T))/2;
    A12 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3+2^1.5)*T*sin(2* ...
                                                      cf*pi*T)./exp(B*T))/2;
    A13 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3-2^1.5)*T*sin(2* ...
                                                      cf*pi*T)./exp(B*T))/2;
    A14 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3-2^1.5)*T*sin(2* ...
                                                      cf*pi*T)./exp(B*T))/2;
    zros = -[A11 A12 A13 A14]/T;

    gain(i) =  abs((-2*exp(4*j*cf*pi*T)*T + ...
                2*exp(-(B*T) + 2*j*cf*pi*T).*T.* ...
                (cos(2*cf*pi*T) - sqrt(3 - 2^(3/2))* ...
                 sin(2*cf*pi*T))) .* ...
               (-2*exp(4*j*cf*pi*T)*T + ...
                2*exp(-(B*T) + 2*j*cf*pi*T).*T.* ...
                (cos(2*cf*pi*T) + sqrt(3 - 2^(3/2)) * ...
                 sin(2*cf*pi*T))).* ...
               (-2*exp(4*j*cf*pi*T)*T + ...
                2*exp(-(B*T) + 2*j*cf*pi*T).*T.* ...
                (cos(2*cf*pi*T) - ...
                 sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) .* ...
               (-2*exp(4*j*cf*pi*T)*T + 2*exp(-(B*T) + 2*j*cf*pi*T).*T.* ...
                (cos(2*cf*pi*T) + sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) ./ ...
               (-2 ./ exp(2*B*T) - 2*exp(4*j*cf*pi*T) +  ...
                2*(1 + exp(4*j*cf*pi*T))./exp(B*T)).^4);
    wts(i,1:(nfft/2+1)) = ((T^4)/gain(i)) ...
        * abs(ucirc-zros(1)).*abs(ucirc-zros(2))...
        .*abs(ucirc-zros(3)).*abs(ucirc-zros(4))...
        .*(abs((pole-ucirc).*(pole'-ucirc)).^-GTord);
  end
end

wts = wts(:,1:maxlen);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function movePTPEllipse_1(iiwa,c,dir,ratio,theta,velocity,accel,TefTool)
%% This funciton is used to
% Move end-effector of the robot on an ellipse.
% this function is used from within the KST class

%% Syntax:
% movePTPEllipse_1(iiwa,c,dir,ratio,TefTool)

%% Arreguments:
% iiwa: KST obejct.
% c: 3x1 vector, is a vector defining the displacment of the the center
% point of the ellipse in relation to the start point of the ellipse,
% position units in (mm).
% dir: the direction vector of the (a) axis of the ellipse.
% ratio: the radious ratio (a/b) of the ellipse.
% theta: is the angle of the ellipe part , for drawing a complete
% ellipse this angle is equal to 2*pi.
% TefTool: 4X4 Transform matrix of the end-efector in the reference frame
% of the flange of the KUKA iiwa, position units in (mm).
% accel: is the acceleration (mm/sec2).
% velocity: is the motion velocity (mm/sec).

% Copy right: Mohammad SAFEEA
% 25th-June-2018

%% Convert inputs to a column vectors
c=colVec(c);
dir=colVec(dir);
ratio=colVec(ratio);
theta=colVec(theta);
velocity=colVec(velocity);
accel=colVec(accel);

%% Check input variables
if(size(c,1)~=3)
fprintf('Error: the vector (c) shall be a 3x1 vector \n');
return;
end

if(size(dir,1)~=3)
fprintf('Error: the vector (dir) shall be a 3x1 vector \n');
return;
end

if(size(ratio,1)~=1)
fprintf('Error: the ratio shall be a scalar  \n');
return;
end

if(size(theta,1)~=1)
fprintf('Error: the variable (theta) shall be a scalar  \n');
return;
end

if(size(TefTool,1)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(TefTool,2)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(velocity,1)~=1)
fprintf('Error: the velocity on the path shall be a scalar  \n');
return;
end

if(size(accel,1)~=1)
fprintf('Error: the acceleration on the path shall be a scalar  \n');
return;
end

%% Convert the units to (meter)
unitConverter=1000;
c=c/unitConverter;
velocity=velocity/unitConverter;
accel=accel/unitConverter;
TefTool(1:3,4)=TefTool(1:3,4)/unitConverter;

%% start the direct servo
     iiwa.realTime_startDirectServoJoints();

%% get current joints angles of robot
    jPos  = iiwa.getJointsPos();

%% Calculate the current position of TCP point of the robot (p0)
    qs=zeros(7,1);
for i=1:7
    qs(i)=jPos{i};
end
        T0=iiwa.directKinematics(qs); % TCP point transformation matrix
        p0=T0(1:3,4); % Current position of TCP point of the robot
        Tt=T0;

%% Calculate the parameters of the ellipse,
p=p0; % this is the starting point of the ellipse
c=c+p; % calculate the center coordinates from the center displacment vector
[R,theta0,a,b,c,errorFlag]=getEllipseParameters(p,c,ratio,dir);
if errorFlag==true % this happens when the start point is on the (dir) vector direction.
    % in such case the ellipse plane can not be specified in space.
    iiwa.realTime_stopDirectServoJoints();
    beep;
    fprintf('Error executing the ellipse function:\n');
    fprintf('Starting point of the ellipse shall not be on the line aligned with the (dir) vector and passing from the center (c)\n');
    return;
end

%% Get the length of the ellipse arc as a function of the parametric angle (theta)
theta1=theta0+theta;
[ thetaVec,sVec ] = getEllipseLengthVector( a,b, theta0,theta1 );
sizesVec=max(size(sVec));
L=sVec(end); % L is the length of the elliptic curve

%% Calculate the times:
[t0,t1,t2]=calculateInterpolationTimes(L,velocity,accel);

%% Joint space control

    [Ttemp,J]=iiwa.directKinematics(qs);
    vec=Ttemp(1:3,4);

%% dls solver parameters
 numberOfIterationForSolver=100;
 lambda=0.5;
 TefTool=eye(4);

    sCoordinate=0;
    interpolationCounter=1;
    thetaVar=theta0;

%% start of time
        dateVector0=datevec(now);
        time0=dateVector0(6)+dateVector0(5)*60+dateVector0(4)*60*60; % calculate time at this instant
while true

    % Calculate the elapsed time
        dateVector=datevec(now);
        timeNow=dateVector(6)+dateVector(5)*60+dateVector(4)*60*60; % calculate time at this instant
        deltaT=timeNow-time0; % elapsed is zero at first excution

    % calculate position of servo point
    if deltaT<t0
        sCoordinate=0.5*accel*deltaT*deltaT;
    elseif(deltaT<t1)
        s0=0.5*accel*t0*t0;
        sCoordinate=velocity*(deltaT-t0)+s0;
    elseif(deltaT<t2)
           s0=0.5*accel*t0*t0;
           s1=velocity*(t1-t0)+s0;
           sCoordinate=velocity*(deltaT-t1)-0.5*accel*(deltaT-t1)*(deltaT-t1)+s1;
    end
 %% when time ends break the loop
    if deltaT>t2
        break;
    end
%% Interpolate theta from the sVec,thetaVec
    for counter=interpolationCounter:(sizesVec-1)
        if(sCoordinate>sVec(counter))
            tetaRange=thetaVec(counter+1)-thetaVec(counter);
            sRange=sVec(counter+1)-sVec(counter);
            thetaVar=(sCoordinate-sVec(counter))...
                *tetaRange/sRange+thetaVec(counter);
            interpolationCounter=counter;
        end
    end

    pPrime=R*ellipseParametricFunction( a,b,thetaVar );
    p=c+pPrime;

    % calculate target transform
    Tt(1:3,4)=p;
    [ qs ] = iiwa.gen_InverseKinematics( qs, Tt,numberOfIterationForSolver,lambda );

    for i=1:7
        jPos{i}=qs(i);
    end

	%% Send joint positions to robot
	iiwa.sendJointsPositions( jPos);

end

    % Stop the reltime motion
    iiwa.realTime_stopDirectServoJoints();

end

function [R,theta0,a,b,c,errorFlag]=getEllipseParameters(p,c,ratio,dir)

%% About
% This function returns the parameters of the ellipse

%% Arreguments:
% p: the starting point of the ellipse
% c: the center of the ellipse
% ratio: the ratio (a/b) of the ellipse
% dir: direction of the

%% Return values:
% R, the rotation matrix from the ellipse frame to the base frame of the
% robot
% theta0: the starting angle of the ellipse (given by the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse
% c: is the center of the ellipse
% error flag: returns true if an error happens, this error is when the
% vector (p-c) is coinsident with the vector (dir).

p=colVec(p);
c=colVec(c);
dir=colVec(dir);

% X dirction vector
i=dir/norm(dir);

u=p-c;
%
k=cross(u,i);
normk=norm(k);

if(normk==0)
    errorFlag=true;
else
    errorFlag=false;
end

k=k/norm(k);
j=cross(k,i);

% rotation matrix from ellipse frame to base frame
R=[i,j,k];
% calculate a,b minimum,maximum radious of the ellipse
x=u'*i; % x coordinate of the starting point of the ellipse in the ellipse frame
y=u'*j; % y coordinate of the starting point of the ellipse in the ellipse frame
a=(x*x+y*y/(ratio*ratio))^0.5;
b=ratio*a;
% calculate theta0, begining angle of the ellipse
theta0=atan2(y,x);

end

function y=colVec(x)
    if size(x,2)==1
        y=x;
    else
        y=x';
    end
end

function [ thetaVec,sVec ] = getEllipseLengthVector( a,b, theta0,theta1 )
%% About:
% This function is used to return the length of the ellipse as a function
% of the angle theta.

%% Inputs:
% theta0, theta1: start/end angle of the arc of the ellipse (from the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse.

thetaSpan = [theta0 theta1];
s0 = 0;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[thetaVec,sVec] = ode45(@(theta,s) ((a*a*sin(theta)*sin(theta)+b*b*cos(theta)*cos(theta))^0.5), thetaSpan, s0,opts);

end

function [ vec ] = ellipseParametricFunction( a,b,theta )
%% About:
% Calculate the X and Y coordinates of the ellipse point corrsponding to
% angle (theta) from the parametric equation

%% inputs:
% theta: angle of point on the ellipse (from the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse

%% Return value:
% The return value is a column vector, vec=[x;y;0] where:
% x: the x coordiante of the point on the ellipse correspodning to angle
% (theta) as described in the frame of the ellipse
% y: the y coordiante of the point on the ellipse  correspodning to angle
% (theta) as described in the frame of the ellipse

x=a*cos(theta);
y=b*sin(theta);
vec=[x;y;0];

end

function [t0,t1,t2]=calculateInterpolationTimes(L,v,a)
%% About:
% calculates the times of the acceleration, stable motion, deceleration

%% Arreguments:
% L: length of the curve (mm).
% v: linear velocity of eef on the the curve (mm/sec).
% a: linear acceleration of eef on the the curve (mm/sec2).

%% Return values:
% t0; time for acceleration.
% t1: time for acceleration, constant velocity.
% t2: time for acceleration, constant velocity and for deceleration

taw=v/a;

if(L>a*taw*taw)
    t0=taw;
    deltat=(L-a*taw*taw)/v;
    t1=t0+deltat;
    t2=t1+taw;
else
    taw1=(L/a)^0.5;
    t0=taw1;
    t1=taw1;
    t2=2*taw1;
end

end


text:
function r = foldr(list, identity, func)
% VAL = FOLDR(ARRAY, ID, FUNC)
%   Perform a right fold on ARRAY using FUNC with identity ID.  FOLDR is
%   the generic recursive list operator.  It applies FUNC to every element
%   of ARRAY in reverse order, along with ID.  The result is an `update' to
%   ID.  The final result of ID is returned.
%
%   FOLDR is most simply viewed as the haskell code
%     foldr []   id f = id
%     foldr x:xs id f = f (foldr xs id f) x
%
%   FOLDL and FOLDR are implemented as a loop for optimization purposes in
%   these libraries, whereas usually FOLDL is tail recursive and thus more
%   efficient on most platforms.
%

  if isempty(list)
    r = identity;
    return;
  end

  if iscell(list)
    for i = numel(list):-1:1
      identity = func(identity, list{i});
    end
  else
    for i = numel(list):-1:1
      identity = func(identity, list(i));
    end
  end

  r = identity;
end


text:
%TIME_7
%
%   The purpose of this test is to compare the FB routines.
%
%   fb_1 is the original, using two passes through memory and no integer
%   optimizations.
%
%   fb_2 use only one pass and integer optimizations.

Lr=[480000*sf^2,480000*sf^2,262144*sf^2,262144*sf^2,900*sf^2];
ar=[     600*sf,     600*sf,        512,        512,       2];
Mr=[     800*sf,     800*sf,       1024,       1024,  600*sf];
gr=[     800*sf,  40*600*sf,       1024,     40*512,  600*sf];

for ii=1:length(Lr)

  L=Lr(ii);

  M=Mr(ii);
  a=ar(ii);

  gl=gr(ii);

  g=rand(L,1);
  f=rand(L,1);
  gfir=rand(gl,1);

  [L, a, M, gl]

  c1=mex_dgt_fb_1(f,gl,a,M,1);
  c2=mex_dgt_fb_2(f,gl,a,M,1);

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function C = cs_permute (A,p,q)                                             %#ok
%CS_PERMUTE permute a sparse matrix.
%   C = cs_permute(A,p,q) computes C = A(p,q)
%
%   Example:
%       Prob = UFget ('HB/arc130') ; A = Prob.A ; [m n] = size (A) ;
%       p = randperm (m) ; q = randperm (n) ;
%       C = cs_permute (A,p,q) ;    % C = A(p,q)
%
%   See also CS_SYMPERM, SUBSREF.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_permute mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function i4vec_transpose_print_test (  )

%*****************************************************************************80
%
%% I4VEC_TRANSPOSE_PRINT_TEST tests I4VEC_TRANSPOSE_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
  n = 12;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_TRANSPOSE_PRINT_TEST\n' );
  fprintf ( 1, '  I4VEC_TRANSPOSE_PRINT prints an integer vector\n' );
  fprintf ( 1, '  with 5 entries to a row, and a title.\n' );

  a = i4vec_indicator1 ( n );

  i4vec_print ( n, a, '  Output from I4VEC_PRINT:' );

  i4vec_transpose_print ( n, a, '  My array:  ' );

  return
end


text:
function ft_create_buffer(port)

% FT_CREATE_BUFFER starts the thread with the TCP server attached to the local
% Matlab instance. The TCP server will listen to the specified network
% port, and accept incoming read and write requests.
%
% Use as
%   ft_create_buffer(port)
% where port is the TCP port to which the server listens. The default port
% number is 1972.
%
% See also FT_DESTROY_BUFFER

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

if nargin<1
  port = 1972;
end

try
  buffer('tcpserver', 'init', 'localhost', port);
  pause(1);
catch
  if ~isempty(strfind(lasterr, 'thread is already running'))
    warning('thread is already running');
  else
    rethrow(lasterror);
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function out = DN_pleft(y,th)
% DN_pleft  Distance from the mean at which a given proportion of data are
%                   more distant.
%
% Measures the maximum distance from the mean at which a given fixed proportion,
% p, of the time-series data points are further.
% Normalizes by the standard deviation of the time series
% (could generalize to separate positive and negative deviations in future)
% Uses the quantile function from Matlab's Statistics Toolbox
%
%---INPUTS:
% y, the input data vector
% th, the proportion of data further than p from the mean
%           (output p, normalized by standard deviation)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if nargin < 2 || isempty(th)
    th = 0.1; % default
end

p = quantile(abs(y-mean(y)),1-th);

% A proportion, th, of the data lie further than p from the mean
out = p/std(y);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%% GPS
% for i = 1:32
%     codes_L1CA(:,i) = GNSScodegen(i,'L1CA',0);
%     codes_L2CM(:,i) = GNSScodegen(i,'L2CM',0);
%     codes_L2CL(:,i) = GNSScodegen(i,'L2CL',0);
%     codes_L5I(:,i) = GNSScodegen(i,'L5I',0);
%     codes_L5Q(:,i) = GNSScodegen(i,'L5Q',0);
% end

% save('codes_L1CA.mat','codes_L1CA');
% save('codes_L2CM.mat','codes_L2CM');
% save('codes_L2CL.mat','codes_L2CL');
% save('codes_L5I.mat','codes_L5I');
% save('codes_L5Q.mat','codes_L5Q');

%% Galileo
% for i = 1:50
%     codes_E5aI(:,i) = GNSScodegen(i,'E5aI',0);
%     codes_E5aQ(:,i) = GNSScodegen(i,'E5aQ',0);
%     codes_E5bI(:,i) = GNSScodegen(i,'E5bI',0);
%     codes_E5bQ(:,i) = GNSScodegen(i,'E5bQ',0);
%     codes_E1B(:,i) = GNSScodegen(i,'E1B',0);
%     codes_E1C(:,i) = GNSScodegen(i,'E1C',0);
% end
%
% save('codes_E5aI.mat','codes_E5aI');
% save('codes_E5aQ.mat','codes_E5aQ');
% save('codes_E5bI.mat','codes_E5bI');
% save('codes_E5bQ.mat','codes_E5bQ');
% save('codes_E1B.mat','codes_E1B');
% save('codes_E1C.mat','codes_E1C');

%% BeiDou
% for i = 1:37
%     codes_B1I(:,i) = GNSScodegen(i,'B1I',0);
% end
%
% save('codes_B1I.mat','codes_B1I');

author: danipascual

repo: GNSS-matlab


text:
function cellcell = BF_cell2cellcell(cellin,delimiter)
% BF_cell2cellcell  Turn a cell of strings into a cell of cells, splitting on a delimiter.
%
% Inputs a cell with some delimiter, and outputs a cell of cells using this
% delimiter.
%
%---INPUTS:
% cellin, the cell
% delimiter, the delimiter
%
%---OUTPUT:
% cellcell, the cell of cells, using this delimiter

% (Used for some tasks involving mySQL)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

if nargin < 2 || isempty(delimiter)
    delimiter = ','; % comma as the default delimiter
end

numElements = length(cellin); % number of elements in the input cell
cellcell = cell(numElements,1);

for i = 1:numElements
    cellcell{i} = regexp(cellin{i}, delimiter, 'split');
end

end


text:
function [bads] = fiff_read_bad_channels(fid,node,ch_rename)
%
% [bads] = fiff_read_bad_channels(fid,node[,ch_rename])
%
% Reas the bad channel list from a node if it exists
%
% fid       - The file id
% node      - The node of interes
% ch_rename - Short-to-long channel name mapping
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.5  2006/05/03 18:53:04  msh
%   Approaching Matlab 6.5 backward compatibility
%
%   Revision 1.4  2006/05/03 15:51:17  msh
%   Re-added fiff_read_bad_channels
%
%   Revision 1.2  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.1  2006/04/12 17:09:28  msh
%   Added routines for reading noise-covariance matrices
%
%

me='MNE:fiff_read_bad_channels';

if nargin == 2
    ch_rename = {};
elseif nargin ~= 3
    error(me,'Incorrect number of arguments');
end

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

node = fiff_dir_tree_find(node,FIFF.FIFFB_MNE_BAD_CHANNELS);

bads = [];
if ~isempty(node)
    tag = find_tag(node,FIFF.FIFF_MNE_CH_NAME_LIST);
    if ~isempty(tag)
        bads = fiff_split_name_list(tag.data);
    end
end
bads = fiff_rename_list(bads, ch_rename);

return;

    function [tag] = find_tag(node,findkind)

        for p = 1:node.nent
            kind = node.dir(p).kind;
            pos  = node.dir(p).pos;
            if kind == findkind
                tag = fiff_read_tag(fid,pos);
                return;
            end
        end
        tag = [];
        return
    end

end

author: brainstorm-tools