text:
% Defines constants
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 5/7/2009
classdef ConstantUtil
    properties (Constant=true)
        nAvogadro=6.0221417930e23;
        elements=struct(...
            'H',1.0079,...
            'He',4.0026,...
            'Li',6.941,...
            'Be',9.0122,...
            'B',10.811,...
            'C',12.0107,...
            'N',14.0067,...
            'O',15.9994,...
            'F',18.9984,...
            'Ne',20.1797,...
            'Na',22.9897,...
            'Mg',24.305,...
            'Al',26.9815,...
            'Si',28.0855,...
            'P',30.9738,...
            'S',32.065,...
            'Cl',35.453,...
            'Ar',39.948,...
            'K',39.0983,...
            'Ca',40.078,...
            'Sc',44.9559,...
            'Ti',47.867,...
            'V',50.9415,...
            'Cr',51.9961,...
            'Mn',54.938,...
            'Fe',55.845,...
            'Co',58.9332,...
            'Ni',58.6934,...
            'Cu',63.546,...
            'Zn',65.39,...
            'Ga',69.723,...
            'Ge',72.64,...
            'As',74.9216,...
            'Se',78.96,...
            'Br',79.904,...
            'Kr',83.8,...
            'Rb',85.4678,...
            'Sr',87.62,...
            'Y',88.9059,...
            'Zr',91.224,...
            'Nb',92.9064,...
            'Mo',95.94,...
            'Tc',98,...
            'Ru',101.07,...
            'Rh',102.9055,...
            'Pd',106.42,...
            'Ag',107.8682,...
            'Cd',112.411,...
            'In',114.818,...
            'Sn',118.71,...
            'Sb',121.76,...
            'Te',127.6,...
            'I',126.9045,...
            'Xe',131.293,...
            'Cs',132.9055,...
            'Ba',137.327,...
            'La',138.9055,...
            'Ce',140.116,...
            'Pr',140.9077,...
            'Nd',144.24,...
            'Pm',145,...
            'Sm',150.36,...
            'Eu',151.964,...
            'Gd',157.25,...
            'Tb',158.9253,...
            'Dy',162.5,...
            'Ho',164.9303,...
            'Er',167.259,...
            'Tm',168.9342,...
            'Yb',173.04,...
            'Lu',174.967,...
            'Hf',178.49,...
            'Ta',180.9479,...
            'W',183.84,...
            'Re',186.207,...
            'Os',190.23,...
            'Ir',192.217,...
            'Pt',195.078,...
            'Au',196.9665,...
            'Hg',200.59,...
            'Tl',204.3833,...
            'Pb',207.2,...
            'Bi',208.9804,...
            'Po',209,...
            'At',210,...
            'Rn',222,...
            'Fr',223,...
            'Ra',226,...
            'Ac',227,...
            'Th',232.0381,...
            'Pa',231.0359,...
            'U',238.0289,...
            'Np',237,...
            'Pu',244,...
            'Am',243,...
            'Cm',247,...
            'Bk',247,...
            'Cf',251,...
            'Es',252,...
            'Fm',257,...
            'Md',258,...
            'No',259,...
            'Lr',262,...
            'Rf',261,...
            'Db',262,...
            'Sg',266,...
            'Bh',264,...
            'Hs',277,...
            'Mt',268);
        elementNumbers=struct(...
            'H', 1,...
            'He',4,...
            'Li',7,...
            'Be',9,...
            'B', 11,...
            'C', 12,...
            'N', 14,...
            'O', 16,...
            'F', 19,...
            'Ne',20,...
            'Na',23,...
            'Mg',24,...
            'Al',27,...
            'Si',28,...
            'P', 31,...
            'S', 32,...
            'Cl',35,...
            'Ar',40,...
            'K', 39,...
            'Ca',40,...
            'Sc',45,...
            'Ti',48,...
            'V', 51,...
            'Cr',52,...
            'Mn',55,...
            'Fe',56,...
            'Co',59,...
            'Ni',59,...
            'Cu',64,...
            'Zn',65,...
            'Ga',70,...
            'Ge',73,...
            'As',75,...
            'Se',79,...
            'Br',80,...
            'Kr',84,...
            'Rb',85,...
            'Sr',87,...
            'Y', 89,...
            'Zr',91,...
            'Nb',93,...
            'Mo',96,...
            'Tc',98,...
            'Ru',101,...
            'Rh',103,...
            'Pd',107,...
            'Ag',108,...
            'Cd',112,...
            'In',115,...
            'Sn',119,...
            'Sb',122,...
            'Te',128,...
            'I', 127,...
            'Xe',131,...
            'Cs',133,...
            'Ba',137,...
            'La',139,...
            'Ce',140,...
            'Pr',141,...
            'Nd',144,...
            'Pm',145,...
            'Sm',150,...
            'Eu',152,...
            'Gd',157,...
            'Tb',159,...
            'Dy',163,...
            'Ho',165,...
            'Er',167,...
            'Tm',169,...
            'Yb',173,...
            'Lu',175,...
            'Hf',178,...
            'Ta',181,...
            'W', 184,...
            'Re',186,...
            'Os',190,...
            'Ir',192,...
            'Pt',195,...
            'Au',197,...
            'Hg',201,...
            'Tl',204,...
            'Pb',207,...
            'Bi',209,...
            'Po',209,...
            'At',210,...
            'Rn',222,...
            'Fr',223,...
            'Ra',226,...
            'Ac',227,...
            'Th',232,...
            'Pa',231,...
            'U', 238,...
            'Np',237,...
            'Pu',244,...
            'Am',243,...
            'Cm',247,...
            'Bk',247,...
            'Cf',251,...
            'Es',252,...
            'Fm',257,...
            'Md',258,...
            'No',259,...
            'Lr',262,...
            'Rf',261,...
            'Db',262,...
            'Sg',266,...
            'Bh',264,...
            'Hs',277,...
            'Mt',268);
        secondsPerMinute = 60;
        secondsPerHour=3600;
        minutesPerHour=60;

        realmax = 1e10;
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function kern = polyardKernParamInit(kern)

% POLYARDKERNPARAMINIT POLYARD kernel parameter initialisation.
% The automatic relevance determination version of the polynomial
% kernel is included for completeness, but its use is generally not
% recommended.
%
%  k(x_i, x_j) = sigma2*(w*x_i'*A*x_j+b)^d
%
% The kernel parameters are sigma2 (kern.variance), w
% (kern.weightVariance), b (kern.biasVariance), A
% (kern.inputScales) a diagonal matrix of input scales and d
% (kern.degree). Only gradients of the first four are provided for
% kernel optimisation, it is assumed that polynomial degree would
% be set by hand.
%
% The kernel is not recommended as it is badly behaved when the
% w*x_i'*A*x_j + b has a magnitude greater than one.
%
% SEEALSO : polyKernParamInit
%
% FORMAT
% DESC initialises the automatic relevance determination polynomial
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% KERN

kern.weightVariance = 1;
kern.biasVariance = 1;
kern.variance = 1;
% These parameters are restricted to lie between 0 and 1.
kern.inputScales = 0.999*ones(1, kern.inputDimension);
kern.nParams = 3 + kern.inputDimension;

kern.degree = 2;

kern.transforms(1).index = [1 2 3];
kern.transforms(1).type = optimiDefaultConstraint('positive');
kern.transforms(2).index = [4:kern.nParams];
kern.transforms(2).type = optimiDefaultConstraint('zeroone');

kern.isStationary = false;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function kts = mph2kts(mph)
%MPH2KTS Convert speed from miles per hour to knots
%
%  kts = MPH2KTS(mph) convert speeds from miles per hour to knots.
%
%  See also MPH2KMPH, MPH2FTPS, MPH2MPS, KTS2MPH.

% Jonathan Sullivan
% Original: May 2011
% jonathan.sullivan@ll.mit.edu

kts = mph*0.868976;

author: Sable

repo: mcbench-benchmarks


text:
function InitializeAppData(hfig)
%% Import paths

data_masterdir = GetCurrentDataDir();%'C:\Janelia2014'; %[pwd '\example data']; % 'F:\Janelia2014';%
save_dir = GetOutputDataDir();

global VAR; %#ok<NUSED>
load(fullfile(save_dir,'VAR_new.mat'),'VAR'); % stores all clustering indices

name_MASKs = 'MaskDatabase.mat';
name_ReferenceBrain = 'ReferenceBrain.mat';

setappdata(hfig,'data_masterdir',data_masterdir);% directory for full fish data (.mat)
setappdata(hfig,'save_dir',save_dir);% directory for full fish data (.mat)

%% load ZBrain Atlas
if exist(fullfile(data_masterdir,name_MASKs),'file') ...
        && exist(fullfile(data_masterdir,name_ReferenceBrain),'file'),

    % load masks for ZBrain Atlas
    MASKs = load(fullfile(data_masterdir,name_MASKs));
    setappdata(hfig,'MASKs',MASKs);
    % load reference brain image stack with the 3 projections
    load(fullfile(data_masterdir,name_ReferenceBrain));
    setappdata(hfig,'anat_stack_norm',anat_stack_norm);
    setappdata(hfig,'anat_yx_norm',anat_yx_norm);
    setappdata(hfig,'anat_yz_norm',anat_yz_norm);
    setappdata(hfig,'anat_zx_norm',anat_zx_norm);

    FishOutline = load(fullfile(data_masterdir,'FishOutline.mat'));
    setappdata(hfig,'FishOutline',FishOutline);
end

%% preload
% fish protocol sets (different sets have different parameters)
M_fish_set = GetFishStimset(); % M = Matrix
setappdata(hfig,'M_fish_set',M_fish_set);

% constant
setappdata(hfig,'z_res',19.7); % resoltion in z, um per slice

%% load target fish
% need to reload fish if changed
setappdata(hfig,'i_fish',0);
setappdata(hfig,'isFullData',1);

%% display params
% for general display
setappdata(hfig,'isPopout',1);
setappdata(hfig,'isPlotAnatomyOnly',0);
setappdata(hfig,'clrmap_name','hsv_new');

% for functional display: need to UpdateTimeIndex if changed
setappdata(hfig,'isStimAvr',1); % show average/full stimulus
setappdata(hfig,'isRawtime',0); % show stimulus in original order or sorted
setappdata(hfig,'isZscore',1); % show normalized (z-score) version of fluorescent
setappdata(hfig,'isTrialRes',0); % show trial residuals (subtracting trial average) of fluorescent
% % temp:
setappdata(hfig,'isClusRes',0);

% for functional display: display only
setappdata(hfig,'isCentroid',0);
setappdata(hfig,'isPlotLines',0);
setappdata(hfig,'isPlotBehavior',1);
setappdata(hfig,'isPlotRegWithTS',0);
setappdata(hfig,'rankscore',[]); %
setappdata(hfig,'rankID',0); % whether to plot scores or not, toggled away

% for anat display:
setappdata(hfig,'isRefAnat',0);
setappdata(hfig,'isShowFishOutline',0);
setappdata(hfig,'isShowMasks',1);
setappdata(hfig,'isShowMskOutline',0);
setappdata(hfig,'isWeighAlpha',0); % for regression plots, toggled away

%% for behavior / motor regression:
setappdata(hfig,'isMotorseed',1);

%% all below: associated with specialized functions
% set operations
setappdata(hfig,'opID',0);

% regression
thres_merge = 0.9;
thres_split = 0.7;
thres_reg = 0.7;
thres_size = 10;
thres_ttest = 0.001;
setappdata(hfig,'thres_merge',thres_merge);
setappdata(hfig,'thres_split',thres_split); % for function 'pushbutton_iter_split'
setappdata(hfig,'thres_reg',thres_reg); % regression threshold, ~correlation coeff
setappdata(hfig,'thres_size',thres_size); % min size for clusters

setappdata(hfig,'thres_ttest',thres_ttest); % min size for clusters

setappdata(hfig,'regchoice',{1,1}); % regressor choice; first cell,1=stim,2=motor,3=centroid
setappdata(hfig,'isPlotReg',1); % plot regressor when selecting it
setappdata(hfig,'isPlotCorrHist',0); % option for regression

setappdata(hfig,'isRegIndividualCells',1);
setappdata(hfig,'isRegCurrentCells',1);

setappdata(hfig,'isAutoclusWithAllCells',1);

% clustering
setappdata(hfig,'isWkmeans',1); % in autoclustering, with/without kmeans
setappdata(hfig,'hierinplace',1); % hier. partitioning, no reordering
setappdata(hfig,'isCorrplotWeighted',0);

% Atlas
setappdata(hfig,'isFindMaskNorm',1);
setappdata(hfig,'isPlotMskHist',1);
setappdata(hfig,'isScreenMskFromAllCells',0);

% GUI interface helpers
setappdata(hfig,'clusgroupID_view',1);
setappdata(hfig,'clusID_view',0); % set in UpdateClusGroupGUI, 0 for placeholder '(choose)'

end

author: xiuyechen

repo: FishExplorer


text:
function ret = my_gpu_conv2(l_prev_layer, filter)
    ret = conv2(l_prev_layer,filter,'valid');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function test_suite = test_createRotation
%TESTCREATEROTATION  One-line description here, please.
%   output = testCreateRotation(input)
%
%   Example
%   testCreateRotation
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.
% Licensed under the terms of the LGPL, see the file "license.txt"

test_suite = functiontests(localfunctions);

function testCreateCentered(testCase) %#ok<*DEFNU>

trans = createRotation(0);
testCase.assertEqual(trans, [1 0 0;0 1 0;0 0 1], 'AbsTol', .01);

trans = createRotation(pi/2);
testCase.assertEqual(trans, [0 -1 0; 1 0 0; 0 0 1], 'AbsTol', .01);

trans = createRotation(pi);
testCase.assertEqual(trans, [-1 0 0;0 -1 0;0 0 1], 'AbsTol', .01);

trans = createRotation(3*pi/2);
testCase.assertEqual(trans, [0 1 0; -1 0 0; 0 0 1], 'AbsTol', .01);

function testCreateShifted(testCase)

p0 = [3 5];
theta = pi/3;

trans1 = createRotation(p0, theta);
t1 = createTranslation(-p0);
rot = createRotation(theta);
t2 = createTranslation(p0);
trans2 = t2*rot*t1;

testCase.assertEqual(trans1, trans2, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function varargout = subscript(m,range,wantc)

% function varargout = subscript(m,range,wantc)
%
% <m> is a matrix (of numbers) or a string referring to a matrix
%   variable in the base workspace.  cell matrices are okay.
% <range> is:
%   (1) vector index range that does not use
%       the 'end' keyword, e.g., 4, [5 6]
%   (2) a logical indexing matrix
%   (3) the string ':'
%   (4) a cell vector where elements can be of types (1), (2), (3),
%       e.g., {':' 4:5}
% <wantc> (optional) is whether to perform cell de-referencing.  default: 0.
%
% return something like m(range), or m{range} when <wantc>.
%
% this function is useful for cases where you want
% to get a range of elements from something that
% already has parentheses or brackets.
% it is also useful for working with string-variable
% representations of matrices.  it is also useful
% for writing functions that make no assumption about
% the size of the matrices used as input.
%
% example:
% isequal(size(subscript(randn(10,10),{1:2 ':'})),[2 10])
% a = [1 2; 3 4];
% subscript('a',1)
% subscript('a',a<1.5)

% input
if ~exist('wantc','var') || isempty(wantc)
  wantc = 0;
end

% do it
if wantc
  if iscell(range)
    if ischar(m)
      varargout = evalin('base',['subscript(',m,',',cell2str(range),',1);']);
    else
      varargout = m(range{:});
    end
  else
    if ischar(m)
      if ischar(range) && isequal(range,':')
        varargout = evalin('base',[m,'(:);']);
      else
        varargout = evalin('base',[m,'(',mat2str(range),');']);
      end
    else
      if ischar(range) && isequal(range,':')
        varargout = m(:);
      else
        varargout = m(range);
      end
    end
  end
else
  varargout = cell(1,1);
  if iscell(range)
    if ischar(m)
      varargout{1} = evalin('base',['subscript(',m,',',cell2str(range),');']);
    else
      varargout{1} = m(range{:});
    end
  else
    if ischar(m)
      if ischar(range) && isequal(range,':')
        varargout{1} = evalin('base',[m,'(:);']);
      else
        varargout{1} = evalin('base',[m,'(',mat2str(range),');']);
      end
    else
      if ischar(range) && isequal(range,':')
        varargout{1} = m(:);
      else
        varargout{1} = m(range);
      end
    end
  end
end

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function tr=tau_rest(veh,p)

% tr=tau_rest(veh,p); calculates restoring forces from
% vehicle variables and generalized position p

% Hydrostatic force and moment
FB_e=-veh.vol*veh.rho*veh.g_e;
FB_b=rpy2R_eb(p(4:6))*FB_e;
MB_b=vp(veh.B_b,FB_b);
tb=[FB_b;MB_b];

% Gravitational force and moment
FG_e=veh.m*veh.g_e;
FG_b=rpy2R_eb(p(4:6))*FG_e;
MG_b=vp(veh.G_b,FG_b);
tg=[FG_b;MG_b];

tr=tb+tg;


text:
%% The Class EBSDhex
%
%%
%

author: mtex-toolbox

repo: mtex

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


text:
function [cMap]=rgbImage2ColorMap(X,optStruct)

% function [cMap]=rgbImage2ColorMap(X,optStruct)
% ------------------------------------------------------------------------
%
% This function creates a colormap based on an input RGB (red green blue)
% color image (X). Colormap colors are based on the intensities occuring in
% the image. Therefore monotonic color images probably yeild the best
% results.
% The colormap is harvested from the input image using the settings defined
% in optStruct. The latter contains the fields:
% n: The number of colormap levels (default 250)
% normFactor: The normalisation factor e.g. 255 (default is maximum value
% occuring in image).
% numBins: Number of bins to sample the image with (colors are averaged for
% each bin, default=n)
%
% See also: |colormap|
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2015/04/01 Added to GIBBON
%------------------------------------------------------------------------

%% Parse input

%Get options from structure or use defaults

if isfield(optStruct,'n')
   n=optStruct.n;
else
    n=250;
end

if isfield(optStruct,'normFactor')
   normFactor=optStruct.normFactor;
else
    normFactor=max(X(:));
end

if isfield(optStruct,'numBins')
   numBins=optStruct.numBins;
else
    numBins=n;
end

if isfield(optStruct,'imageIntensityLabel')
   imageIntensityLabel=optStruct.imageIntensityLabel;
else
    imageIntensityLabel=mean(X,3);
end

%Convert image to a double if it isn't already
if ~isa(X,'double')
    X=double(X);
end

%% Sample colors for desired number of intensity levels

%Compute normalised intensity in numBins+1 integers
imageIntensityLabel=imageIntensityLabel-min(imageIntensityLabel(:));
imageIntensityLabel=imageIntensityLabel./max(imageIntensityLabel(:));
imageIntensityLabel=imageIntensityLabel.*numBins;
imageIntensityLabel=round(imageIntensityLabel);

intAll=unique(imageIntensityLabel(:))';

%Get color layers
R=X(:,:,1);
G=X(:,:,2);
B=X(:,:,3);

cMap_sub=zeros(numel(intAll),3);
for q=1:1:numel(intAll)
    indNow=find(imageIntensityLabel==intAll(q));
    cMap_sub(q,:)=[mean(R(indNow)) mean(G(indNow)) mean(B(indNow))];
end

%% Interpolate to evenly across intensities with n levels

intRange=linspace(0,numBins,n)';

cMap=zeros(n,3);
for q=1:1:3;
    cMap(:,q)=interp1(intAll,cMap_sub(:,q),intRange,'linear');
end

%% Normalise colors
cMap=cMap./normFactor;

%%
if any(cMap(:)>1) || any(cMap(:))<0
    error('Colormapped values should be in the range [0 1]. Alter input image and/or normFactor');
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ node_num, node_dim, node_att_num, node_marker_num ] = ...
  triangle_node_size_example ( )

%*****************************************************************************80
%
%% TRIANGLE_NODE_SIZE_EXAMPLE returns the sizes of node information for the example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer NODE_DIM, the spatial dimension.
%
%    Output, integer NODE_ATT_NUM, number of node attributes
%    listed on each node record.
%
%    Output, integer NODE_MARKER_NUM, 1 if every node record
%    includes a final boundary marker value.
%
  node_num = 21;
  node_dim = 2;
  node_att_num = 0;
  node_marker_num = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function handle = imageModify(handle, imageValues, imageSize, transpose, negative, ...
		     scale)

% IMAGEMODIFY Helper code for visualisation of image data.
% FORMAT
% DESC is a helper function for visualising image data using latent
% variable models.
% ARG handle : the handle of the image data.
% ARG imageValues : the values to set the image data to.
% ARG imageSize : the size of the image.
% ARG transpose : whether the resized image needs to be transposed
% (default 1, which is yes).
% ARG negative : whether to display the negative of the image
% (default 0, which is no).
% ARG scale : dummy input, to maintain compatability with
% IMAGEVISUALISE.
% RETURN handle : a the handle to the image data.
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2006
%
% SEEALSO : imageVisualise, fgplvmResultsDynamic

% MLTOOLS

if nargin < 4
  transpose = 1;
end
if nargin< 5
  negative = 0;
end
if negative
  imageValues = -imageValues;
end
if transpose
  set(handle, 'CData', reshape(imageValues(1:imageSize(1)*imageSize(2)), imageSize(1), imageSize(2))');
else
  set(handle, 'CData', reshape(imageValues(1:imageSize(1)*imageSize(2)), imageSize(1), imageSize(2)));
end


text:
% TRIMMED_GRASSMANN_AVERAGE    Robustly estimate average subspace spanned by data by element-wise trimmed averages.
%
%     TRIMMED_GRASSMANN_AVERAGE(X, P) returns a basis vector for the
%     average one-dimensional subspace spanned by the data X. This is
%     a N-by-D matrix containing N observations in R^D. Robustness is
%     achieved by element-wise trimming at P percent.
%
%     TRIMMED_GRASSMANN_AVERAGE(X, P, K) returns a D-by-K matrix of
%     orthogonal basis vectors spanning a K-dimensional average subspace.
%
%     Reference:
%     "Grassmann Averages for Scalable Robust PCA".
%     S. Hauberg, A. Feragen and M.J. Black. In CVPR 2014.
%     http://ps.is.tue.mpg.de/project/Robust_PCA

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

function vectors = trimmed_grassmann_average(X, percent, K)

  %% Check input
  if (nargin < 2)
    error ('trimmed_grassmann_average: not enough input arguments');
  end % if
  if (nargin < 3)
    K = 1;
  end % if
  if (percent < 0 || percent > 50)
    error('trimmed_grassmann_average: trimming percentage must be between 0 and 50');
  end % if

  %% Initialisation
  [N, D] = size(X);
  vectors = NaN(D, K, class(X));

  %% Options
  epsilon = 1e-5;

  %% Create output structure
  for k = 1:K
    %% Compute k'th principal component
    mu = rand(D, 1) - 0.5;
    mu = mu / norm(mu);

    %% Initialize using a few EM iterations
    for iter = 1:3
      dots = X * mu; % Nx1
      mu = dots.' * X; % 1xD
      mu = mu(:) / norm(mu); % Dx1
    end % for
    dots = []; % clear up memory

    %% Now the Grassmann average
    for iter = 1:N
      prev_mu = mu;

      %% Compute angles and flip
      dot_signs = sign(X * mu); % Nx1

      %% Compute weighted Grassmannian mean
      mu = trimmed_mean(bsxfun(@times, X, dot_signs), percent); % 1xD
      mu = mu(:) / norm(mu); % Dx1

      %% Check for convergence
      if (max(abs(mu - prev_mu)) < epsilon)
        iter
        break;
      end % if
    end % for
    dot_signs = []; % clear up memory
    prev_mu = []; % clear up memory

    %% Store the estimated vector (and possibly subtract it from data, and perform reorthonormalisation)
    if (k == 1)
      vectors(:, k) = mu;
      X = X - (X * mu) * mu.';
    elseif (k < K)
      mu = reorth(vectors(:, 1:k-1), mu, 1);
      mu = mu / norm(mu);
      vectors(:, k) = mu;

      X = X - (X * mu) * mu.';
    else % k == K
      mu = reorth(vectors(:, 1:k-1), mu, 1);
      mu = mu / norm(mu);
      vectors(:, k) = mu;
    end % if
  end % for
end % function

author: andrewssobral

repo: lrslibrary


text:
function output = recall(tp,fn)
%
%
% Recall = tp/(tp+fn)
% (see page 268 of Manning and Schutze)
%
% Inputs
%    tp: True Positive
%    fn: False Negative
% Outputs
%    recall measure

output = tp/(tp+fn);

end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [coords, data] = dtiRoiGrow2d(data, xformToAcpc, seedAcpc, ax, tol, thick)
%
% [selectedCoords, data] = dtiRoiGrow2d(data, xformToAcpc, seedAcpc, ax, tol, [thick=1])
%
% data is either an image volume (scalar or dt6) or the data struct
% returned after calling this function once (useful if you want to iterate
% with the same image volume, just chaning the tolerance).
%
% seedAcpc: the seed point, in ac-pc coords
%
% ax: the axis to grow in (1,2 or 3). E.g., ax=1 will grow in the slice
% X=seedAcpc(1), ax=2 will grow in Y=seedAcpc(2), and ax=3 in
% Z=seedAcpc(3).
%
% tol: the tolerance (between 0 and 1, with lower values including fewer
% voxels).
%
% thick: the thickness of the resulting coordinates. The algorithm grows in
% a single slice. thick > 1 will replicate the selected coords in adjacent
% slices for you.
%
% E.g.:
% dt = dtiLoadDt6('dti06trilinrt/dt6');
% seedAcpc = [-36 -28 26];
% ax = 3;
% tol = .20;
% [coords, data] = dtiRoiGrow2d(dt.dt6, dt.xformToAcpc, seedAcpc, ax, tol);
% % Show the delta image with the selected coords overlaid
% figure; image(data.x(:), data.y(:), data.img); axis equal tight xy;
% hold on; plot(coords(:,1),coords(:,2),'k.'); hold off;
%
% % do it again with stricter tol and passing the pre-processed data:
% coords = dtiRoiGrow2d(data, dt.xformToAcpc, seedAcpc, ax, 0.10);
% figure; image(data.x(:), data.y(:), data.img); axis equal tight xy;
% hold on; plot(coords(:,1),coords(:,2),'k.'); hold off;
%
% HISTORY
% 2009.08.26 RFD: pulled code from dtiFiberUI so that we can easily scrip it.
%

if(~exist('thick','var')||isempty(thick))
    thick = 1;
end

if(~isstruct(data))
    if(ndims(data)==3)
        [img,x,y,z] = dtiGetSlice(xformToAcpc, data, ax, seedAcpc(ax), [], 'n');
        if(ax==2), yseed = find(x(:,1)>seedAcpc(1)-.5&x(:,1)<=seedAcpc(1)+.5);
        else yseed = find(y(:,1)>seedAcpc(2)-.5&y(:,1)<=seedAcpc(2)+.5); end
        if(ax==3), xseed = find(x(1,:)>seedAcpc(1)-.5&x(1,:)<=seedAcpc(1)+.5);
        else xseed = find(z(1,:)>seedAcpc(3)-.5&z(1,:)<=seedAcpc(3)+.5); end
    elseif(ndims(data)==4)
        % Then it's a dt6 (tensor). We'll create a map of the angle between the
        % seed point PDD and each image voxel PDD.
        [dt6,x,y,z] = dtiGetSlice(xformToAcpc, data, ax, seedAcpc(ax), [], 'n');
        if(ax==2), yseed = find(x(:,1)>seedAcpc(1)-.5&x(:,1)<=seedAcpc(1)+.5);
        else yseed = find(y(:,1)>seedAcpc(2)-.5&y(:,1)<=seedAcpc(2)+.5); end
        if(ax==3), xseed = find(x(1,:)>seedAcpc(1)-.5&x(1,:)<=seedAcpc(1)+.5);
        else xseed = find(z(1,:)>seedAcpc(3)-.5&z(1,:)<=seedAcpc(3)+.5); end
        [eigVec,eigVal] = dtiEig(dt6);
        seedPdd = squeeze(eigVec(yseed,xseed,:,1));
        % dot(eigVec(1,:),[1 0 0]) = eigVec(1,1) and dot(eigVec(1,:),[0 1 0]) = eigVec(1,2)
        img = eigVec(:,:,1,1).*seedPdd(1)+eigVec(:,:,2,1).*seedPdd(2)+eigVec(:,:,3,1).*seedPdd(3);
        img(img>1) = 1; img(img<-1) = -1;
        img = acos(img);
        % Reflect about pi/2 for angles > pi/2 (diffusion is symmetric along the eigenvector axis)
        img(img>pi/2) = pi/2-(img(img>pi/2)-pi/2);
        img = img./(pi/2);
        % FA an MD should matter too, since we usually don't want to grow
        % into regions with very different anisotropy from our seed. But-
        % should it get equal weighting to our angle measure? Maybe make it a
        % parameter?
        [fa,md] = dtiComputeFA(eigVal);
        maxMd = md(yseed,xseed)*3;
        md(md>maxMd) = maxMd;
        md = md./maxMd;
        %img = 0.5*img + 0.5*abs(fa-fa(yseed,xseed));
        img = cat(3, img, abs(fa-fa(yseed,xseed)), abs(md-md(yseed,xseed)));
        %figure;imagesc(img);axis image xy; colormap gray
    end
    clear data;
    data.img = img;
    data.x = x;
    data.y = y;
    data.z = z;
    data.xseed = xseed;
    data.yseed = yseed;
end

binImg = magicwand1(data.img, data.yseed, data.xseed, tol);
ind = find(binImg);
coords = [data.x(ind), data.y(ind), data.z(ind)];
for(ii=[1:thick-1,-1:-1:-(thick-1)])
	newLayer = coords;
    newLayer(:,ax) = newLayer(:,ax)+ii;
    coords = [coords; newLayer];
end
coords = unique(coords, 'rows');
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out = trigcoeffs(f, N)
%TRIGCOEFFS   Trigonometric Fourier coefficients of a CLASSICFUN.
%   C = TRIGCOEFFS(F) returns the trigonometric Fourier coefficients of F
%   using complex-exponential form.  Specifically, for N = length(F)
%   If N is odd
%       F(x) = C(1)*z^(N-1)/2 + C(2)*z^((N-1)/2-1) + ... + C((N+1)/2) + ...
%                + C(N)*z^(-(N-1)/2)
%   If N is even
%       F(x) = C(1)*z^(N/2-1) + C(2)*z^(N/2-2) + ... + C(N/2) + ...
%                + C(N-1)*z^(-N/2-1) + 1/2*C(N)*(z^(N/2) + z^(-N/2))
%   where z = exp(1i*pi*x).
%
%   A = TRIGCOEFFS(F, N) truncates or pads the vector C so that N coefficients
%   of F are returned.
%
%   If F is array-valued with M columns, then C is an MxN matrix.
%
% See also LEGCOEFFS, CHEBCOEFFS.

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

if ( nargin == 1 )
    N = length(f);
end

% Call TRIGCOEFFS() of the .ONEFUN:
out = trigcoeffs(f.onefun, N);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [depth_map_in_meters, is_disparity_invalid] =...
    depth_in_meters_cityscapes_with_invalid_parts(input_disparity,...
    camera_parameters_file)
%DEPTH_IN_METERS_CITYSCAPES_WITH_INVALID_PARTS  Compute depth map in meters from
%provided Cityscapes disparity map and relevant camera parameters.
%
%   INPUTS:
%
%   -|input_disparity|: matrix of uint16 format with same resolution as
%   Cityscapes RGB images.
%
%   -|camera_parameters_file|: full path to JSON file where camera parameters
%   are stored.
%
%   OUTPUTS:
%
%   -|depth_map_in_meters|: matrix of double format containing depth in meters.
%   It contains invalid measurements that are known to be such and are set to
%   infinity by convention, and may also contain more erroneous values that are
%   not known to be such and can assume any positive value.
%
%   -|is_disparity_invalid|: boolean matrix, where true indicates that the
%   corresponding pixel in |depth_map_in_meters| contains a wrong infinity
%   value.

% Identify known wrong values for disparity, based on specifications provided in
% the Cityscapes README.
is_disparity_invalid = input_disparity == 0;

% Compute the disparity in pixels, based on specifications provided in the
% Cityscapes README.
disparity_in_pixels = disparity_in_pixels_cityscapes(input_disparity);
is_disparity_zero = disparity_in_pixels == 0;

% Retrieve baseline and focal length in x-axis, which are both required in order
% to get the absolute scale of the depth map.
[B, f_x] = camera_parameters_cityscapes(camera_parameters_file);

% Compute the depth as inversely proportional to disparity. Wherever the
% disparity is zero, the depth is equal to infinity. Convention: any known
% invalid disparity is assigned infinite depth.
depth_map_in_meters = zeros(size(disparity_in_pixels));
depth_map_in_meters(~is_disparity_invalid & ~is_disparity_zero) =...
    B * f_x ./ disparity_in_pixels(~is_disparity_invalid & ~is_disparity_zero);
depth_map_in_meters(is_disparity_zero | is_disparity_invalid) = Inf;

end

author: sakaridis


text:
function [ res ] = crossprod( a, b )
%CROSSPROD Summary of this function goes here
%   Detailed explanation goes here

    Ax = [  0,   -a(3),   a(2);
           a(3),     0,  -a(1);
          -a(2),  a(1),      0;];

   res = Ax * b(:);

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function varargout = wavf_gen(varargin)
% WAVF_GEN M-file for wavf_gen.fig
%      WAVF_GEN, by itself, creates a new WAVF_GEN or raises the existing
%      singleton*.
%
%      H = WAVF_GEN returns the handle to a new WAVF_GEN or the handle to
%      the existing singleton*.
%
%      WAVF_GEN('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in WAVF_GEN.M with the given input arguments.
%
%      WAVF_GEN('Property','Value',...) creates a new WAVF_GEN or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before wavf_gen_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to wavf_gen_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help wavf_gen

% Last Modified by GUIDE v2.5 02-Apr-2009 18:21:03

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @wavf_gen_OpeningFcn, ...
                   'gui_OutputFcn',  @wavf_gen_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before wavf_gen is made visible.
function wavf_gen_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to wavf_gen (see VARARGIN)

% Choose default command line output for wavf_gen
handles.output = hObject;

set(handles.axes1,'NextPlot','add');

global r xys sp iscnt p mth res hc f

f=440; % initial frequency
T=1/f;

mth='spline'; % intial method is spline

res=200; % curve resolution

iscnt=get(handles.cnt,'value');

r=0.015*T;

sp=1*r; % span between markers

% xys(:,1)=[0.5;
%           0.5];
% xys(:,2)=[0.5;
%           0.7];

% p(1)=create_marker(handles.axes1,xys(1,1),xys(2,1),r,[0 1 0],'a',1);
% p(2)=create_marker(handles.axes1,xys(1,2),xys(2,2),r,[1 0 0],'a',2);
nm=str2num(get(handles.nm,'string'));
p=zeros(1,nm);
%p=zeros(1,2);
xys=zeros(2,nm);
for n=1:nm
    stp=T/nm; % step
    xys(1,n)=stp/2+stp*(n-1);
    p(n)=create_marker(handles.axes1,xys(1,n),xys(2,n),r,[0 1 0],'a',n);
end

set(handles.axes1,'XLim',[0 T]);
set(handles.axes1,'YLim',[-1 1]);

% % axis equal:
% set(handles.axes1,'DataAspectRatio',[2 1 1]);
% set(handles.axes1,'DataAspectRatioMode','manual');
set(handles.axes1,'PlotBoxAspectRatio',[6 4 4]);
set(handles.axes1,'PlotBoxAspectRatioMode','manual');

for n=1:nm
    make_round(p(n),handles.axes1);
end

% get(handles.axes1,'DataAspectRatio')

% draw curve:
xl=get(handles.axes1,'Xlim');
yl=get(handles.axes1,'Ylim');
dxl=xl(2)-xl(1);
dyl=yl(2)-yl(1);

xc=xl(1):(dxl/res):xl(2);
yc=interp1(xys(1,:),xys(2,:),xc,mth,'extrap');
handles.cr=plot(xc,yc,'k-');
hc=handles.cr;
set(handles.cr,'hittest','off');

drawnow;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes wavf_gen wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = wavf_gen_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function nm_Callback(hObject, eventdata, handles)

global p xys r iscnt sp res hc mth

% delete old:
for pc=1:length(p)
    delete(p(pc));
end

ha=handles.axes1;
xl=get(ha,'Xlim');
dxl=xl(2)-xl(1);
yl=get(ha,'Ylim');

% create new:
nm=str2num(get(handles.nm,'string'));
p=zeros(1,nm);
%p=zeros(1,2);
xys=zeros(2,nm);
if iscnt
    stp=dxl/nm; % step
    for n=1:nm
        xys(1,n)=stp/2+stp*(n-1);
        p(n)=create_marker(handles.axes1,xys(1,n),xys(2,n),r,[0 1 0],'a',n);
    end
else
    stp=dxl/(nm-1); % step
    for n=1:nm
        xys(1,n)=stp*(n-1);
        p(n)=create_marker(handles.axes1,xys(1,n),xys(2,n),r,[0 1 0],'a',n);
    end
end

% set(handles.axes1,'XLim',[0 3]);
% set(handles.axes1,'YLim',[-1 1]);

% % axis equal:
% set(handles.axes1,'DataAspectRatio',[2 1 1]);
% set(handles.axes1,'DataAspectRatioMode','manual');
% set(handles.axes1,'PlotBoxAspectRatio',[6 4 4]);
% set(handles.axes1,'PlotBoxAspectRatioMode','manual');

for n=1:nm
    make_round(p(n),handles.axes1);
end

% update curve:
update_curve;

% --- Executes during object creation, after setting all properties.
function nm_CreateFcn(hObject, eventdata, handles)
% hObject    handle to nm (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in cnt.
function cnt_Callback(hObject, eventdata, handles)

global iscnt xys p sp res hc mth

iscnt=get(hObject,'Value');

ha=handles.axes1;
xl=get(ha,'Xlim');
dxl=xl(2)-xl(1);
yl=get(ha,'Ylim');

if iscnt
    % if from not continius to continius

    xdt=get(p(end),'XData');

    nm=length(xys(1,:));
    if (dxl-xys(1,nm-1))>sp
        dx=xl(2)-sp-xys(1,end);
    else
        sp1=(xys(1,nm)-xys(1,nm-1))/2;
        dx=xl(2)-sp1-xys(1,end);
    end

    xdt=xdt+dx;
    set(p(end),'XData',xdt);
    xys(1,end)=xys(1,end)+dx;

else
    % if from continius to not continius

    xdt=get(p(1),'XData');
    dx=xl(1)-xys(1,1);
    xdt=xdt+dx;
    set(p(1),'XData',xdt);
    xys(1,1)=xl(1);

    xdt=get(p(end),'XData');
    dx=xl(2)-xys(1,end);
    xdt=xdt+dx;
    set(p(end),'XData',xdt);
    xys(1,end)=xl(end);

    %xys(1,end)=xl(2);

end

update_curve;

% --- Executes when uipanel1 is resized.
function uipanel1_ResizeFcn(hObject, eventdata, handles)
% hObject    handle to uipanel1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes when selected object is changed in uipanel1.
function uipanel1_SelectionChangeFcn(hObject, eventdata, handles)

global mth hc xys res iscnt

if get(handles.nearest,'value')
    mth='nearest';
end

if get(handles.linear,'value')
    mth='linear';
end

if get(handles.spline,'value')
    mth='spline';
end

xl=get(handles.axes1,'Xlim');
yl=get(handles.axes1,'Ylim');
dxl=xl(2)-xl(1);
dyl=yl(2)-yl(1);
update_curve;

function fr_Callback(hObject, eventdata, handles)

global f p xys r sp res mth hc iscnt

f_old=f; % old value of f
T_old=1/f_old;

frs=get(handles.fr,'string');
f=str2num(frs);
T=1/f;

xys(1,:)=T*xys(1,:)/T_old; % scretch data as T

r=0.015*T;

sp=1*r; % span between markers

set(handles.axes1,'XLim',[0 T]);
set(handles.axes1,'YLim',[-1 1]);

nm=length(xys(1,:));

for n=1:nm
    make_round(p(n),handles.axes1);
end

xl=get(handles.axes1,'Xlim');
yl=get(handles.axes1,'Ylim');
dxl=xl(2)-xl(1);
dyl=yl(2)-yl(1);
update_curve;

% --- Executes during object creation, after setting all properties.
function fr_CreateFcn(hObject, eventdata, handles)
% hObject    handle to fr (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function sl_Callback(hObject, eventdata, handles)
% hObject    handle to sl (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of sl as text
%        str2double(get(hObject,'String')) returns contents of sl as a double

% --- Executes during object creation, after setting all properties.
function sl_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sl (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pl.
function pl_Callback(hObject, eventdata, handles)

global xys iscnt mth f

prepare_signal;

sound(yc,Fs); % play

% --- Executes on button press in sv.
function sv_Callback(hObject, eventdata, handles)

[filename, pathname] = uiputfile(...
 {'*.wav'},...
 'Save sound as wav file');

filename1=[pathname filename];

global xys iscnt mth f

prepare_signal;

wavwrite(yc,Fs,filename1);

% --- Executes on selection change in prs.
function prs_Callback(hObject, eventdata, handles)
% hObject    handle to prs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns prs contents as cell array
%        contents{get(hObject,'Value')} returns selected item from prs

global p xys r iscnt sp res hc mth f

spr=get(hObject,'Value'); % selected preset

T=1/f;

if spr~=1
    % delete old:
    for pc=1:length(p)
        delete(p(pc));
    end

    ha=handles.axes1;
    xl=get(ha,'Xlim');
    dxl=xl(2)-xl(1);
    yl=get(ha,'Ylim');
end

% create new:

switch spr
    case 1 % nothing

    case 2 % zero
        np=5; % number of points
        set(handles.nm,'string',num2str(np));
        stp=T/np;
        p=zeros(1,np);
        xys=zeros(2,np);
        for n=1:np
            xys(1,n)=(n-1)*stp;
            xys(2,n)=0;
        end

        % set continius:
        iscnt=true;
        set(handles.cnt,'value',iscnt);

        % set spline method:
        mth='spline';
        set(handles.nearest,'value',0);
        set(handles.linear,'value',0);
        set(handles.spline,'value',1);
    case 3 % sine
        snrs=8; % sine resolution
        set(handles.nm,'string',num2str(snrs));
        stp=T/snrs;
        p=zeros(1,snrs);
        xys=zeros(2,snrs);
        for n=1:snrs
            xys(1,n)=(n-1)*stp;
            xys(2,n)=0.9*sin(2*pi*f*xys(1,n));
        end

        % set continius:
        iscnt=true;
        set(handles.cnt,'value',iscnt);

        % set spline method:
        mth='spline';
        set(handles.nearest,'value',0);
        set(handles.linear,'value',0);
        set(handles.spline,'value',1);

    case 4 % square
        np=2; % number of points
        set(handles.nm,'string',num2str(np));
        stp=T/np;
        p=zeros(1,np);
        xys=zeros(2,np);
        for n=1:np
            xys(1,n)=(n-1)*stp+T/4;
            if n==1
                xys(2,n)=-0.9;
            else
                xys(2,n)=0.9;
            end
        end

        % set continius:
        iscnt=true;
        set(handles.cnt,'value',iscnt);

        % set nearest method:
        mth='nearest';
        set(handles.nearest,'value',1);
        set(handles.linear,'value',0);
        set(handles.spline,'value',0);

    case 5 % triangle
        np=2; % number of points
        set(handles.nm,'string',num2str(np));
        stp=T/np;
        p=zeros(1,np);
        xys=zeros(2,np);
        for n=1:np
            xys(1,n)=(n-1)*stp+T/4;
            if n==1
                xys(2,n)=-0.9;
            else
                xys(2,n)=0.9;
            end
        end

        % set continius:
        iscnt=true;
        set(handles.cnt,'value',iscnt);

        % set linear method:
        mth='linear';
        set(handles.nearest,'value',0);
        set(handles.linear,'value',1);
        set(handles.spline,'value',0);

    case 6 % sawtooth
        np=2; % number of points
        set(handles.nm,'string',num2str(np));
        stp=T/np;
        p=zeros(1,np);
        xys=zeros(2,np);

        xys(1,1)=0;
        xys(2,1)=-0.9;

        xys(1,2)=T;
        xys(2,2)=0.9;

        % set not continius:
        iscnt=false;
        set(handles.cnt,'value',iscnt);

        % set linear method:
        mth='linear';
        set(handles.nearest,'value',0);
        set(handles.linear,'value',1);
        set(handles.spline,'value',0);
end

if spr~=1

    nm=length(xys(1,:));
    for n=1:nm
        p(n)=create_marker(handles.axes1,xys(1,n),xys(2,n),r,[0 1 0],'a',n);
    end

    for n=1:nm
        make_round(p(n),handles.axes1); % use make_round to update circles positions
    end

    update_curve;
end

% --- Executes during object creation, after setting all properties.
function prs_CreateFcn(hObject, eventdata, handles)
% hObject    handle to prs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in p3p.
function p3p_Callback(hObject, eventdata, handles)

global xys f res mth iscnt

plot_3p;


text:
function [sens] = ft_appendsens(cfg, varargin)

% FT_APPENDSENS concatenates multiple sensor definitions that have been processed
% separately.
%
% Use as
%   combined = ft_appendsens(cfg, sens1, sens2, ...)
%
% A call to FT_APPENDSENS results in the label, pos and ori fields to be
% concatenated, and the tra matrix to be merged. Any duplicate electrodes
% will be removed. The labelold and chanposold fields are kept under the
% condition that they are identical across the inputs.
%
% See also FT_ELECTRODEPLACEMENT, FT_ELECTRODEREALIGN, FT_DATAYPE_SENS,
% FT_APPENDDATA, FT_APPENDTIMELOCK, FT_APPENDFREQ, FT_APPENDSOURCE

% Copyright (C) 2017-2018, Arjen Stolk
%
% 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    varargin
ft_preamble provenance varargin
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
% and ensure it is up to the latest standards
for i=1:length(varargin)
  assert(ft_datatype(varargin{i}, 'sens'), 'incorrect input, should be a grad/elec/opto structure');
  varargin{i} = ft_datatype_sens(varargin{i});
end

% do a basic check whether the senstype, units, and coordinate systems match
senstype = cell(1,length(varargin));
for i=1:length(varargin)
  senstype{i} = ft_senstype(varargin{i});
end
typematch = all(strcmp(senstype{1}, senstype));

unitmatch = 1;
if isfield(varargin{1}, 'unit')
  unit = cell(1,length(varargin));
  for i=1:length(varargin)
    unit{i} = varargin{i}.unit;
  end
  unitmatch = all(strcmp(unit{1}, unit));
end

coordsysmatch = 1;
if isfield(varargin{1}, 'coordsys')
  coordsys = cell(1,length(varargin));
  for i=1:length(varargin)
    coordsys{i} = varargin{i}.coordsys;
  end
  coordsysmatch = all(strcmp(coordsys{1}, coordsys));
end

if ~typematch || ~unitmatch || ~coordsysmatch
  ft_error('the senstype, units, or coordinate systems of the inputs do not match');
end

% keep these fields (when present) in the output
sens = keepfields(varargin{1}, {'type', 'unit', 'coordsys'});

% make inventory
haselecpos = 0;
hascoilpos = 0;
hascoilori = 0;
haschanori = 0;
hasoptopos = 0;
hastra     = 0;
haslabelold = 0;
haschanposold = 0;
for i=1:length(varargin)
  % the following fields should be present in any sens structure
  if isfield(varargin{i}, 'label')
    label{i} = varargin{i}.label(:); % ensure column orientation
  end
  if isfield(varargin{i}, 'chanpos')
    chanpos{i} = varargin{i}.chanpos;
  end

  % some the following fields are likely present in a sens structure
  if isfield(varargin{i}, 'elecpos') % EEG
    elecpos{i} = varargin{i}.elecpos;
    haselecpos = 1;
  end
  if isfield(varargin{i}, 'coilpos') % MEG
    coilpos{i} = varargin{i}.coilpos;
    hascoilpos = 1;
  end
  if isfield(varargin{i}, 'coilori') % MEG
    coilori{i} = varargin{i}.coilori;
    hascoilori = 1;
  end
  if isfield(varargin{i}, 'chanori') % MEG
    chanori{i} = varargin{i}.chanori;
    haschanori = 1;
  end
  if isfield(varargin{i}, 'optopos') % NIRS
    optopos{i} = varargin{i}.optopos;
    hasoptopos = 1;
  end
  if isfield(varargin{i}, 'tra') % tra
    tra{i} = varargin{i}.tra;
    hastra = 1;
  end

  % the following fields might be present in a sens structure
  if isfield(varargin{i}, 'labelold')
    labelold{i} = varargin{i}.labelold(:); % ensure column orientation
    haslabelold = 1;
  end
  if isfield(varargin{i}, 'chanposold')
    chanposold{i} = varargin{i}.chanposold;
    haschanposold = 1;
  end
end

% concatenate the main fields and remove duplicates
sens.label = cat(1,label{:});
[dum, labidx] = unique(sens.label, 'stable');
if ~isequal(numel(labidx), numel(sens.label))
  fprintf('removing duplicate labels\n')
  sens.label = sens.label(labidx);
end

sens.chanpos = cat(1,chanpos{:});
[dum, chanidx] = unique(sens.chanpos, 'rows', 'stable');
if ~isequal(numel(chanidx), size(sens.chanpos,1))
  fprintf('removing duplicate channels\n')
  sens.chanpos = sens.chanpos(chanidx,:);
  if ~isequal(labidx, chanidx) % check for matching order
    ft_error('inconsistent order or number of channel labels and positions')
  end
end
if ~isequal(numel(sens.label), size(sens.chanpos,1)) % check for matching number
  ft_error('inconsistent number of channel labels and positions')
end
if haschanori
  sens.chanori = cat(1,chanori{:});
  if ~isequal(numel(chanidx), size(sens.chanpos,1))
    sens.chanori = sens.chanori(chanidx,:); % chanori should match chanpos
  end
end

if hastra && ~any(cellfun(@isempty, tra))
  sens.tra = [];
  for t = 1:numel(tra)
    trarow = [1:size(tra{t},1)]+size(sens.tra,1);
    tracol = [1:size(tra{t},2)]+size(sens.tra,2);
    sens.tra(trarow, tracol) = tra{t};
  end
end

if haselecpos
  sens.elecpos = cat(1,elecpos{:});
  [dum, elecidx, elecidx2] = unique(sens.elecpos, 'rows', 'stable');
  if ~isequal(numel(elecidx), size(sens.elecpos,1))
    fprintf('removing duplicate electrodes\n')
    sens.elecpos = sens.elecpos(elecidx,:);
  end
  if isfield(sens, 'tra')
    % shape duplicates into a single column, if necessary
    for idx = 1:numel(elecidx)
      tmp(:, idx) = sum(sens.tra(chanidx, find(idx==elecidx2)),2); % find and take sum over duplicates
    end
    sens.tra = tmp;
    % check for expected size and non-empty rows
    if ~isequal(size(sens.tra,1), size(sens.chanpos,1)) || ~isequal(size(sens.tra,2), size(sens.elecpos,1)) ...
        || any(any(sens.tra,2)==0) % all channels need to be linked to their origins
      fprintf('removing inconsistent tra matrix\n')
      sens = rmfield(sens, 'tra');
    end
  end
end

if hasoptopos
  sens.optopos = cat(1,optopos{:});
  [dum, optoidx, optoidx2] = unique(sens.optopos, 'rows', 'stable');
  if ~isequal(numel(optoidx), size(sens.optopos,1))
    fprintf('removing duplicate optodes\n')
    sens.optopos = sens.optopos(optoidx,:);
  end
  if isfield(sens, 'tra')
    % shape duplicates into a single column, if necessary
    for idx = 1:numel(optoidx)
      tmp(:, idx) = sum(sens.tra(chanidx, find(idx==optoidx2)),2); % find and take sum over duplicates
    end
    sens.tra = tmp;
    % check for expected size and non-empty rows
    if ~isequal(size(sens.tra,1), size(sens.chanpos,1)) || ~isequal(size(sens.tra,2), size(sens.optopos,1)) ...
        || any(any(sens.tra,2)==0) % all channels need to be linked to their origins
      fprintf('removing inconsistent tra matrix\n')
      sens = rmfield(sens, 'tra');
    end
  end
end

if hascoilpos
  sens.coilpos = cat(1,coilpos{:});
  [dum, coilidx, coilidx2] = unique(sens.coilpos, 'rows', 'stable');
  if ~isequal(numel(coilidx), size(sens.coilpos,1))
    fprintf('removing duplicate coils\n')
    sens.coilpos = sens.coilpos(coilidx,:);
  end
  if isfield(sens, 'tra')
    % shape duplicates into a single column, if necessary
    for idx = 1:numel(coilidx)
      tmp(:, idx) = sum(sens.tra(chanidx, find(idx==coilidx2)),2); % find and take sum over duplicates
    end
    sens.tra = tmp;
    % check for expected size and non-empty rows
    if ~isequal(size(sens.tra,1), size(sens.chanpos,1)) || ~isequal(size(sens.tra,2), size(sens.coilpos,1)) ...
        || any(any(sens.tra,2)==0) % all channels need to be linked to their origins
      fprintf('removing inconsistent tra matrix\n')
      sens = rmfield(sens, 'tra');
    end
  end
end
if hascoilori
  sens.coilori = cat(1,coilori{:});
  if ~isequal(numel(coilidx), size(sens.coilori,1))
    sens.coilori = sens.coilori(coilidx,:); % coilori should match coilpos
  end
end

% keep the following fields only when identical across inputs
if haslabelold && all(isequal(labelold{1}, labelold{:})) % labeloldmatch
  sens.labelold = labelold{1};
end
if haschanposold && all(isequal(chanposold{1}, chanposold{:})) % chanposoldmatch
  sens.chanposold = chanposold{1};
end

% ensure up-to-date and consistent output sensor description
sens = ft_datatype_sens(sens);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous varargin
ft_postamble provenance sens
ft_postamble history sens
ft_postamble savevar sens

author: spm


text:
%% divisors
% Below is a demonstration of the features of the |divisors| function

%%
clear; close all; clc;

%% Syntax
% |d=divisors(x,arg);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function [mbar] = dynpcm22mbar(dynpcm2)
% Convert pressure from dynes per square centimeter to millibar
% Chad Greene 2012
mbar = dynpcm2*0.00100000;


text:
function [OffDec,OffVel] = LCSA_CoefficientSMPSOOperator(Particle,Pbest,Gbest,xlower,xupper)
% -----------------------------------------------------------------------
%  Copyright (C) 2020 Heiner Zille
%
%  This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%  International License. (CC BY-NC-SA 4.0). To view a copy of this license,
%  visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or see the
%  pdf-file "License-CC-BY-NC-SA-4.0.pdf" that came with this code.
%
%  You are free to:
%  * Share ? copy and redistribute the material in any medium or format
%  * Adapt ? remix, transform, and build upon the material
%  Under the following terms:
%  * Attribution ? You must give appropriate credit, provide a link to the
%     license, and indicate if changes were made. You may do so in any reasonable
%     manner, but not in any way that suggests the licensor endorses you or your use.
%  * NonCommercial ? You may not use the material for commercial purposes.
%  * ShareAlike ? If you remix, transform, or build upon the material, you must
%    distribute your contributions under the same license as the original.
%  * No additional restrictions ? You may not apply legal terms or technological
%    measures that legally restrict others from doing anything the license permits.
%
%  Author of this Code:
%   Heiner Zille <heiner.zille@ovgu.de> or <heiner.zille@gmail.com>
%
%  This code is based on the following publications:
%
%  1) Heiner Zille
%     "Large-scale Multi-objective Optimisation: New Approaches and a Classification of the State-of-the-Art"
%     PhD Thesis, Otto von Guericke University Magdeburg, 2019
%     http://dx.doi.org/10.25673/32063
%
%  2) Heiner Zille and Sanaz Mostaghim
%     "Linear Search Mechanism for Multi- and Many-Objective Optimisation"
%     10th International Conference on Evolutionary Multi-Criterion Optimization (EMO 2019),
%        Lecture Notes in Computer Science, vol 11411.
%        Deb K. et al. (eds), Springer, Cham, East Lansing, Michigan, USA, March 2019
%     https://doi.org/10.1007/978-3-030-12598-1_32.
%
%  This file is intended to work with the PlatEMO framework version 2.5.
%  Date of publication of this code: 06.04.2020
%  Last Update of this code: 06.04.2020
%  A newer version of this algorithm may be available. Please contact the author
%  or see http://www.ci.ovgu.de/Research/Codes.html.
%
% The files may have been modified in Feb 2021 by the authors of the Platemo framework to work with the Platemo 3.0 release.
% -----------------------------------------------------------------------
% This file is derived from its original version containied in the PlatEMO
% framework.
% -----------------------------------------------------------------------

    %% Parameter setting
    [ParticleDec,ParticleVel] = unpackDecAndVel(Particle.adds);
    [PbestDec,~] = unpackDecAndVel(Pbest.adds);
    [GbestDec,~] = unpackDecAndVel(Gbest.adds);
    [N,D]        = size(ParticleDec);

    %% Particle swarm optimization
    W  = repmat(unifrnd(0.1,0.5,N,1),1,D);
    r1 = repmat(rand(N,1),1,D);
    r2 = repmat(rand(N,1),1,D);
    C1 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    C2 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    OffVel = W.*ParticleVel + C1.*r1.*(PbestDec-ParticleDec) + C2.*r2.*(GbestDec-ParticleDec);
    phi    = max(4,C1+C2);
    OffVel = OffVel.*2./abs(2-phi-sqrt(phi.^2-4*phi));
    delta  = repmat((xupper-xlower)/2,N,1);
    OffVel = max(min(OffVel,delta),-delta);
    OffDec = ParticleDec + OffVel;

    %% Deterministic back
    Lower  = repmat(xlower,N,1);
    Upper  = repmat(xupper,N,1);
    repair = OffDec < Lower | OffDec > Upper;
    OffVel(repair) = 0.001*OffVel(repair);
    OffDec = max(min(OffDec,Upper),Lower);

    %% Polynomial mutation
    disM  = 20;
    Site1 = repmat(rand(N,1)<0.15,1,D);
    Site2 = rand(N,D) < 1/D;
    mu    = rand(N,D);
    temp  = Site1 & Site2 & mu<=0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                   (1-(OffDec(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp  = Site1 & Site2 & mu>0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                   (1-(Upper(temp)-OffDec(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
end

function [positions, velocities] = unpackDecAndVel(pack)
   noOfSolutions = size(pack,1);
   P = arrayfun(@(K) pack(K).xDecs, 1:noOfSolutions, 'UniformOutput',0);
   positions = cell2mat(transpose(P));
   V = arrayfun(@(K) pack(K).xDecs, 1:noOfSolutions, 'UniformOutput',0);
   velocities = cell2mat(transpose(V));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function table_io_test01 ( )

%*****************************************************************************80
%
%% TABLE_IO_TEST01 tests R8MAT_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2010
%
%  Author:
%
%    John Burkardt
%
  n = 20;
  m = 5;
  output_filename = 'r8mat_05_00020.txt';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  R8MAT_WRITE writes an R8MAT file.\n' );

  table = r8mat_indicator ( m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );
  fprintf ( 1, '  Number of points N  = %d\n', n );

  r8mat_print_some ( m, n, table, 1, 1, 5, 5, ...
    '  5x5 portion of the data written to file:' );

  r8mat_transpose_print_some ( m, n, table, 1, 1, 5, 5, ...
    '  5x5 portion of the TRANSPOSED data:' );

  r8mat_write ( output_filename, m, n, table );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote the file "%s".\n', output_filename );

  return
end


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Globally Optimal Solutions,
%           Solvability/Observability Analysis and Uncertainty Description.
%           IEEE Transactions on Robotics.
%           https://doi.org/10.1109/TRO.2022.3155880
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com
%
%
% syms_hand_eye.m: Generating symbolic expressions for solving hand-eye
%                  calibration problem

clear all
close all
clc

addpath('func_files');
addpath('utils');

syms A11 A12 A13 A14 real
syms A21 A22 A23 A24 real
syms A31 A32 A33 A34 real

A = [
    A11, A12, A13, A14;
    A21, A22, A23, A24;
    A31, A32, A33, A34;
      0,   0,   0,   1;
      ];

syms B11 B12 B13 B14 real
syms B21 B22 B23 B24 real
syms B31 B32 B33 B34 real

B = [
    B11, B12, B13, B14;
    B21, B22, B23, B24;
    B31, B32, B33, B34;
      0,   0,   0,   1;
      ];

syms q0 q1 q2 q3 real
syms t1 t2 t3 lambda real
q = [q0; q1; q2; q3];
tt = [t1; t2; t3];
RR = q2R(q);
XX = [
    RR, tt;
    zeros(1, 3), 1];
J_pure = J_func_hand_eye(XX, A, B);
[coef_J_pure, mon_J_pure] = coeffs(J_pure, [q; tt]);
generateFuncFile(coef_J_pure, fullfile('func_files', 'coef_J_pure_hand_eye_func_new.m'), {A(1 : 3, :), B(1 : 3, :)});
generateFuncFile(mon_J_pure, fullfile('func_files', 'mon_J_pure_hand_eye_func_new.m'), {q, tt});
J = J_pure - 0.5 * lambda * (q.' * q - 1);
x = [q; tt; lambda];
Jacob = jacobian(J, x);
Jacob = expand(Jacob.');

coef_len = length(coef_J_pure);
str = sprintf('coef_J = sym(''coef_J'', [1, %d]);', coef_len);
eval(str);
showSyms(symvar(coef_J));
J = coef_J * mon_J_pure.' - 0.5 * lambda * (q.' * q - 1);
x = [q; tt; lambda];
Jacob = jacobian(J, x);
Jacob = expand(Jacob.');

[coeft1, mont1] = coeffs(Jacob(5), tt);
g1 = coeft1(1 : 3).';
[coeftq1, montq1] = coeffs(expand(Jacob(5) - g1.' * tt ), q);
generateFuncFile(coeftq1, fullfile('func_files', 'coeftq1_hand_eye_func_new.m'), {coef_J});
[coeft2, mont2] = coeffs(Jacob(6), tt);
g2 = coeft2(1 : 3).';
[coeftq2, montq2] = coeffs(expand(Jacob(6) - g2.' * tt ), q);
generateFuncFile(coeftq2, fullfile('func_files', 'coeftq2_hand_eye_func_new.m'), {coef_J});
[coeft3, mont3] = coeffs(Jacob(7), tt);
g3 = coeft3(1 : 3).';
[coeftq3, montq3] = coeffs(expand(Jacob(7) - g3.' * tt ), q);
generateFuncFile(coeftq3, fullfile('func_files', 'coeftq3_hand_eye_func_new.m'), {coef_J});
G = - [g1.'; g2.'; g3.'];
generateFuncFile(G, fullfile('func_files', 'G_hand_eye_func_new.m'), {coef_J});

pinvG = sym('pinvG', [3, 3]);
showSyms(symvar(pinvG));
coefs_tq = sym('coefs_tq', [3, 11]);
showSyms(symvar(coefs_tq));
ts = pinvG * coefs_tq * montq1.';
t1 = ts(1);
t2 = ts(2);
t3 = ts(3);
generateFuncFile(t1, fullfile('func_files', 't1_hand_eye_func_new.m'), {pinvG, coefs_tq, q});
generateFuncFile(t2, fullfile('func_files', 't2_hand_eye_func_new.m'), {pinvG, coefs_tq, q});
generateFuncFile(t3, fullfile('func_files', 't3_hand_eye_func_new.m'), {pinvG, coefs_tq, q});

[coef_Jacob1_qt, mon_Jacob1_qt] = coeffs(Jacob(1) + lambda * q0, [q; tt]);
generateFuncFile(coef_Jacob1_qt, fullfile('func_files', 'coef_Jacob1_qt_hand_eye_func_new.m'), {coef_J});
[coef_Jacob2_qt, mon_Jacob2_qt] = coeffs(Jacob(2) + lambda * q1, [q; tt]);
generateFuncFile(coef_Jacob2_qt, fullfile('func_files', 'coef_Jacob2_qt_hand_eye_func_new.m'), {coef_J});
[coef_Jacob3_qt, mon_Jacob3_qt] = coeffs(Jacob(3) + lambda * q2, [q; tt]);
generateFuncFile(coef_Jacob3_qt, fullfile('func_files', 'coef_Jacob3_qt_hand_eye_func_new.m'), {coef_J});
[coef_Jacob4_qt, mon_Jacob4_qt] = coeffs(Jacob(4) + lambda * q3, [q; tt]);
generateFuncFile(coef_Jacob4_qt, fullfile('func_files', 'coef_Jacob4_qt_hand_eye_func_new.m'), {coef_J});

coef_Jacob1_qt_syms = sym('coef_Jacob1_qt_syms', [1, 36]);
coef_Jacob2_qt_syms = sym('coef_Jacob2_qt_syms', [1, 36]);
coef_Jacob3_qt_syms = sym('coef_Jacob3_qt_syms', [1, 36]);
coef_Jacob4_qt_syms = sym('coef_Jacob4_qt_syms', [1, 36]);
showSyms(symvar(coef_Jacob1_qt_syms));
showSyms(symvar(coef_Jacob2_qt_syms));
showSyms(symvar(coef_Jacob3_qt_syms));
showSyms(symvar(coef_Jacob4_qt_syms));
coef_Jacob_qt_syms = [
    coef_Jacob1_qt_syms;
    coef_Jacob2_qt_syms;
    coef_Jacob3_qt_syms;
    coef_Jacob4_qt_syms;
    ];
Jacob_ = coef_Jacob_qt_syms * mon_Jacob1_qt.';

eqs = expand(eval([Jacob_; Jacob(8)]));

f0 = eqs(1);
f1 = eqs(2);
f2 = eqs(3);
f3 = eqs(4);
[coef_f0_q, mon_f0_q] = coeffs(f0, q);
[coef_f1_q, mon_f1_q] = coeffs(f1, q);
[coef_f2_q, mon_f2_q] = coeffs(f2, q);
[coef_f3_q, mon_f3_q] = coeffs(f3, q);
coef_f0_q_sym = sym('coef_f0_q_sym', [1, 24]);
coef_f1_q_sym = sym('coef_f1_q_sym', [1, 24]);
coef_f2_q_sym = sym('coef_f2_q_sym', [1, 24]);
coef_f3_q_sym = sym('coef_f3_q_sym', [1, 24]);
showSyms(symvar(coef_f0_q_sym));
showSyms(symvar(coef_f1_q_sym));
showSyms(symvar(coef_f2_q_sym));
showSyms(symvar(coef_f3_q_sym));

coef_f_q_sym = [
    coef_f0_q;
    coef_f1_q;
    coef_f2_q;
    coef_f3_q;
    ];
generateFuncFile(coef_f_q_sym, fullfile('func_files', 'coef_f_q_sym_hand_eye_func_new.m'), {pinvG, coefs_tq, coef_Jacob_qt_syms});

f0 = coef_f0_q_sym * mon_f0_q.';
f1 = coef_f1_q_sym * mon_f1_q.';
f2 = coef_f2_q_sym * mon_f2_q.';
f3 = coef_f3_q_sym * mon_f3_q.';
eq = expand([
     f0 * q1 - f1 * q0;
     f0 * q2 - f2 * q0;
     f0 * q3 - f3 * q0;
     q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3 - 1;
     ]);
generateFuncFile(eq, fullfile('func_files', 'eq_hand_eye_func_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym, q});

eq_ = eq;
eq_ = expand(eval(subs(eq_, q0^2, (1 - q1^2 - q2^2 - q3^2))));
[coefeq1, moneq1] = coeffs(eq_(1), q);
mon2th1 = [ moneq1(10), moneq1(14), moneq1(16), moneq1(23), moneq1(27), moneq1(29), moneq1(33), moneq1(35), moneq1(37)];
mon4th1 = [ moneq1(1 : 9), moneq1(11 : 13), moneq1(15), moneq1(17 : 22), moneq1(24 : 26), moneq1(28), moneq1(30 : 32), moneq1(34), moneq1(36), moneq1(38)];

[coefeq2, moneq2] = coeffs(eq_(2), q);
mon2th2 = [ moneq2(10), moneq2(14), moneq2(16), moneq2(23), moneq2(27), moneq2(29), moneq2(33), moneq2(35), moneq2(37)];
mon4th2 = [ moneq2(1 : 9), moneq2(11 : 13), moneq2(15), moneq2(17 : 22), moneq2(24 : 26), moneq2(28), moneq2(30 : 32), moneq2(34), moneq2(36), moneq2(38)];

[coefeq3, moneq3] = coeffs(eq_(3), q);
mon2th3 = [ moneq3(10), moneq3(14), moneq3(16), moneq3(23), moneq3(27), moneq3(29), moneq3(33), moneq3(35), moneq3(37)];
mon4th3 = [ moneq3(1 : 9), moneq3(11 : 13), moneq3(15), moneq3(17 : 22), moneq3(24 : 26), moneq3(28), moneq3(30 : 32), moneq3(34), moneq3(36), moneq3(38)];

coef2th1 = [ coefeq1(10), coefeq1(14), coefeq1(16), coefeq1(23), coefeq1(27), coefeq1(29), coefeq1(33), coefeq1(35), coefeq1(37)];
coef4th1 = [ coefeq1(1 : 9), coefeq1(11 : 13), coefeq1(15), coefeq1(17 : 22), coefeq1(24 : 26), coefeq1(28), coefeq1(30 : 32), coefeq1(34), coefeq1(36), coefeq1(38)];

coef2th2 = [ coefeq2(10), coefeq2(14), coefeq2(16), coefeq2(23), coefeq2(27), coefeq2(29), coefeq2(33), coefeq2(35), coefeq2(37)];
coef4th2 = [ coefeq2(1 : 9), coefeq2(11 : 13), coefeq2(15), coefeq2(17 : 22), coefeq2(24 : 26), coefeq2(28), coefeq2(30 : 32), coefeq2(34), coefeq2(36), coefeq2(38)];

coef2th3 = [ coefeq3(10), coefeq3(14), coefeq3(16), coefeq3(23), coefeq3(27), coefeq3(29), coefeq3(33), coefeq3(35), coefeq3(37)];
coef4th3 = [ coefeq3(1 : 9), coefeq3(11 : 13), coefeq3(15), coefeq3(17 : 22), coefeq3(24 : 26), coefeq3(28), coefeq3(30 : 32), coefeq3(34), coefeq3(36), coefeq3(38)];

y = mon2th1.';
generateFuncFile(y, fullfile('func_files', 'y_func_hand_eye_new.m'), {q});
v = mon4th1(1 : end - 1).';
generateFuncFile(v, fullfile('func_files', 'v_func_hand_eye_new.m'), {q});

D = [
    coef4th1(1 : end - 1);
    coef4th2(1 : end - 1);
    coef4th3(1 : end - 1);
    ];

G = [
    coef2th1;
    coef2th2;
    coef2th3;
    ];

c = [
    coef4th1(end);
    coef4th2(end);
    coef4th3(end);
    ];

generateFuncFile(D, fullfile('func_files', 'D_func_hand_eye_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym});
generateFuncFile(G, fullfile('func_files', 'G_func_hand_eye_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym});
generateFuncFile(c, fullfile('func_files', 'c_func_hand_eye_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym});

HH = jacobian(eq, q);
generateFuncFile(HH, fullfile('func_files', 'Jacob_hand_eye_func_new.m'), {coef_f0_q_sym, coef_f1_q_sym, coef_f2_q_sym, coef_f3_q_sym, q});
gradient = expand(HH.' * eq);

[coef1, mon1] = coeffs(eqs(1), x);
[coef2, mon2] = coeffs(eqs(2), x);
[coef3, mon3] = coeffs(eqs(3), x);
[coef4, mon4] = coeffs(eqs(4), x);

Q_sym = [
    coef1(11), coef1(18), coef1(22), coef1(24);
    coef2(11), coef2(18), coef2(22), coef2(24);
    coef3(11), coef3(18), coef3(22), coef3(24);
    coef4(11), coef4(18), coef4(22), coef4(24);
    ];
generateFuncFile(Q_sym, fullfile('func_files', 'Q_hand_eye_func_new.m'), {pinvG, coefs_tq, coef_Jacob_qt_syms});

res = expand(eqs(1 : 4) - Q_sym * q);
H = jacobian(res, q);
h1 = H(:, 1);
h2 = H(:, 2);
h3 = H(:, 3);
h4 = H(:, 4);
P1 = jacobian(h1, q);
P2 = jacobian(h2, q);
P3 = jacobian(h3, q);
P4 = jacobian(h4, q);

for i = 1 : 4
    for j = 1 : 4
        str = sprintf('W%d%d = [jacobian(P%d(1, :), q%d).''; jacobian(P%d(2, :), q%d).''; jacobian(P%d(3, :), q%d).''; jacobian(P%d(4, :), q%d).''];', ...
                        i, j, i, j - 1, i, j - 1, i, j - 1, i, j - 1);
        eval(str);
    end
end

W1 = [W11, W12, W13, W14];
W2 = [W21, W22, W23, W24];
W3 = [W31, W32, W33, W34];
W4 = [W41, W42, W43, W44];
v = kron(q, q);
u = kron(v, q);
W = - [W1, W2, W3, W4] / 6;
generateFuncFile(W, fullfile('func_files', 'W_hand_eye_func_new.m'), {pinvG, coefs_tq, coef_Jacob_qt_syms});

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function f=comp_iwpfbt(c,wtNodes,pOutIdxs,chOutIdxs,Ls,ext,interscaling)
%COMP_IWFBT Compute Inverse Wavelet Packet Filter-Bank Tree
%   Usage:  f=comp_iwpfbt(c,wtNodes,pOutIdxs,chOutIdxs,Ls,ext)
%
%   Input parameters:
%         c          : Coefficients stored in cell array.
%         wtNodes    : Filterbank tree nodes (elementary filterbans) in
%                      reverse BF order. Cell array of structures of length *nodeNo*.
%         pOutIdxs   : Idx of each node's parent. Array of length *nodeNo*.
%         chOutIdxs  : Idxs of each node children. Cell array of vectors of
%                      length *nodeNo*.
%         ext        : Type of the forward transform boundary handling.
%
%   Output parameters:
%         f          : Reconstructed data in L*W array.
%

% Do non-expansve transform if ext=='per'
doPer = strcmp(ext,'per');

interscalingfac = 1;
if strcmp('intscale',interscaling)
    interscalingfac = 1/2;
elseif strcmp('intsqrt',interscaling)
    interscalingfac = 1/sqrt(2);
end

% For each node in tree in the BF order...
 for jj=1:length(wtNodes)
    % Node filters to a cell array
    %gCell = cellfun(@(gEl)conj(flipud(gEl.h(:))),wtNodes{jj}.g(:),'UniformOutput',0);
    gCell = cellfun(@(gEl)gEl.h(:),wtNodes{jj}.g(:),'UniformOutput',0);
    % Node filters subs. factors
    a = wtNodes{jj}.a;
    % Node filters initial skips
    if(doPer)
       %offset = cellfun(@(gEl) 1-numel(gEl.h)-gEl.offset,wtNodes{jj}.g);
       offset = cellfun(@(gEl) gEl.offset,wtNodes{jj}.g);
    else
       offset = -(a-1);
    end

    if(pOutIdxs(jj))
       % Run filterbank and add to the existing subband.
       ctmp = comp_ifilterbank_td(c(chOutIdxs{jj}),gCell,a,size(c{pOutIdxs(jj)},1),offset,ext);
       c{pOutIdxs(jj)} = c{pOutIdxs(jj)}+ctmp;

       c{pOutIdxs(jj)} = interscalingfac*c{pOutIdxs(jj)};

    else
       % We are at the root.
       f = comp_ifilterbank_td(c(chOutIdxs{jj}),gCell,a,Ls,offset,ext);
    end
 end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

author: Allopart

repo: rbpf-gmapping


text:
% clear all; close all; clc;

img1=imread('lena.jpg');
imshow(img1,[]);
[m n]=size(img1);

img=sqrt(double(img1));      %伽马校正
% figure;hold on;
figure;
a = min(min(min(img)));
b = max(max(max(img)));
% can not directly operator on a matrix more than three dimensions
c = 255/(b-a);
for i = 1:3
    img2(:,:,i) =img(:,:,i) * c;
end

imshow(uint8(img2));

%下面是求边缘
fy=[-1 0 1];        %定义竖直模板,用于计算y方向上的梯度值，向右为正方向
fx=fy';             %定义水平模板,用于计算x方向上的梯度值
Iy=imfilter(img,fy,'replicate');    %竖直边缘/梯度
Ix=imfilter(img,fx,'replicate');    %水平边缘/梯度
Ied=sqrt(Ix.^2+Iy.^2);              %边缘强度/梯度
Iphase=Iy./Ix;              %边缘斜率，有些为inf,-inf,nan，其中nan需要再处理一下

%下面是求cell
step=16;                %step*step个像素作为一个cell
orient=9;               %方向直方图的方向个数
jiao=360/orient;        %每个方向包含的角度数
Cell=cell(1,1);              %所有的角度直方图,cell是可以动态增加的，所以先设了一个
ii=1;
jj=1;
for i=1:step:m          %如果处理的m/step不是整数，最好是i=1:step:m-step， no overlapping
    ii=1;
    for j=1:step:n      %注释同上
        tmpx=Ix(i:i+step-1,j:j+step-1);
        tmped=Ied(i:i+step-1,j:j+step-1);
        tmped=tmped/sum(sum(tmped));        %局部边缘强度归一化
        tmpphase=Iphase(i:i+step-1,j:j+step-1);
        Hist=zeros(1,orient);               %当前step*step像素块统计角度直方图,就是cell
        for p=1:step
            for q=1:step
                if isnan(tmpphase(p,q))==1  %0/0会得到nan，如果像素是nan，重设为0
                    tmpphase(p,q)=0;
                end
                ang=atan(tmpphase(p,q));    %atan求的是[-90 90]度之间
                ang=mod(ang*180/pi,360);    %全部变正，[0,360), -90变270
                if tmpx(p,q)<0              %根据x方向确定真正的角度
                    if ang<90               %如果是第一象限
                        ang=ang+180;        %移到第三象限
                    end
                    if ang>270              %如果是第四象限
                        ang=ang-180;        %移到第二象限
                    end
                end
                ang=ang+0.0000001;          %防止ang为0
                Hist(ceil(ang/jiao))=Hist(ceil(ang/jiao))+tmped(p,q);   %ceil向上取整，使用边缘强度加权
            end
        end
        Hist=Hist/sum(Hist);    %方向直方图归一化
        Cell{ii,jj}=Hist;       %放入Cell中
        ii=ii+1;                %针对Cell的y坐标循环变量
    end
    jj=jj+1;                    %针对Cell的x坐标循环变量
end

%下面是求feature,2*2个cell合成一个block,没有显式的求block
[m n]=size(Cell);
feature=cell(1,(m-1)*(n-1));%step = size(Cell,1),而且每个block由2*2的cell组成
for i=1:m-1
   for j=1:n-1
        f=[];
        f=[f Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+1,j}(:)' Cell{i+1,j+1}(:)'];
        feature{(i-1)*(n-1)+j}=f;
   end
end

%到此结束，feature即为所求
%下面是为了显示而写的
l=length(feature);
f=[];
for i=1:l
    f=[f;feature{i}(:)'];
end
figure
mesh(f)

author: zhangqianqianQQ


text:
function demoNcutImage;
% demoNcutImage
%
% demo for NcutImage
% also initialize matlab paths to subfolders
% Timothee Cour, Stella Yu, Jianbo Shi, 2004.

disp('Ncut Image Segmentation demo');

%% read image, change color image to brightness image, resize to 160x160
I = imread_ncut('jpg_images/3.jpg',160,160);

%% display the image
figure(1);clf; imagesc(I);colormap(gray);axis off;
disp('This is the input image to segment, press Enter to continue...');
pause;

%% compute the edges imageEdges, the similarity matrix W based on
%% Intervening Contours, the Ncut eigenvectors and discrete segmentation
nbSegments = 5;
disp('computing Ncut eigenvectors ...');
tic;
[SegLabel,NcutDiscrete,NcutEigenvectors,NcutEigenvalues,W,imageEdges]= NcutImage(I,nbSegments);
disp(['The computation took ' num2str(toc) ' seconds on the ' num2str(size(I,1)) 'x' num2str(size(I,2)) ' image']);

%% display the edges
figure(2);clf; imagesc(imageEdges); axis off
disp('This is the edges computed, press Enter to continue...');
pause;

%% display the segmentation
figure(3);clf
bw = edge(SegLabel,0.01);
J1=showmask(I,imdilate(bw,ones(2,2))); imagesc(J1);axis off
disp('This is the segmentation, press Enter to continue...');
pause;

%% display Ncut eigenvectors
figure(4);clf;set(gcf,'Position',[100,500,200*(nbSegments+1),200]);
[nr,nc,nb] = size(I);
for i=1:nbSegments
    subplot(1,nbSegments,i);
    imagesc(reshape(NcutEigenvectors(:,i) , nr,nc));axis('image');axis off;
end
disp('This is the Ncut eigenvectors...');
disp('The demo is finished.');

author: jwyang


text:
function parrec = spm_cfg_parrec
% SPM Configuration file for Philips PAR/REC Import
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_cfg_parrec.m 6315 2015-01-23 17:09:06Z guillaume $

%--------------------------------------------------------------------------
% data PAR files
%--------------------------------------------------------------------------
data         = cfg_files;
data.tag     = 'data';
data.name    = 'PAR files';
data.help    = {'Select the PAR files to convert.'};
data.filter  = 'any';
data.ufilter = '.*(par|PAR)';
data.num     = [1 Inf];

%--------------------------------------------------------------------------
% outdir Output directory
%--------------------------------------------------------------------------
outdir         = cfg_files;
outdir.tag     = 'outdir';
outdir.name    = 'Output directory';
outdir.help    = {'Select a directory where files are written.'};
outdir.filter  = 'dir';
outdir.ufilter = '.*';
outdir.num     = [1 1];

%--------------------------------------------------------------------------
% ext Output image format
%--------------------------------------------------------------------------
ext         = cfg_menu;
ext.tag     = 'ext';
ext.name    = 'Output image format';
ext.help    = {'Output files can be written as .img + .hdr, or the two can be combined into a .nii file.'};
ext.labels = {
              'Two file (img+hdr) NIfTI'
              'Single file (nii) NIfTI'
}';
ext.values = {'img', 'nii'};
ext.def    = @(val)spm_get_defaults('images.format', val{:});

%--------------------------------------------------------------------------
% opts Options
%--------------------------------------------------------------------------
opts         = cfg_branch;
opts.tag     = 'opts';
opts.name    = 'Options';
opts.val     = {outdir ext};
opts.help    = {'Conversion options'};

%--------------------------------------------------------------------------
% parrec ECAT Import
%--------------------------------------------------------------------------
parrec         = cfg_exbranch;
parrec.tag     = 'parrec';
parrec.name    = 'PAR/REC Import';
parrec.val     = {data opts};
parrec.help    = {'Philips PAR/REC Import.'};
parrec.prog    = @convert_parrec;
parrec.vout    = @vout;

%==========================================================================
function out = convert_parrec(job)
for i=1:numel(job.data)
    spm_parrec2nifti(job.data{i},job.opts);
end

out.files = cell(size(job.data));
for i=1:numel(job.data)
    out.files{i} = spm_file(job.data{i}, 'path',pwd, 'ext',job.opts.ext);
end

%==========================================================================
function dep = vout(job)
dep            = cfg_dep;
dep.sname      = 'Converted Images';
dep.src_output = substruct('.','files');
dep.tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});