text:
function [LapCnt, TrajCnt, s_loc_m, s_glob_m, x_m, y_m, psi_rad, kappa_radpm, v_mps, ax_mps2, ...
    banking_rad, ax_lim_mps2, ay_lim_mps2] =...
  learningTrack(sendMessages, ActualPos, InitialPos, ...
  P_VDC_AccLimits_v_mps, P_VDC_AccLimits_ax_mps2, P_VDC_AccLimits_ay_mps2)
% Authors:     Alexander Wischnewski
%
% Description:
%   creates a local trajectory which covers all relevant lateral accelerations and velocities

nTrajPoints = 50;
nTrajHighResPoints = 300;

% create speed and acceleration array to be tested
TargetArray_v_mps = 10:5:80;
TargetArray_ay_mps2 = 0:2:30;
% maximum curvature to prevent crazy accelerations at low speeds
kappa_lim_radpm = 0.12;
% number of cycles for one operating point
nTestCycles = 150;
% path length specified as foresight in seconds
foresight_s = 2.5;

persistent TrajOld;
persistent currentIdx_v currentIdx_ay OperatingPoint_cnt done backToStraight ayOld vOld

if(isempty(TrajOld))
    TrajOld.LapCnt = uint32(0);
    TrajOld.TrajCnt = uint32(0);
    TrajOld.s_loc_m = zeros(nTrajPoints, 1);
    TrajOld.s_glob_m = zeros(nTrajPoints, 1);
    TrajOld.x_m = zeros(nTrajPoints, 1);
    TrajOld.y_m = zeros(nTrajPoints, 1);
    TrajOld.psi_rad = zeros(nTrajPoints, 1);
    TrajOld.kappa_radpm = zeros(nTrajPoints, 1);
    TrajOld.v_mps = zeros(nTrajPoints, 1);
    TrajOld.ax_mps2 = zeros(nTrajPoints, 1);
    TrajOld.banking_rad = zeros(nTrajPoints, 1);
    TrajOld.ax_lim_mps2 = zeros(nTrajPoints, 1);
    TrajOld.ay_lim_mps2 = zeros(nTrajPoints, 1);
    TrajOld.tube_r_m = zeros(nTrajPoints, 1);
    TrajOld.tube_l_m = zeros(nTrajPoints, 1);
    currentIdx_v = 1;
    currentIdx_ay = 1;
    OperatingPoint_cnt = 0;
    done = 0;
    backToStraight = 0;
    ayOld = 0;
    vOld = 0;
end

if(sendMessages && done == 0)
    % if this is the very first trajectory, start at the initial pose
    if(TrajOld.TrajCnt == 0)
        StartPose = InitialPos;
    else
        % find current trajectory point
        VehicleDynamicState.Pos.x_m = ActualPos(1);
        VehicleDynamicState.Pos.y_m = ActualPos(2);
        VehicleDynamicState.Pos.psi_rad = ActualPos(3);
        [PathPos, ~] = localTrajectoryMatching(VehicleDynamicState, TrajOld);
        % read out trajectory point
        [ActualTrajectoryPoint] = trajectoryInterpolation(TrajOld, PathPos.s_m);
        StartPose = [ActualTrajectoryPoint.x_m, ActualTrajectoryPoint.y_m, ActualTrajectoryPoint.psi_rad];
    end
    % create a path with constant radius and speed in front of the vehicle
    if backToStraight
        % decrease radius smoothly over the time window to prevent spinning
        w_old = (nTestCycles - OperatingPoint_cnt)/nTestCycles;
        TargetAy_mps2 = ayOld*w_old;
        % increase speed smoothly over the time window to prevent spinning
        Target_v_mps_loc = vOld*w_old + TargetArray_v_mps(currentIdx_v)*(1-w_old);
        TargetKappa_radpm = TargetAy_mps2/Target_v_mps_loc^2;
        TargetSpeed_v_mps = Target_v_mps_loc;
    else
        TargetKappa_radpm = TargetArray_ay_mps2(currentIdx_ay)/TargetArray_v_mps(currentIdx_v)^2;
        TargetSpeed_v_mps = TargetArray_v_mps(currentIdx_v);
    end
    % high detail extrapolation of the path
    s_loc_high_res = linspace(0, foresight_s*TargetSpeed_v_mps, nTrajHighResPoints)';
    psi_rad_high_res = normalizeAngle(cumsum([StartPose(3); TargetKappa_radpm*diff(s_loc_high_res)]));
    x_m_high_res = cumsum([StartPose(1); -sin(psi_rad_high_res(1:end-1)).*diff(s_loc_high_res)]);
    y_m_high_res = cumsum([StartPose(2); cos(psi_rad_high_res(1:end-1)).*diff(s_loc_high_res)]);
    % map everything to trajectory structure
    TrajOld.s_loc_m = linspace(0, foresight_s*TargetSpeed_v_mps, nTrajPoints)';
    TrajOld.x_m = x_m_high_res(1:(nTrajHighResPoints/nTrajPoints):end-1);
    TrajOld.y_m = y_m_high_res(1:(nTrajHighResPoints/nTrajPoints):end-1);
    TrajOld.psi_rad = psi_rad_high_res(1:(nTrajHighResPoints/nTrajPoints):end-1);
    TrajOld.kappa_radpm = TargetKappa_radpm*ones(nTrajPoints, 1);
    TrajOld.v_mps = TargetSpeed_v_mps*ones(nTrajPoints, 1);
    % update acceleration limits
    TrajOld.ax_lim_mps2 = interp1(P_VDC_AccLimits_v_mps, P_VDC_AccLimits_ax_mps2, TargetSpeed_v_mps).*ones(nTrajPoints, 1);
    TrajOld.ay_lim_mps2 = interp1(P_VDC_AccLimits_v_mps, P_VDC_AccLimits_ay_mps2, TargetSpeed_v_mps).*ones(nTrajPoints, 1);
    % increase trajectory counter
    TrajOld.TrajCnt = TrajOld.TrajCnt + uint32(1);
    % map this to the outputs
    LapCnt = uint32(1);
    TrajCnt = TrajOld.TrajCnt;
    s_loc_m = TrajOld.s_loc_m;
    s_glob_m = TrajOld.s_glob_m;
    x_m = TrajOld.x_m;
    y_m = TrajOld.y_m;
    psi_rad = TrajOld.psi_rad;
    kappa_radpm = TrajOld.kappa_radpm;
    v_mps = TrajOld.v_mps;
    ax_mps2 = TrajOld.ax_mps2;
    banking_rad = TrajOld.banking_rad;
    ax_lim_mps2 = TrajOld.ax_lim_mps2;
    ay_lim_mps2 = TrajOld.ay_lim_mps2;
    % update operating point
    OperatingPoint_cnt = OperatingPoint_cnt + 1;
    % every six seconds a new operating point is targeted
    if(OperatingPoint_cnt > nTestCycles)
        % reset counter
        OperatingPoint_cnt = 0;
        % check whether we are in normal operation
        if backToStraight == 0
            % go to next lateral acceleration
            currentIdx_ay = currentIdx_ay + 1;
            % if last one was reached go back to straight line driving and increase speed
            % do the same if curvature limit or acceleration limit is exceeded
            if(currentIdx_ay > length(TargetArray_ay_mps2) ...
                    || TargetArray_ay_mps2(currentIdx_ay)/TargetArray_v_mps(currentIdx_v)^2 > kappa_lim_radpm ...
                    || TargetArray_ay_mps2(currentIdx_ay) > interp1(P_VDC_AccLimits_v_mps, P_VDC_AccLimits_ay_mps2, TargetArray_v_mps(currentIdx_v)))
                % store last ay and last speed
                ayOld = TargetArray_ay_mps2(currentIdx_ay - 1);
                vOld = TargetArray_v_mps(currentIdx_v);
                % reset ay to zero
                currentIdx_ay = 1;
                % increase speed
                currentIdx_v = currentIdx_v + 1;
                % go to back to straight mode (smoothly reduce ay to zero)
                backToStraight = 1;
                % check if all operating points have been tested
                if(currentIdx_v > length(TargetArray_v_mps))
                    done = 1;
                end
            end
        else
            % reset to normal operation after one cycle with back to straight
           backToStraight = 0;
        end
    end
else
  % output a zero trajectory as long as the component is not considered active
  LapCnt = uint32(0);
  TrajCnt = uint32(0);
  s_loc_m = zeros(nTrajPoints, 1);
  s_glob_m = zeros(nTrajPoints, 1);
  x_m = zeros(nTrajPoints, 1);
  y_m = zeros(nTrajPoints, 1);
  psi_rad = zeros(nTrajPoints, 1);
  kappa_radpm = zeros(nTrajPoints, 1);
  v_mps = zeros(nTrajPoints, 1);
  ax_mps2 = zeros(nTrajPoints, 1);
  banking_rad = zeros(nTrajPoints, 1);
  ax_lim_mps2 = zeros(nTrajPoints, 1);
  ay_lim_mps2 = zeros(nTrajPoints, 1);
end

author: TUMFTM


text:
function er3OC_sym
%EG3OC_sym    Example 3 of optimal control tutorial.
%    This example is from D.S.Naidu, "Optimal contorl systems"
%    page 77-80, Example 2.14
%    Symbolic toolbox is used to get the analytical solution

sol = dsolve('Dx1 = x2, Dx2 = -p2, Dp1 = 0, Dp2 = -p1',...
           'x1(0) = 1, x2(0) = 2, x1(tf) = 3, p2(tf) = 0');
eq1 = subs(sol.x1) - 'x1tf';
eq2 = subs(sol.x2) - 'x2tf';
eq3 = subs(sol.p1) - 'p1tf';
eq4 = subs(sol.p2) - 'p2tf';
eq5 = sym('p1tf*x2tf - 0.5*p2tf^2');
%%
sol_2 = solve(eq1, eq2, eq3, eq4, eq5);
tf = sol_2.tf;
x1tf = sol_2.x1tf;
x2tf = sol_2.x2tf;

x1 = subs(sol.x1);
x2 = subs(sol.x2);
p1 = subs(sol.p1);
p2 = subs(sol.p2);

%%
sol_book = {@(t)((4/54).*t.^3-(2/3)*t.^2+2.*t+1),...
            @(t)((4/18).*t.^2-(4/3).*t+2)};
u = @(t)((4/9).*t-(4/3));
t = double(tf);
time = linspace(0,t,20);
s_book = [sol_book{1}(time);sol_book{2}(time)];
ut = u(time);

figure(1);
ezplot(x1,[0 t]); hold on;
ezplot(x2,[0 t]);
plot(time, s_book,'*');
plot(time, ut, 'k:');
axis([0 t -1.5 3]);
text(1.3,2.5,'x_1(t)');
text(1.3,1,'x_2(t)');
text(1.3,-.5,'u(t)');
xlabel('time');
ylabel('states');
title('Analytical solution');
hold off;
% print -djpeg90 -r300 eg3a.jpg

%% ------------------------------------------------------------------------
% % The same procedure applied on example problem 1 with tf unknown.
% % no explicit solution could be found
% sol = dsolve('Dx1 = x2, Dx2 = -x2 - p2, Dp1 = 0, Dp2 = p1 - p2',...
%            'x1(0) = 1, x2(0) = 2, x1(tf) = 5, x2(tf) = 2');
% eq1 = subs(sol.x1) - 'x1tf';
% eq2 = subs(sol.x2) - 'x2tf';
% eq3 = subs(sol.p1) - 'p1tf';
% eq4 = subs(sol.p2) - 'p2tf';
% eq5 = sym('p1tf*x2tf-p2tf*x2tf-0.5*p2tf^2');
%
% sol_2 = solve(eq1, eq2, eq3, eq4, eq5);
% tf = sol_2.tf;
% x1tf = sol_2.x1tf;
% x2tf = sol_2.x2tf;
% clear t;
% x1 = subs(sol.x1);
% x2 = subs(sol.x2);
% p1 = subs(sol.p1);
% p2 = subs(sol.p2);

author: Sable

repo: mcbench-benchmarks


text:
function month_name = month_to_month_name_hebrew ( y, m )

%*****************************************************************************80
%
%% MONTH_TO_MONTH_NAME_HEBREW returns the name of a Hebrew month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, the year and month.
%
%    Output, string MONTH_NAME, the month name.
%
  name = [ ...
    'Tishri   '; 'Heshvan  '; 'Kislev   '; 'Tebet    '; 'Shebat   '; ...
    'Adar     '; 'Veadar   '; 'Nisan    '; 'Iyar     '; 'Sivan    '; ...
    'Tammuz   '; 'Ab       '; 'Elul     ' ];

  if ( year_is_embolismic_hebrew ( y ) )

    if ( m < 1 || 13 < m )
      month_name = '?????';
    else
      month_name = name(m,:);
    end

  else

    if ( m < 1 || 12 < m )
      month_name = '?????';
    elseif ( m <= 6 )
      month_name = name(m,:);
    else
      month_name = name(m+1,:);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Objectness measure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Top-level routine implementation
% of the system described in the paper:
%
% B. Alexe, T. Deselaers, V. Ferrari
% What is an object?
% CVPR 2010
%
% The author and copyright holder is Bogdan Alexe.
% You might redistribute, reuse, and mofidy this code for
% research purposes. Any commercial application
% is forbidden without prior agreement in written
% by the author.

function boxes = runObjectness(img,numberProposals,params)
%This function computes the objectness measure and samples boxes from it.

%INPUT
%img - input image
%numberProposals - number of windows sampled from the objectness measure
%params - struct containing parameters of the function (loaded in startup.m)

%OUTPUT
%boxes - samples windows from the objectness measure
%      - each row contains a window using the format [xmin ymin xmax ymax score]

%dir_root = rootpath;%change this to an absolute path

img = gray2rgb(img);
%{
if nargin < 3
    try
        struct = load([dir_root '/Data/params.mat']);
        params = struct.params;
        clear struct;
    catch
        params = defaultParams(dir_root);
        save([dir_root '/Data/params.mat'],'params');
    end
end
%params = updatePath(dir_root,params);
%}
if length(params.cues)==1
    %single cues

    distributionBoxes = computeScores(img,params.cues{1},params);

    switch lower(params.sampling)
        case 'nms'
            %nms sampling

            %consider only params.distribution_windows (= 100k windows)
            if size(distributionBoxes,1) > params.distribution_windows
                indexSamples = scoreSampling(distributionBoxes(:,5),params.distribution_windows,1);
                distributionBoxes = distributionBoxes(indexSamples,:);
            end

            %sampling
            boxes = nms_pascal(distributionBoxes, 0.5,numberProposals);

        case 'multinomial'
            %multinomial sampling

            %sample from the distribution of the scores
            indexSamples = scoreSampling(distributionBoxes(:,end),numberProposals,1);
            boxes = distributionBoxes(indexSamples,:);

        otherwise
            display('sampling procedure unknown')
    end

else
    %combination of cues

    if not(ismember('MS',params.cues))
        display('ERROR: combinations have to include MS');
        boxes = [];
        return
    end

    if length(unique(params.cues)) ~= length(params.cues)
        display('ERROR: repetead cues in the combination');
        boxes = [];
        return
    end

    distributionBoxes = computeScores(img,'MS',params);
    %rearrange the cues such that 'MS' is the first cue
    if ~strcmp(params.cues{1},'MS')
        params.cues{strcmp(params.cues,'MS')} = params.cues{1};
        params.cues{1} ='MS';
    end

    score = zeros(size(distributionBoxes,1),length(params.cues));
    score(:,1) = distributionBoxes(:,end);
    windows = distributionBoxes(:,1:4);
    for idx = 2:length(params.cues)
        temp = computeScores(img,params.cues{idx},params,windows);
        score(:,idx) = temp(:,end);
    end
    scoreBayes = integrateBayes(params.cues,score,params);

    switch lower(params.sampling)
        case 'nms'
            %nms sampling

            distributionBoxes(:,5) = scoreBayes;
            boxes = nms_pascal(distributionBoxes, 0.5, numberProposals);

        case 'multinomial'
            %multinomial sampling

            %sample from the distribution of the scores
            indexSamples = scoreSampling(scoreBayes,numberProposals,1);
        boxes = [windows(indexSamples,:) scoreBayes(indexSamples,:)];

        otherwise
            display('sampling procedure unknown')
    end

end

author: Cloud-CV

repo: object-proposals


text:
%RBSVC Trainable automatic radial basis Support Vector Classifier
%
%   [W,KERNEL,NU,C] = RBSVC(A)
%   [W,KERNEL,NU,C] = A*RBSVC
%
% INPUT
%   A	      Dataset
%
% OUTPUT
%   W       Mapping: Radial Basis Support Vector Classifier
%   KERNEL  Untrained mapping, representing the optimised kernel
%   NU      Resulting value for NU from NUSVC (W = NUSVC(A,KERNEL,C)
%   C       Resulting value for C (W = SVC(A,KERNEL,C)
%
% DESCRIPTION
% This routine computes a classifier by NUSVC using a radial basis kernel
% with an optimised standard deviation by REGOPTC. The resulting classifier
% W is identical to NUSVC(A,KERNEL,NU). As the kernel optimisation is based
% on internal cross-validation the dataset A should be sufficiently large.
% Moreover it is very time-consuming as the kernel optimisation needs
% about 100 calls to SVC.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, PROXM, SVC, NUSVC, REGOPTC

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function [w,kernel,nu,c] = rbsvc(a,sig)

if nargin < 2 | isempty(sig)
	sig = NaN;
end

if nargin < 1 | isempty(a)

	w = prmapping(mfilename,{sig});

else

	islabtype(a,'crisp');
	isvaldfile(a,2,2);              % at least 1 object per class, 2 classes
	a = testdatasize(a,'objects');
  c = getsize(a,3);

	if (c > 2)

    % Compute c classifiers: each class against all others.
		w = mclassc(a,prmapping(mfilename,{sig}));

  else

	  if isnan(sig) % optimise sigma

		  % find upper bound
		  d = sqrt(+distm(a));
		  sigmax = min(max(d)); % max: smallest furthest neighbor distance
      % find lower bound
		  d = d + 1e100*eye(size(a,1));
		  sigmin = max(min(d)); % min: largest nearest neighbor distance
		  % call optimiser
		  defs = {1};
		  parmin_max = [sigmin,sigmax];
		  [w,kernel,nu,c] = regoptc(a,mfilename,{sig},defs,[1],parmin_max,testc([],'soft'));

	  else % kernel is given

		  kernel = proxm([],'r',sig);
		  [w,J,nu,c] = nusvc(a,kernel);

	  end

  end

end

w = setname(w,'RB-SVM');
return

author: marianux

repo: ecg-kit

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


text:
function amino_name = ch_to_amino_name ( c )

%*****************************************************************************80
%
%% CH_TO_AMINO_NAME converts a character to an amino acid name.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl Branden, John Tooze,
%    Introduction to Protein Structure,
%    Garland Publishing, 1991.
%
%  Parameters:
%
%    Input, character C, the one letter code for an amino acid.
%    Lower and upper case letters are treated the same.
%
%    Output, string AMINO_NAME, the full name of the
%    corresponding amino acid.  The longest name is 27 characters.
%    If the input code is not recognized, then AMINO_NAME will be set to '???'.
%
  n = 23;

  amino_table = [
    'Alanine                    ', ...
    'Aspartic acid or Asparagine', ...
    'Cysteine                   ', ...
    'Aspartic acid              ', ...
    'Glutamic acid              ', ...
    'Phenylalanine              ', ...
    'Glycine                    ', ...
    'Histidine                  ', ...
    'Isoleucine                 ', ...
    'Lysine                     ', ...
    'Leucine                    ', ...
    'Methionine                 ', ...
    'Asparagine                 ', ...
    'Proline                    ', ...
    'Glutamine                  ', ...
    'Arginine                   ', ...
    'Serine                     ', ...
    'Threonine                  ', ...
    'Valine                     ', ...
    'Tryptophan                 ', ...
    'Undetermined amino acid    ', ...
    'Tyrosine                   ', ...
    'Glutamic acid or Glutamine ' ];

  ch_table = [...
    'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', ...
    'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', ...
    'X', 'Y', 'Z' ];

  for i = 1 : n
    if ( ch_eqi ( c, ch_table(i) ) )
      amino_name = amino_table(i);
      return
    end
  end

  amino_name = '???';

  return
end


text:
function fx1 = p18_fx1 ( x )

%*****************************************************************************80
%
%% P18_FX1 evaluates the derivative of the function for problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX1(*), the value of the derivative at X.
%
  fx1 = 10.0^14 * ( ...
          7 * x.^6 ...
      -  42 * x.^5 ...
      + 105 * x.^4 ...
      - 140 * x.^3 ...
      + 105 * x.^2 ...
      -  42 * x.^1 ...
      +   7 );

  return
end


text:
% pop_delset() - Delete a dataset from the variable containing
%                all datasets.
%
% Usage: >> ALLEEG = pop_delset(ALLEEG, indices);
%
% Inputs:
%   ALLEEG   - array of EEG datasets
%   indices  - indices of datasets to delete. None -> a pop_up window asks
%              the user to choose. Index < 0 -> it's positive is given as
%              the default in the pop-up window (ex: -3 -> default 3).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: pop_copyset(), eeglab()

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

% load a set and store it in the current set
% ------------------------------------------
function [ALLSET, command] = pop_delset(ALLSET, set_in);

command = '';
if nargin < 1
	help pop_delset;
	return;
end;
if isempty( ALLSET )
	error('Cannot delete dataset. Restart eeglab to clear all dataset information');
    return;
end;

if nargin < 2 | set_in < 0
	% which set to delete
	% -----------------
	promptstr    = { 'Dataset(s) to delete:' };
	if nargin == 2
		inistr       = { int2str(-set_in) };
	else
		inistr       = { '1' };
	end;
	result       = inputdlg2( promptstr, 'Delete dataset -- pop_delset()', 1,  inistr, 'pop_delset');
	size_result  = size( result );
	if size_result(1) == 0 return; end;
	set_in   	 = eval( [ '[' result{1} ']' ] );
end;

if isempty(set_in)
	return;
end;

A = fieldnames( ALLSET );
A(:,2) = cell(size(A));
A = A';
for i = set_in
   try
   		ALLSET(i) = struct(A{:});
		%ALLSET = setfield(ALLSET, {set_in}, A{:}, cell(size(A)));
	catch
		error('Error: no such dataset');
		return;
	end;
end;

command = sprintf('%s = pop_delset( %s, [%s] );', inputname(1), inputname(1), int2str(set_in));
return;

author: PatternRecognition

repo: OpenBMI


text:
function LM2OpenCV(database, HOMEIMAGES, base_folder, num_neg, patch_size)
% LM2OpenCV Outputs query in OpenCV Haar-detector format.
%   LM2OpenCV(database, HOMEIMAGES, base_folder, num_neg, patch_size)
%   converts and stores the object images specified by the LabelMe database
%   struct to the OpenCV haar-detector format. base_folder specifies the
%   target directory for which positive and negative training examples are
%   stored for input to the OpenCV haar-detector training utility
%   (createsamples.exe and haartraining.exe).
%
%   To form positive training examples, object images are copied from the
%   LabelMe database into the folder base_folder/positives and the file
%   positives.txt is created and  stored in base_folder. The positives.txt
%   file lists for each image in database the object bounding boxes and is
%   provided as input to the OpenCV createsamples routine.
%
%   The negative examples are generated by taking num_neg subwindows of
%   size patch_size=[width height] from the images in database that do not
%   contain an object instance. These windows are sampled at random
%   locations in the images. The cropped subwindows are stored in the
%   directory base_folder/negatives and the file negatives.txt is created
%   and stored in base_folder, which lists the image filenames stored in
%   the negatives directory. The negatives.txt file is used by the OpenCV
%   haartraining routine.
%
%   The variable HOMEIMAGES specifies the physical location of the LabelMe
%   database used by the database index database.
%
%   To understand how to train an OpenCV haar-detector with the training
%   samples and files  generated by this script use example command line
%   calls are provided below.
%
%   Example (Pedestrian Detector):
%
%      (In MATLAB)
%
%      HOMEIMAGES = 'C:/LabelMe/Images';
%      HOMEANNOTATIONS = 'C:/LabelMe/Annotations';
%      base_folder = 'D:/dtd';
%      num_neg = 10000;
%      patch_size = [16 32]; % width=16, height=32
%
%      LMdb = LMdatabase(HOMEANNOTATIONS);
%      query = 'pedestrian,person,human,man,woman';
%      database = LMquery(LMdb, 'object.name', query);
%      LM2OpenCV(database, HOMEIMAGES, base_folder, num_neg, patch_size);
%      counts = LMcountobject(database, query);
%
%      (Command Prompt, Let counts=900)
%
%      createsamples -info D:\dtd\positive.txt -vec D:\dtd\positives.vec
%      -num 900 -w 16 -h 32
%
%      haartraining -data D:\dtd\peddetector\ -vec D:\dtd\positives.vec -bg
%      D:\dtd\negatives.txt -npos 900 -nneg 10000 -w 16 -h 32
%
%   In the above example, the call to createsamples generates the file
%   positives.vec used by the haartraining routine. The haartraining
%   routine then creates the OpenCV Haar pedestrian detector and saves it
%   in the director 'D:\dtd\peddetector'.
%

% get number of images in the database
Nimages = length(database);

% sample negatives evenly from each image
nneg_per_image = ceil(num_neg/Nimages);
if(nneg_per_image<1)
    nneg_per_image = 1;
end

% create directories positives and negatives
pos_res = mkdir(base_folder, 'positives');
neg_res = mkdir(base_folder, 'negatives');
if ~(pos_res & neg_res)
    error('LM2OpenCV error: Unable to create positives and/or negatives directories.');
end

% create files positives.txt and negatives.txt
pbdir = sprintf('%s/positives',base_folder);
nbdir = sprintf('%s/negatives',base_folder);
pfilename = sprintf('%s/positives.txt', base_folder);
nfilename = sprintf('%s/negatives.txt', base_folder);
pos_file = fopen(pfilename, 'w+t');
neg_file = fopen(nfilename, 'w+t');
if (pos_file<0 || neg_file<0)
    error('LM2OpenCV error: Unable to create files positive.txt and/or negatives.txt.');
end

% traverse images:
%   - copy over positive image files to the positives directory
%   - crop negative examples and store in negatives directory
%   - create files positives.txt and negatives.txt
neg_ctr = 0;
obj_ctr = 0;
pw = patch_size(1);
ph = patch_size(2);
for i = 1:Nimages
    if isfield(database(i).annotation, 'object')
        Nobjects = length(database(i).annotation.object);
        try
            % load image
            img = LMimread(database, i, HOMEIMAGES); % Load image
            [nrows ncols c] = size(img);
            if (c==3) % convert to grayscale? (I'm not sure if this is required by OpenCV...)
                img = rgb2gray(img);
            end;

            % crop and add positive examples from image i to positive image
            % directory
            bboxes = []; % store bounding boxes for negative patch generation
            for j = 1:Nobjects
                [X,Y] = getLMpolygon(database(i).annotation.object(j).polygon);

                % compute bounding box of object
                x = min(X)-2;
                y = min(Y)-2;
                w = max(X)+2 - x;
                h = max(Y)+2 - y;
                ctr_x = x+w/2;
                ctr_y = y+h/2;

                % round width and height to be a multiple of the patch size
                % width and height, and re-center bounding box
                sfactor = max(w/pw,h/ph);
                w = floor(sfactor*pw);
                h = floor(sfactor*ph);
                x = floor(ctr_x - w/2);
                y = floor(ctr_y - h/2);

                % make sure that patch fits inside image (otherwise skip
                % this positive)
                if (w>ncols || h>nrows)
                    continue;
                end

                % check boundaries
                if(x<1); x = 1; end;
                if(y<1); y = 1; end;
                if((x+w)>ncols); x = ncols-w; end;
                if((y+h)>nrows); y = nrows-h; end;

                % crop and save image in positive image directory
                img_filename = sprintf('%s/pos%06d.jpg', pbdir, obj_ctr);
                cimg = img(y:(y+h-1),x:(x+w-1));
                imwrite(cimg, img_filename, 'JPG', 'Quality', 100);

                % add entry into positives file:
                % OpenCV format is upper-left corner (x,y), width (w) and
                % height (h)
                fprintf(pos_file, 'positives\\pos%06d.jpg\t1\t0 0 %d %d\n', obj_ctr, w, h);

                % use bboxes below...
                bboxes = [bboxes; x y w h];

                % next object
                obj_ctr = obj_ctr + 1;
            end

            % generate negative patches
            x = [];
            y = [];
            for j = 1:nneg_per_image
                if neg_ctr > num_neg
                    break;
                end

                cx = floor(rand(1)*(ncols-pw-1));
                cy = floor(rand(1)*(nrows-ph-1));

                % check boundaries
                if(cx<1); cx = 1; end;
                if(cy<1); cy = 1; end;
                if((cx+pw)>ncols); cx = ncols-pw; end;
                if((cy+ph)>nrows); cy = nrows-ph; end;

                while (sum(cx==x) || sum(cy==y) || ...
                        isOverlap(cx,cy,pw,ph,bboxes))
                    cx = floor(rand(1)*(ncols-pw-1));
                    cy = floor(rand(1)*(nrows-ph-1));

                    % check boundaries
                    if(cx<1); cx = 1; end;
                    if(cy<1); cy = 1; end;
                    if((cx+pw)>ncols); cx = ncols-pw; end;
                    if((cy+ph)>nrows); cy = nrows-ph; end;
                end

                x = [x cx];
                y = [y cy];

                % crop and save negative patch
                cimg = img(cy:(cy+ph-1),cx:(cx+pw-1));
                img_filename = sprintf('%s/neg%06d.jpg', nbdir, neg_ctr);
                imwrite(cimg, img_filename, 'JPG', 'Quality', 100);

                % update negatives file
                fprintf(neg_file, 'negatives\\neg%06d.jpg\n', neg_ctr);

                neg_ctr = neg_ctr + 1;
            end
        catch
            i
            'dimensions (x, y, w, h)'
            x
            y
            w
            h
        end
    end
end

% done.
fclose(pos_file);
fclose(neg_file);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% isOverlap
%

function ans = isOverlap(x, y, w, h, bboxes)

ans = 0;
n = size(bboxes,1);
for idx = 1:n
    cx = bboxes(idx,1);
    cy = bboxes(idx,2);
    cw = bboxes(idx,3);
    ch = bboxes(idx,4);

    if ~(x>(cx+cw) || y>(cy+ch) || (x+w)<cx || (y+h)<cy)
        ans = 1;
        return;
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
function [tests,pass,perf]=test_wfdbdemo(varargin)

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end
%Test the examples
test_string={'wfdbdemo;close all;'};

[tests,pass,perf]=test_wrapper(test_string,[],verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function value = tetrahedron_unit_monomial ( expon )

%*****************************************************************************80
%
%% TETRAHEDRON_UNIT_MONOMIAL integrates a monomial over the unit tetrahedron.
%
%  Discussion:
%
%    This routine integrates a monomial of the form
%
%      product ( 1 <= dim <= 3 ) x(dim)^expon(dim)
%
%    where the exponents are nonnegative integers.  Note that
%    if the combination 0^0 is encountered, it should be treated
%    as 1.
%
%    Integral ( over unit tetrahedron ) x^l y^m z^n dx dy =
%    l% * m% * n% / ( m + n + 3 )%
%
%    The integration region is defined as:
%
%      0 <= X
%      0 <= Y
%      0 <= Z
%      0 <= X + Y + Z <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer EXPON(3), the exponents.
%
%    Output, real VALUE, the integral of the monomial.
%

%
%  The first computation ends with VALUE = 1.0;
%
  value = 1.0;
%
%  The first loop simply calculates 1, so we short circuit it.
%
% k = 0;
%
% for i = 1 : expon(1)
%   k = k + 1;
%   value = value * i / k;
% end

  k = expon(1);
  for i = 1 : expon(2)
    k = k + 1;
    value = value * i / k;
  end

  for i = 1 : expon(3)
    k = k + 1;
    value = value * i / k;
  end

  k = k + 1;
  value = value / k;

  k = k + 1;
  value = value / k;

  k = k + 1;
  value = value / k;

  return
end


text:
function [ v1, v2, v3, v4 ] = tetrahedron_ref ( )

%*****************************************************************************80
%
%% TETRAHEDRON_REF returns the vertices of the reference tetrahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 July 2014
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Output, real V1(3), V2(3), V3(3), V4(3), the vertices.
%
  v1(1) = - 1.0;
  v1(2) = - 1.0 / sqrt ( 3.0 );
  v1(3) = - 1.0 / sqrt ( 6.0 );

  v2(1) =   0.0;
  v2(2) =   2.0 / sqrt ( 3.0 );
  v2(3) = - 1.0 / sqrt ( 6.0 );

  v3(1) =   1.0;
  v3(2) = - 1.0 / sqrt ( 3.0 );
  v3(3) = - 1.0 / sqrt ( 6.0 );

  v4(1) =   0.0;
  v4(2) =   0.0;
  v4(3) =   3.0 / sqrt ( 6.0 );

  return
end


text:
function r = get_rand(N, M, seed)

% r = get_rand(N, M, seed)
%
% DESC:
% returns a logical vector of length M with N trues uniformly ditributed
%
% AUTHOR
% Marco Zuliani - marco.zuliani@gmail.com
%
% VERSION
% 1.0.1
%
% INPUT:
% N             = number of ones
% M             = number of elements
% seed          = seed of the random number generator
%
% OUTPUT:
% r             = M-dimensional vector with N ones and M-N zeros
%
% HISTORY
% 1.0.0         - ??/??/04 - Initial version
% 1.0.1         - ??/??/06 - Uses logical indexing
% 1.0.2         - 06/25/08 - Fixes the seed of the random number generator

% fix the seed of the random number generator
if (nargin == 3) && ~isempty(seed)
    rand('twister', seed);
end;

if (N > M)
    error('RANSACToolbox:get_rand', 'N should be less or equal than M');
end;

r = [true(1, N) false(1, M-N+1)];

temp = rand(1, M);

[dummy ind] = sort(temp);
r = r(ind);

return

author: RANSAC

repo: RANSAC-Toolbox


text:
%
% In order to obtain the source code for Fingerprint Recognition System Release 5.1
% ( you are using this last release ) please visit my webpage
% http://utenti.lycos.it/matlab/fingerprint4.htm
%
% For any question please email me luigi.rosa@tiscali.it
%
% Release 5.1
% Major features:
% - An improved algorithm used in fingerprint matching, highly recommended for high-performance applications
% - A generalized and optimized version of the algorithm for core point localization
% - Improved database
% - Demo code (protected P-files) available for performance evaluation
%
% Release 5.0
% Major features:
% - The merging technique is greatly improved. This new algorithm was tested on FVC2004 training fingerprint images. Test results are available on request. Please email me in order to obtain them.
% - Parameters for image segmentation are estimated automatically
% - An improved algorithm is used when the rotated FingerCode is added to database. Now this procedure does not introduce any additional noise
% - A faster fingerprint image acquisition when a new fingerprint image is added to database
% - Implementation of 1D and 2D recursive Gabor filtering
% - Optimized pixel-wise orientation field estimation ( 80% faster than Release 4.0 )
% - List of fingerprint databases available on the web
%
% Release 4.0
% Major features:
% - An improved algorithm for core point detection, based on a novel hybrid technique
% - A better fingerprint segmentation which makes use of morphological operations (binary erosion and dilation)
% - Exhaustive documentation of the implemented algorithms
% - Improved GUI
% - Better error management
% - Image enhancement
% - Orientation field estimation
%
% For more informations please visit
% http://utenti.lycos.it/matlab/fingerprint4.htm
%
% Luigi ROSA
% Via Centrale 35
% 67042 Civita Di Bagno
% L'Aquila - ITALY
% mobile +39 3207214179
% email luigi.rosa@tiscali.it
%
%


text:
% Fig. 5.38  Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%script to generate figure 5.38, Piper Dakota autopilot design

clf
numG = 160*conv ([1 2.5],[1 0.7]);
denG = conv([1 5 40],[1  .03  .06]);
sysG = tf(numG,denG);
sysD=tf([1   3],[1   20]);
sysDG=sysD*sysG;
figure(1)
rlocus(sysG)
axis([-30 2 -12 12])
hold on
pause;
rlocus(sysDG)
hold off
pause;
figure(2)
rlocus(sysG)
axis([-1.5 .1    -.6 .6])
hold on
pause;
rlocus(sysDG)
hold off
pause;
Kp = 0.3;
sysH=tf(1,1);
sysT = feedback (Kp*sysG,sysH);
Kc = 1.5;
sysTD=feedback(Kc*sysDG,sysH);
figure(3)
step(sysT)
hold on
step(sysTD)
nicegrid;
hold off


text:
%MTRAJ Multi-axis trajectory between two points
%
% [Q,QD,QDD] = MTRAJ(TFUNC, Q0, QF, M) is a multi-axis trajectory (MxN) varying
% from configuration Q0 (1xN) to QF (1xN) according to the scalar trajectory function
% TFUNC in M steps. Joint velocity and acceleration can be optionally returned as
% QD (MxN) and QDD (MxN) respectively.  The trajectory outputs have one row per
% time step, and one column per axis.
%
% The shape of the trajectory is given by the scalar trajectory function
% TFUNC which is applied to each axis:
%      [S,SD,SDD] = TFUNC(S0, SF, M);
% and possible values of TFUNC include @lspb for a trapezoidal trajectory, or
% @tpoly for a polynomial trajectory.
%
% [Q,QD,QDD] = MTRAJ(TFUNC, Q0, QF, T) as above but T (Mx1) is a time
% vector which dictates the number of points on the trajectory.
%
% Notes::
% - If no output arguments are specified Q, QD, and QDD are plotted.
% - When TFUNC is @tpoly the result is functionally equivalent to JTRAJ except
%   that no initial velocities can be specified. JTRAJ is computationally a little
%   more efficient.
%
% See also JTRAJ, MSTRAJ, LSPB, TPOLY.

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

function [S,Sd,Sdd] = mtraj(tfunc, q0, qf, M)

    if ~isa(tfunc, 'function_handle')
        error('first argument must be a function handle');
    end

    M0 = M;
    if ~isscalar(M)
        M = length(M);
    end
    if numcols(q0) ~= numcols(qf)
        error('must be same number of columns in q0 and qf')
    end

    s = zeros(M, numcols(q0));
    sd = zeros(M, numcols(q0));
    sdd = zeros(M, numcols(q0));

    for i=1:numcols(q0)
        % for each axis
        [s(:,i),sd(:,i),sdd(:,i)] = tfunc(q0(i), qf(i), M);
    end

% - If no output arguments are specified S, SD, and SDD are plotted
%   against time.

    switch nargout
        case 0
            clf

            if isscalar(M0)
                t = [1:M0]';
            else
                t = M0;
            end
            subplot(311)
            plot(t, s); grid; ylabel('s');

            subplot(312)
            plot(t, sd); grid; ylabel('sd');

            subplot(313)
            plot(t, sdd); grid; ylabel('sdd');
            if ~isscalar(M0)
                xlabel('time')
            else
                for c=get(gcf, 'Children');
                    set(c, 'XLim', [1 M0]);
                end
            end
            shg
        case 1
            S = s;
        case 2
            S = s;
            Sd = sd;
        case 3
            S = s;
            Sd = sd;
            Sdd = sdd;
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
function filter = tapas_physio_gui_filter()
% Creates sub-part for preproc filter options in Batch Editor GUI
%
%  tapas_physio_gui_filter
%
%
%   See also

% Author:   Lars Kasper
% Created:  2019-07-05
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under
% the terms of the GNU General Public License (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either
% version 3 or, at your option, any later version). For further details,
% see the file COPYING or <http://www.gnu.org/licenses/>.

%--------------------------------------------------------------------------
% filter_type
%--------------------------------------------------------------------------
filter_type        = cfg_menu;
filter_type.tag    = 'type';
filter_type.name   = 'Filter Type';
filter_type.help   = {'Which infinite impulse response filter shall be used?'
    '   ''Chebychev Type II (cheby2)''    Chebychev Type II filter, use for steep transition from'
    '               start to stop band'
    '   ''Butterworth (butter)''    Butterworth filter, standard filter with maximally flat'
    '               passband (Infinite impulse response), but stronger'
    '               ripples in transition band'
};
filter_type.labels = {'Chebychev Type II', 'Butterworth'};
filter_type.values = {'cheby2', 'butter'};
filter_type.val    = {'cheby2'};

%--------------------------------------------------------------------------
% filter_stopband
%--------------------------------------------------------------------------
filter_stopband         = cfg_entry;
filter_stopband.tag     = 'stopband';
filter_stopband.name    = 'Stopband';
filter_stopband.help    = {
  '[f_min, f_max] frequency interval in Hz of all frequencies, s.th. frequencies'
  '                 outside this band should definitely *NOT* pass the filter'
  '                 Default: [] '
  '                 NOTE: only relevant for ''cheby2'' filter type'
  '                 if empty, and passband is empty, no filtering is performed'
  '                 if empty, but passband exists, stopband interval is'
  '                 10% increased passband interval'
    };
filter_stopband.strtype = 'e';
filter_stopband.num     = [0 Inf];
filter_stopband.val     = {[]};

%--------------------------------------------------------------------------
% filter_passband
%--------------------------------------------------------------------------
filter_passband         = cfg_entry;
filter_passband.tag     = 'passband';
filter_passband.name    = 'Passband';
filter_passband.help    = {
  '[f_min, f_max] frequency interval in Hz of all frequency that should'
  '                 pass the passband filter'
  '                   default: [0.3 9] (to remove slow drifts, breathing'
  '                                    and slice sampling artifacts)'
  '                 if empty, no filtering is performed'
      };
filter_passband.strtype = 'e';
filter_passband.num     = [0 Inf];
filter_passband.val     = {[0.3 9]};

%--------------------------------------------------------------------------
%% filter_no
%--------------------------------------------------------------------------

filter_no         = cfg_branch;
filter_no.tag  = 'no';
filter_no.name = 'No';
filter_no.val  = {};
filter_no.help = {'Cardiac data remains unfiltered before pulse detection.'};

%--------------------------------------------------------------------------
%% filter_yes
%--------------------------------------------------------------------------

filter_yes      = cfg_branch;
filter_yes.tag  = 'yes';
filter_yes.name = 'Yes';
filter_yes.val  = {filter_type, filter_passband, ...
filter_stopband};
filter_yes.help = {''};

%--------------------------------------------------------------------------
%% filter
%--------------------------------------------------------------------------

filter      = cfg_choice;
filter.tag  = 'filter';
filter.name = 'Filter Raw Cardiac Time Series';
filter.val  = {filter_no};
filter.values  = {filter_no, filter_yes};
filter.help = {
    'Filter properties for bandpass-filtering of cardiac signal before peak'
    'detection, phase extraction, and other physiological traces'
 };

author: translationalneuromodeling

repo: tapas


text:
function [U] = hyperNfindr(M, q)
% HYPERNFINDR Performs the N-FINDR (endmember extraction) algorithm
%   Performs the N-FINDR algorithm to find q endmembers. If only M is
%  given as input, this function calls hyperHfcVd to estimate the number
%  of endmembers (q) and then hyperPct to reduce dimensionality to (q-1).
%
% Usage
%   [U] = hyperNfindr(M)
%   [U] = hyperNfindr(M, q)
% Inputs
%   M - 2d matrix of HSI data (p x N)
%   q - Number of endmembers to find
%       -- if not given, q is obtained from hyperHfcVd(M, 10^-3)
% Outputs
%   U - Recovered endmembers (p x q)
%
% References
%   M. Winter, "N-findr: an algorithm for fast autonomous
% spectral endmember determination in hyperspectral data," SPIE’s
% International Symposium on Optical Science, Engineering, and
% Instrumentation, pages 266–275. International Society for Optics
% and Photonics, 1999.

% Error trapping
if ndims(M) ~= 2
    warning('WarnTests:dim', ...
            'Input image must be p x N.\n',...
            'Converting with hyperConvert2d.\n')
    M = hyperConvert2d(M);
end

M_orig = M;
[p, N] = size(M);

if nargin == 1
    fprintf('Implementing hyperHfcVd to determine the number of endmembers.\n')
    q = hyperHfcVd(M_orig, [10^-3]);
    fprintf('Reducing dimensionality to (q-1) using hyperPct.\n')
    M = hyperPct(M, q-1);
elseif q < p+1
    warning('WarnTests:dim', ...
    strcat('N-FINDR requires (q-1) spectral bands.\n',...
           'Performing PCA to reduce dimensionality.\n'))
    M = hyperPct(M, q-1);
elseif q > p+1
    warning('WarnTests:dim', ...
    strcat('N-FINDR requires (q-1) spectral bands.\n',...
           'Performing PCA to reduce dimensionality.\n'))
    error('ErrTests:dim', ...
        strcat('N-FINDR cannot find more than (p+1) endmembers (q),\n', ...
               'where p is the number of available spectral bands.\n'))
end

% Initialize
U_idx = randperm(N,q); % Random endmember selection
E     = M(:,U_idx);    % Endmember matrix
V     = abs(det([ones(1,q); E])) / factorial(q-1); % Simplex volume
vols  = zeros(q,1);

% Search for maximum volume simplex
for j = 1:N;
    % Replace each column of E with sample vector M(:,j)
    % and compute the volume for each
    for k = 1:q;
        E_tmp      = E;
        E_tmp(:,k) = M(:,j);
        vols(k)    = abs(det([ones(1,q); E_tmp])) / factorial(q-1);
    end
    % If max volume is greater than previous V, update E and V
    [V_tmp,k_idx] = max(vols);
    if V_tmp > V
        V            = V_tmp;
        E(:,k_idx)   = M(:,j);
        U_idx(k_idx) = j;
    end
end

% Return endmembers
U = M_orig(:, U_idx);

author: davidkun

repo: HyperSpectralToolbox


text:
function tf = evalbool(s)
%EVALBOOL Evaluate argument to a bool

if ~ischar(s) % S may not and cannot be empty
        tf = s;

else % Evaluation of string S
    if strncmpi(s, 'yes', 3) || strncmpi(s, 'on', 2) ...
        || strncmpi(s, 'true', 4) || strncmp(s, '1 ', 2)
            tf = 1;
    elseif strncmpi(s, 'no', 2) || strncmpi(s, 'off', 3) ...
        || strncmpi(s, 'false', 5) || strncmp(s, '0 ', 2)
            tf = 0;
    else
        try tf = evalin('caller', s); catch
            error(['String value "' s '" cannot be evaluated']);
        end
        try tf ~= 0; catch
            error(['String value "' s '" cannot be evaluated reasonably']);
        end
    end

end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function visMissingRecall_spatial(first_idxs,annolist_gt,xaxis,saveTo)

if (~exist(saveTo,'dir'))
    mkdir(saveTo);
end

for j = 1:length(first_idxs)
    clf;
    img = imread(annolist_gt(j).image.name);
    imagesc(img); axis equal; hold on;
    rects = annolist_gt(j).annorect(first_idxs{j});

    points = [];
    for ridx = 1:length(rects)
        p1 = util_get_annopoint_by_id(rects(ridx).annopoints.point,xaxis(2));
        p2 = util_get_annopoint_by_id(rects(ridx).annopoints.point,xaxis(1));
        points = [points; [p1.x p1.y p2.x p2.y]];
%         if (isfield(rects(ridx), 'annopoints') && ~isempty(rects(ridx).annopoints) && ...
%             isfield(rects(ridx).annopoints, 'point') && ~isempty(rects(ridx).annopoints.point))
%             point = util_get_annopoint_by_id(rects(ridx).annopoints.point,jidx);
%             if (~isempty(point))
%                 points = [points; [point.x point.y]];
%             end
%         end
    end
%     plot(points(:,1),points(:,2),'ro','MarkerFaceColor','y','MarkerEdgeColor','k','MarkerSize',5);
    for i = 1:size(points,1)
        plot([points(i,1),points(i,3)],[points(i,2),points(i,4)],'y-','lineWidth',3);
    end
%     lname = sprintf('detection, score: %1.2f, dist: %1.2f\n',first_scores(j),first_dist(j));
%     legend({lname,'GT'},'fontSize',12);
    axis off;
%     fprintf('imgidx: %d\n',first_idxs(j));
    print(gcf, '-dpng', [saveTo '/imgidx_' padZeros(num2str(j),5) '.png']);
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
classdef AnisotropicContributionTensor < Stiffness3DTensor

    properties (Access = private)
        firstTensor
        secondTensor
        thirdTensor
        dir
    end

    methods (Access = public)

        function obj = AnisotropicContributionTensor(A,direction)
            obj.generateTensors(A,direction);
            obj.addTensors()
        end

        function t = clone(obj)
            tens = obj.firstTensor();
            d = obj.dir;
            t = AnisotropicContributionTensor(tens,d);
        end
    end

    methods (Access = private)

        function generateTensors(obj,A,direction)
            obj.dir          = direction;
            obj.firstTensor  = A;
            obj.secondTensor = SecondAnisotropicTensor(A,obj.dir);
            obj.thirdTensor  = ThirdAnisotropicTensor(A,obj.dir);
        end

        function addTensors(obj)
            t1  = obj.firstTensor.getValue();
            t2  = obj.secondTensor.getValue();
            t3  = obj.thirdTensor.getValue();
            tensor = t1 + t2 + t3;
            obj.setValue(tensor);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% Demo of *format_scat*

%% Usage
% [out,meta] = *format_scat*(S,fmt) (see
% <matlab:doc('format_scat') format_scat>).
%
%% Description
% In this demo, we show how format_scat can be used to gather all
% scattering coefficients into one single array, and visualize this array
% using imagesc.
% NB : so as to improve contrast, the colormap only ranges from 0 to 2,
% setting all higher values in white.

N = 65536;
load handel;
y = y(1:N);
T = 4096;
filt_opt = default_filter_options('audio', T);
scat_opt.M = 2;
[Wop, filters] = wavelet_factory_1d(N, filt_opt, scat_opt);
S = scat(y, Wop);
fmt = 'table'; % default format
X = format_scat(S,fmt);

figure;
coefft = (1:size(X,1));
time = (1:size(X,2)) * T/Fs;
imagesc(time,coefft,X);
colormap bone;
caxis([0 2]); % color thresholding improves contrast
colorbar;
xlabel('Time (seconds)');
ylabel('Coefficients');

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [ state ] =movePTPArcYZ_AC(t,theta,c,vel)
%% This function is used for moving the endeffector on an arc in the YZ plane, for the KUKA iiwa 7 R 800.

%% Syntax:
% [ state ] =movePTPArcYZ_AC(t,theta,c,vel)

%% About:
% This function is used to move the end-effector on an arc in the YZ plane,

%% Arreguments:
% t: is the TCP/IP connection
% theta: is the arc angle, in radians
% c: the YZ coordinates of the center of the circle, it is 1x2 vector.
% vel : is a double, defines the motion velocity mm/sec.

% Copyright, Mohammad SAFEEA, 9th of May 2017

    k=[1;0;0];
    pos=getEEFPos( t );
    c=colVec(c);
    c1=[pos{1};c(1);c(2)];
    state=movePTPArc_AC(t,theta,c1,k,vel);

end

function [ y ] = colVec( v)
% Convert a vector to a column vector:
    if(size(v,2)==1)
        y=v;
    else
        y=v';
    end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
%DRAWMARKERARUCO  Draw a canonical marker image
%
%     img = cv.drawMarkerAruco(dictionary, id, sidePixels)
%     img = cv.drawMarkerAruco(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __dictionary__ dictionary of markers indicating the type of markers.
% * __id__ identifier of the marker that will be returned. It has to be a
%   valid id in the specified dictionary (0-based).
% * __sidePixels__ size of the image in pixels.
%
% ## Output
% * __img__ output image with the marker, of size `sidePixels-by-sidePixels`.
%
% ## Options
% * __BorderBits__ width of the marker border. default 1
%
% This function returns a marker image in its canonical form (i.e. ready to be
% printed).
%
% See also: cv.estimatePoseBoard, cv.drawPlanarBoard, cv.drawCharucoBoard
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%kvdetilt 'Remove Illumination Gradient by Subtracting Best-Fit Plane (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vdetilt.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: m 'Masking Image ', optional: 'find tilt operation masking image'
% OutputFile: o 'Output Image ', required: 'resulting detilted image'
%
% Example: o = kvdetilt({i, m}, {'i','';'m','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vdetilt - Remove Illumination Gradient by Subtracting Best-Fit Plane  (K1)
%
%  DESCRIPTION
% .I vdetilt
% computes the best-fit plane for an image, and then subtracts that plane
% from the image to produce the output.  This is very handy for correcting
% illumination gradients in a poorly digitized image.  The "remove mean" option
% also removes the mean from the image, giving it zero-mean statistics.
% If "remove mean" is set to yes (default), then remove mean from image.
%
% The optional mask image must be of the same data type as the input images,
% and is used to gate the operation.  A non-zero mask pixel enables the output
% pixel to contain the adjusted value. A zero mask pixel just transfers
% the value of the corresponding pixel in the first input image to the output
% pixel.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% vgettilt, lvgettilt, vtilt, lvtilt
%
%  RESTRICTIONS
% .I vdetilt
% will not work on BIT, transform or COMPLEX data storage types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function numOut = analyzeWriteHeader(fname, hdr, endianType)
%
% numOut = analyzeWriteHeader(fname, hdr, [endianType])
%
% writes Analyze-format header
%
% fname   - filename 	     (e.g 'spm' or 'spm.img')
%
% numOut       - number of elements successfully written (should be 348)
%
% see also dbh.h (ANALYZE) spm_hread.m and spm_type.m
%

% Default to little-endian
if ~exist('endianType','var') | isempty(endianType)
    endianType = 'ieee-le';
end

% ensure correct suffix {.hdr} and open header file
%
[p,f,e] = fileparts(fname);
fname = fullfile(p, [f '.hdr']);

fid = fopen(fname,'w',endianType);

if (fid == -1),
	error(['Error opening ' fname '. Check that you have write permission.']);
end;

%---------------------------------------------------------------------------
% if hdr.datatype == 1;   hdr.bitpix = 1;  hdr.glmax = 1;        hdr.glmin = 0;	end
% if hdr.datatype == 2;   hdr.bitpix = 8;  hdr.glmax = 255;      hdr.glmin = 0;	end
% if hdr.datatype == 4;   hdr.bitpix = 16; hdr.glmax = 32767;    hdr.glmin = 0;  	end
% if hdr.datatype == 8;   hdr.bitpix = 32; hdr.glmax = (2^31-1); hdr.glmin = 0;	end
% if hdr.datatype == 16;  hdr.bitpix = 32; hdr.glmax = 1;        hdr.glmin = 0;	end
% if hdr.datatype == 64;  hdr.bitpix = 64; hdr.glmax = 1;        hdr.glmin = 0;	end

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

fseek(fid,0,'bof');

% write (struct) header_key
%
fwrite(fid, 348, 'int32');
fwrite(fid, fixLen(hdr.data_type, 10), 'char' );
fwrite(fid, fixLen(hdr.db_name,18), 'char' );
fwrite(fid, fixLen(hdr.extents,1), 'int32');
fwrite(fid, fixLen(hdr.session_error,1), 'int16');
fwrite(fid, fixLen(hdr.regular,1), 'char' );
fwrite(fid, fixLen(hdr.hkey_un0,1), 'char' );

% write (struct) image_dimension
%---------------------------------------------------------------------------
fseek(fid,40,'bof');

fwrite(fid, fixLen(hdr.dim,8), 'int16');
fwrite(fid, fixLen(hdr.vox_units,4), 'char' );
fwrite(fid, fixLen(hdr.cal_units,8), 'char' );
fwrite(fid, fixLen(hdr.unused1,1), 'int16' );
fwrite(fid, fixLen(hdr.datatype,1), 'int16');
fwrite(fid, fixLen(hdr.bitpix,1), 'int16');
fwrite(fid, fixLen(hdr.dim_un0,1), 'int16');
fwrite(fid, fixLen(hdr.pixdim,8), 'float');
fwrite(fid, fixLen(hdr.vox_offset,1), 'float');
fwrite(fid, fixLen(hdr.funused1,1),	'float');
fwrite(fid, fixLen(hdr.funused2,1),	'float');
fwrite(fid, fixLen(hdr.funused3,1), 'float');
fwrite(fid, fixLen(hdr.cal_max,1), 'float');
fwrite(fid, fixLen(hdr.cal_min,1), 'float');
fwrite(fid, fixLen(hdr.compressed,1), 'int32');
fwrite(fid, fixLen(hdr.verified,1), 'int32');
fwrite(fid, fixLen(hdr.glmax,1), 'int32');
fwrite(fid, fixLen(hdr.glmin,1), 'int32');

% write (struct) data_history
%---------------------------------------------------------------------------
fwrite(fid, fixLen(hdr.descrip,80),	'char');
fwrite(fid, fixLen(hdr.aux_file,24), 'char');
fwrite(fid, fixLen(hdr.orient,1), 'char');
fwrite(fid, fixLen(hdr.origin,5), 'int16');
fwrite(fid, fixLen(hdr.generated,10), 'uchar');
fwrite(fid, fixLen(hdr.scannum,10), 'uchar');
fwrite(fid, fixLen(hdr.patient_id,10), 'uchar');
fwrite(fid, fixLen(hdr.exp_date,10), 'uchar');
fwrite(fid, fixLen(hdr.exp_time,10), 'uchar');
fwrite(fid, fixLen(hdr.hist_un0,3), 'uchar');
fwrite(fid, fixLen(hdr.views,1), 'int32');
fwrite(fid, fixLen(hdr.vols_added,1), 'int32');
fwrite(fid, fixLen(hdr.start_field,1), 'int32');
fwrite(fid, fixLen(hdr.field_skip,1), 'int32');
fwrite(fid, fixLen(hdr.omax,1), 'int32');
fwrite(fid, fixLen(hdr.omin,1), 'int32');
fwrite(fid, fixLen(hdr.smax,1), 'int32');
fwrite(fid, fixLen(hdr.smin,1), 'int32');

numOut = ftell(fid);
fclose(fid);
return;

function m = fixLen(m, len)
m = m(:)';
if(length(m) > len)
    m = m(1:len);
elseif(length(m) < len)
    if(ischar(m))
        m = [m char(zeros(1,len-length(m)))];
    else
        m = [m zeros(1,len-length(m))];
    end
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function sDest = struct_copy_fields(sDest, sSrc, override)
% STRUCT_COPY_FIELDS: Copy the fields from sSrc structure to sDest structure

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

if (nargin < 3) || isempty(override)
    override = 1;
end

% No fields to add
if isempty(sSrc)
    return
% No fields in destination structure
elseif isempty(sDest)
    sDest = sSrc;
% Fields in both structures
else
    namesSrc = fieldnames(sSrc);
    for i = 1:length(namesSrc)
        if override || ~isfield(sDest, namesSrc{i})
            sDest.(namesSrc{i}) = sSrc.(namesSrc{i});
        end
    end
end

author: brainstorm-tools


text:
function value = gamma_log_int ( n )

%*****************************************************************************80
%
%% GAMMA_LOG_INT computes the logarithm of Gamma of an integer N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the logarithm of the Gamma function.
%    0 < N.
%
%    Output, real VALUE, the logarithm of
%    the Gamma function of N.
%
  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GAMMA_LOG_INT - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of N = %d\n', n );
    fprintf ( 1, '  But N must be strictly positive.\n' );
    error ( 'GAMMA_LOG_INT - Fatal error!' );
  end

  value = gammaln ( n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [M, F, b, x, e_conn] =  assemb_co ( param )

%  The FEM equation for the steady state temperature distribution is
%
%    M_2 z^{ss} + F + b = 0

%
%% Initialization
%---------------------------------------------------------------------

  n_nodesx   =  param.nodesx;    % number of nodes in the x direction
  n_nodesy   =  param.nodesy;    % number of nodes in the y direction

  n_gauss    = 7;        % number of points used in Gauss integration
%
%%  Geometry Module
%-----------------------------------------------------------------------
%
%  Generate a mesh of the L \times W rectangle with quadratic elements
%
  n_nodes    = n_nodesx*n_nodesy;
  n_elements = 2*(n_nodesx-1)*(n_nodesy-1)/4;

  x_nodes = linspace(0.0, param.L, n_nodesx);
  y_nodes = linspace(0.0, param.W, n_nodesy);
%
%  Generate node coordinates.
%
  x = zeros(n_nodes, 2);
  for ii=1:n_nodesx
    x_i = x_nodes(ii); k = ii - n_nodesx;
    for jj=1:n_nodesy
      k = k + n_nodesx;
      x(k,:) = [ x_i  y_nodes(jj)];
    end
  end
%
%  Generate element connectivity
%
  e_conn = zeros(n_elements, 6);
  ie     = 0;
%
%  Template for odd numbered triangles
%
  elmt_odd = [ 0 ; 2+2*n_nodesx;   2*n_nodesx;
                    1+n_nodesx ; 1+2*n_nodesx; n_nodesx];
%
%  Template for even numbered triangles
%
  elmt_evn = [ 0 ;      2 ;    2+2*n_nodesx;
                         1;    2+n_nodesx ;    1+n_nodesx];

  for jj=1:2:n_nodesy-1
    for ii=1:2:n_nodesx-1
      k = ii + (jj-1)*n_nodesx;
      ie = ie + 1;
      e_conn(ie,:) = k + elmt_odd;
      ie = ie + 1;
      e_conn(ie,:) = k + elmt_evn;
    end
  end
%
%  Set up arrays for elements with faces on left/right boundary.
%
  Omega_left  =            1:(n_nodesx-1):n_elements;
  Omega_right = (n_nodesx-1):(n_nodesx-1):n_elements;
%
%  With Neumann/Robins bc all nodes are unknown.
%
  n_equations = n_nodes;
%
%  Begin the SPMD block:
%
  spmd
%
%  Set up codistributed structure
%
%  Column pointers and such for codistributed arrays.
%
    Vc  = codistributed.colon(1, n_equations);
    lP  = getLocalPart(Vc);
    lP_1 = lP(1); lP_end = lP(end); % first and last columns of M on this lab

    fid  = fopen(['out' int2str(labindex)], 'w');
    fprintf(fid, 'Lab %d : column range %d  %d \n', labindex, lP_1, lP_end);

    co_dist_Vc = getCodistributor(Vc);
    dPM = co_dist_Vc.Partition;
    col_shft = [0 cumsum(dPM(1:end-1))];
%
%  sparse arrays on each lab
%
    M_lab = sparse(n_equations, dPM(labindex));
    b_lab = sparse(dPM(labindex), 1); F_lab = b_lab;

    [n_elements, nel_dof] = size(e_conn);
    [r, s, w] = twod_gauss(n_gauss);
%
%  Build the finite element matrices - Begin loop over elements
%
    n_el_lab = 0;
    fprintf(fid, 'Lab %d : begin n_el loop \n', labindex);

    for n_el=1:n_elements
%
%  Which nodes are in this element?
%
      nodes_local = e_conn(n_el,:);
%
%  subset of nodes/columns on this lab
%
      lab_nodes_local = extract( nodes_local, lP_1, lP_end);

%     fprintf(fid, '\n \n n_el = %d \n', n_el);
%     fprintf(fid, 'lab_nodes_local has %d rows \n', size(lab_nodes_local, 1));

      if ~isempty(lab_nodes_local) % continue the calculation for this elmnt
        n_el_lab = n_el_lab + 1;
%       fprintf(fid, 'Lab %d : n_el loop location 1 \n', labindex);
%
%  coordinates, shape functions & gradients evaluated at the Gauss pts.
%
        x_local                   = x(nodes_local,:);
        [x_g, w_g, phi, p_x, p_y] = twod_shape(x_local, r, s, w);

%-------------------------------------------------------------------
%% Element contributions to the weak form of the equations
%  First we evaluate the surface integral (2D) terms
%-------------------------------------------------------------------
        M_loc =   twod_bilinear(param.k_x, p_x, p_x, w_g) + ...
                  twod_bilinear(param.k_y, p_y, p_y, w_g);
        src_g  =  source(x_g,  param); % distributed source term
        F_loc  =  twod_f_int( src_g , phi, w_g );
%
%  pre-allocate boundary arrays
%
        b_loc  =   zeros(nel_dof,1);
%       fprintf(fid, 'Lab %d : n_el loop location 2 \n', labindex);

%---------------------------------------------------------------------
%%  Add boundary terms
%---------------------------------------------------------------------
%
%  Add boundary-integral terms for left and right bounary elements
%  Left  boundary n_el = [1:2*(n_nodesx-1):n_elements
%  Right boundary n_el = [2*(n_nodesx-1):2*(n_nodesx-1):n_elements
%
        if any(Omega_left == n_el);
%
%  Here we impose Robin boundary term. The node points on the boundary
%  are nodes_local([3 6 1]) and  are at x_local([3 6 1],:)
%
          [phi_g1, ip, x_g1,  w_g1] = boundary(3, x_local, 5); % n_gauss=5
          alpha_g1 = steps(x_g1(:,2), param.yh, param.ep, param.alpha );
          phi_a    = diag(alpha_g1)*phi_g1;  % mult. each row by alpha(y_g)
          beta_g1  = steps(x_g1(:,2), param.yh, param.ep, param.beta );
%
%  1D integrals are from y = w to y = 0
%
          for ii=1:3
            for jj=1:3
              M_loc(ip(ii),ip(jj)) = M_loc(ip(ii),ip(jj)) + ...
                            oned_f_int( phi_g1(:, jj), phi_a(:,ii), w_g1 );
            end
            b_loc(ip(ii),1) = b_loc(ip(ii),1) + ...
                                  oned_f_int( beta_g1, phi_a(:,ii), w_g1 );
          end

        elseif any(Omega_right == n_el);
%
%  Here we impose boundary term for a specified flux
%  The node points on the boundary are nodes_local([2 3 5]) and
%  are at x_local([2 3 5],:)
%  We need \int f(y) \phi(L, y) dy along the right boundary.
%
          [phi_g1, ip,  ~  ,  w_g1] = boundary(2, x_local, 5); % n_gauss=5;
          for ii=1:3
             b_loc(ip(ii),1) = b_loc(ip(ii),1) + ...
                             oned_f_int( param.flux, phi_g1(:,ii), w_g1 );
          end
        end

%-----------------------------------------------------------------
%% Assemble contributions into the global system matrices
%-----------------------------------------------------------------
%  For M_lab we want all the rows (test functions) but only the
%  columns associated with nodes on this lab
%
        for n_t = 1:nel_dof               % local DOF  - test fcn
          t_glb  = nodes_local(n_t);     % global DOF - test fcn
          for n_u = 1:size(lab_nodes_local, 1)
            n_locj = lab_nodes_local(n_u, 1);  % local DOF in current n_el
            n_glbj = lab_nodes_local(n_u, 2) ...
                         - col_shft(labindex); % global DOF
            M_lab(t_glb,n_glbj) = M_lab(t_glb,n_glbj) ...
                                          +  M_loc(n_t,n_locj);
          end
%
%  Some triangles have nodes on more than one lab;
%  we want only the rows associated with nodes on this lab
%
          if t_glb >= lP_1 && t_glb <= lP_end
            t_loc = t_glb - col_shft(labindex);
            b_lab(t_loc,1)  = b_lab(t_loc,1)  + b_loc(n_t,1);
            F_lab(t_loc,1)  = F_lab(t_loc,1)  + F_loc(n_t,1);
          end

        end % for  n_t

      end % if not empty

    end  % n_el

    fprintf(fid,'Fraction of elements evaluated on this lab  is %8.6f \n',...
                                                      n_el_lab/n_elements);
%
%  Assemble the 'lab' parts in a codistributed format.
%  syntax for version R2009b
%
    codist_matrix = codistributor1d( 2, dPM, [n_equations, n_equations]);
    M = codistributed.build(M_lab, codist_matrix );
    codist_vector = codistributor1d( 1, dPM, [n_equations, 1]);
    b  = codistributed.build(b_lab , codist_vector );
    F  = codistributed.build(F_lab , codist_vector );

    fclose(fid);

  end
%
%  This terminates the SPMD block.
%
  return
end


text:
function [mm] = ft2mm(ft)
% Convert length from feet to millimeters.
% Chad Greene 2012
mm = ft*304.8;

author: Sable

repo: mcbench-benchmarks


text:
function prepare_data_gt(input_files, output_file_name, pca_dim)

%load the hyperspectral data
mat_load=load(input_files.data_file);
fields=fieldnames(mat_load);
fields=fields{1};
im=getfield(mat_load,fields);

%load the labels
gt_load=load(input_files.gt_file);
fields=fieldnames(gt_load);
fields=fields{1};
ip_gt=getfield(gt_load,fields);

%get statistics of the data
[h,w,ch]=size(im)
num_pix=length(find(ip_gt>0));%h*w;

im_cen=im;

%pad the im_cen
conv_size=5;
pad_no=(conv_size-1)/2;

%padded image
im_X=zeros(h+2*pad_no,w+2*pad_no,ch);
for i=1:ch
    im_i=im_cen(:,:,i);
    im_X(:,:,i)=padarray(im_i,[pad_no,pad_no],'symmetric');
end

ip_gt_pad=padarray(ip_gt,[pad_no,pad_no],'symmetric');

cnt=1;
X=zeros(num_pix,ch,conv_size,conv_size);
labels=zeros(num_pix,1);
verbose=false;

for y=1:h
    fprintf('\nRow num=%d of %d',y,h);
    for x=1:w
        if ip_gt(y,x)>0
             X(cnt,:,:,:)=permute(im_X(y:y+conv_size-1,x:x+conv_size-1,:),...
                [3,1,2]);
            labels(cnt)=ip_gt(y,x);

            if verbose
                figure(1);
                imagesc(ip_gt_pad);
                rectangle('Position',[x y conv_size conv_size],...
                    'EdgeColor','r','linewidth',2);
                title(sprintf('Label=%d',labels(cnt)));
            end
            drawnow;
            pause(0.01);
            cnt=cnt+1;
        end
    end
end

comp=pcExtract(permute(X,[2,1,3,4]),pca_dim);
X_r=permute(comp,[2,1,3,4]);

save(output_file_name,'X_r','labels','ip_gt','im','-v7.3');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function v = lexicmp(a,b, tol)

% compare using lexicographical ordering the arrays
%
%   v = lexicmp(a,b, tol);
%
%   tol (default 1e-9) is a tolerance of equality.
%
%   return -1 if a<b, +1 if a>b, 0 if a=b.
%
%   Copyright (c) 2008 Gabriel Peyre

tol = 1e-9;
a = a(:);
b = b(:);

n = max(length(a),length(b));
a(n+1:end) = -Inf;
b(n+1:end) = -Inf;

for i=1:n
    if a(i)<b(i)-tol
        v = -1; return;
    end
    if a(i)>b(i)+tol
        v = +1; return;
    end
end
v = 0;

author: gpeyre

repo: matlab-toolboxes


text:
function [mask, fvc]=VObjEllipsoid(p,x,y,z,flag)
%create 3D ellipsoid virtual object

% Initialize parameters
RadiusX=p.RadiusX;
RadiusY=p.RadiusY;
RadiusZ=p.RadiusZ;
CenterX=p.CenterX;
CenterY=p.CenterY;
CenterZ=p.CenterZ;
FaceNum=p.FaceNum;

% Generate ellipsoid coordinate
[X,Y,Z] = ellipsoid(CenterX,CenterY,CenterZ,RadiusX,RadiusY,RadiusZ,FaceNum);

fvc = surf2patch(X,Y,Z);

if flag~=1 % Render object only
    mask=0;
    return;
end

mask = vert2mask(fvc.vertices,x,y,z);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function tests = test_spm_dcm_post_hoc
% Unit Tests for spm_dcm_post_hoc
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_dcm_post_hoc.m 6492 2015-06-26 14:27:40Z guillaume $

tests = functiontests(localfunctions);

%--------------------------------------------------------------------------
function data_path = get_data_path()
data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'test_spm_dcm_post_hoc');

%--------------------------------------------------------------------------
function setup(testCase)
data_path = get_data_path();

% Expected outputs
artefacts = {fullfile(data_path,'DCM_BPA.mat');
             fullfile(data_path,'DCM_opt_fwd_bwd_simulated.mat')};

% Delete if exist
for i = 1:length(artefacts)
    spm_unlink(artefacts{i});
end

% Initialize SPM
spm('defaults','fmri');
spm_get_defaults('cmdline',true);

%--------------------------------------------------------------------------
function test_on_simulated_attention_data(testCase)
import matlab.unittest.constraints.*

data_path = get_data_path();

% Load model matching the generative model
DCM_fwd = load(fullfile(data_path, 'DCM_fwd_simulated.mat'));
DCM_fwd = DCM_fwd.DCM;

% Load model with extra connection
DCM_fwd_bwd = load(fullfile(data_path, 'DCM_fwd_bwd_simulated.mat'));
DCM_fwd_bwd = DCM_fwd_bwd.DCM;

% Run post-hoc
DCM_BPA = spm_dcm_post_hoc(fullfile(data_path, 'DCM_fwd_bwd_simulated.mat'));
DCM_opt = load(fullfile(data_path,'DCM_opt_fwd_bwd_simulated.mat'));
DCM_opt = DCM_opt.DCM;

% Check artefacts were created
testCase.assertEqual( exist(fullfile(data_path,'DCM_BPA.mat'),'file'), 2 );
testCase.assertEqual( exist(fullfile(data_path,'DCM_opt_fwd_bwd_simulated.mat'),'file'), 2 );

% Check reduced priors match the generative model
pC_fwd = full(diag(DCM_fwd.M.pC));
pC_opt = full(diag(DCM_opt.M.pC));
testCase.assertTrue(all(pC_fwd == pC_opt));

% % Plot for visual inspection
% spm_figure('Getwin','BMC'); clf
%
% subplot(2,2,1); spm_plot_ci(DCM_fwd.Ep,DCM_fwd.Cp);
% title('Fwd (generative) model'); xlabel('Parameter')
%
% subplot(2,2,2); spm_plot_ci(DCM_fwd_bwd.Ep,DCM_fwd_bwd.Cp);
% title('Bwd (over-complex) model'); xlabel('Parameter')
%
% subplot(2,2,3); spm_plot_ci(DCM_BPA.Ep,DCM_BPA.Cp);
% title('BPA'); xlabel('Parameter')
%
% subplot(2,2,4); spm_plot_ci(DCM_opt.Ep,DCM_opt.Cp);
% title('Optimal model for subject'); xlabel('Parameter');

author: spm


text:
function nleqs_master_ex2(alg)
if nargin < 1
    alg = 'DEFAULT';
end
x0 = [1; 2];
opt = struct( ...
    'verbose', 2, ...
    'alg', alg, ...
    'fd_opt', struct( ...
        'jac_approx_fcn', @jac_approx_fcn2, ...
        'labels', {{'P','Q'}}), ...
    'gs_opt', struct('x_update_fcn', @x_update_fcn2) );
[x, f, exitflag, output] = nleqs_master(@f2, x0, opt);
fprintf('\nexitflag = %d\n', exitflag);
fprintf('\nx = \n');
fprintf('   %2g\n', x);
fprintf('\nf = \n');
fprintf('   %12g\n', f);

function [f, J] = f2(x)
%% from Christi Patton Luks, https://www.youtube.com/watch?v=pJG4yhtgerg
f = [  x(1)^2 +   x(1)*x(2)   - 10;
       x(2)   + 3*x(1)*x(2)^2 - 57  ];
if nargout > 1
    J = [   2*x(1)+x(2)    x(1);
            3*x(2)^2       6*x(1)*x(2)+1    ];
end

function JJ = jac_approx_fcn2()
%% for use with fast-decoupled Newton's method
J = [7 2; 27 37];
JJ = {J(1,1), J(2,2)};

function x = x_update_fcn2(x, f)
%% for use with Gauss-Seidel method
x(1) = sqrt(10 - x(1)*x(2));
x(2) = sqrt((57-x(2))/3/x(1));

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [pdc, pdcvar, n] = ft_connectivity_pdc(inputdata, varargin)

% FT_CONNECTIVITY_PDC computes partial directed coherence. This function implements
% the metrices described in Baccala et al., Biological Cybernetics 2001, 84(6),
% 463-74. and in Baccala et al., 15th Int.Conf.on DSP 2007, 163-66.
%
% The implemented algorithm has been tested against the implementation in the
% SIFT-toolbox. It yields numerically identical results to what is known there as
% 'nPDC' (for PDC) and 'GPDC' for generalized pdc.
%
% Use as
%   [p, v, n] = ft_connectivity_pdc(inputdata, ...)
%
% The input data should be a spectral transfer matrix organized as
%   Nrpt x Nchan x Nchan x Nfreq (x Ntime),
% where Nrpt can be 1.
%
% Additional optional input arguments come as key-value pairs:
%   'hasjack'   = 0 (default) is a boolean specifying whether the input
%                 contains leave-one-outs, required for correct variance
%                 estimate
%   'invfun'    = 'inv' (default) or 'pinv', the function used to invert the
%                 transfer matrix to obtain the fourier transform of the
%                 MVAR coefficients. Use 'pinv' if the data are
%                 poorly-conditioned.
%   'noisecov'  = matrix containing the covariance of the residuals of the
%                 MVAR model. If this matrix is defined, the function
%                 returns the generalized partial directed coherence.
%   'feedback'  = 'none', 'text', 'textbar', 'dial', 'etf', 'gui' type of feedback showing progress of computation, see FT_PROGRESS
%
% Output arguments:
%   p = partial directed coherence matrix Nchan x Nchan x Nfreq (x Ntime).
%       If multiple observations in the input, the average is returned.
%   v = variance of p across observations.
%   n = number of observations.
%
% Typically, nrpt should be 1 (where the spectral transfer matrix is
% computed across observations. When nrpt>1 and hasjack is true the input
% is assumed to contain the leave-one-out estimates of H, thus a more
% reliable estimate of the relevant quantities.
%
% See also CONNECTIVITY, FT_CONNECTIVITYANALYSIS

% Copyright (C) 2009-2017, Jan-Mathijs Schoffelen
%
% 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$

hasjack  = ft_getopt(varargin, 'hasjack', 0);
feedback = ft_getopt(varargin, 'feedback', 'none');
invfun   = ft_getopt(varargin, 'invfun', 'inv');
noisecov = ft_getopt(varargin, 'noisecov');

switch invfun
  case {'inv' 'pinv'}
    invfun = str2func(invfun);
  otherwise
    ft_error('unknown specification of inversion-function for the transfer matrix');
end

% crossterms are described by chan_chan_therest
siz = [size(inputdata) 1];
n   = siz(1);

if ~isempty(noisecov)
  ft_progress('init', feedback, 'computing generalized partial directed coherence...');
  scale = diag(sqrt(1./diag(noisecov))); % get the 1./sqrt(var) of the signals
else
  ft_progress('init', feedback, 'computing partial directed coherence...');
  scale = eye(siz(2));
end

% pre-allocate some variables
outsum = zeros(siz(2:end));
outssq = zeros(siz(2:end));

% the mathematics for pdc is most straightforward using the inverse of the
% transfer function
pdim     = prod(siz(4:end));
tmpinput = reshape(inputdata, [siz(1:3) pdim]);
ft_progress('init', feedback, 'inverting the transfer function...');
for k = 1:n
  ft_progress(k/n, 'inverting the transfer function for replicate %d from %d\n', k, n);
  tmp = reshape(tmpinput(k,:,:,:), [siz(2:3) pdim]);
  for m = 1:pdim
    tmp(:,:,m) = scale*invfun(tmp(:,:,m));
  end
  tmpinput(k,:,:,:) = tmp;
end
ft_progress('close');
inputdata = reshape(tmpinput, siz);

for j = 1:n
  ft_progress(j/n, 'computing metric for replicate %d from %d\n', j, n);
  invh   = reshape(inputdata(j,:,:,:,:), siz(2:end));

  den    = sum(abs(invh).^2,1);
  tmppdc = abs(invh)./sqrt(repmat(den, [siz(2) 1 1 1 1]));

  outsum = outsum + tmppdc;
  outssq = outssq + tmppdc.^2;
end
ft_progress('close');

pdc = outsum./n;

if n>1
  if hasjack
    bias = (n-1).^2;
  else
    bias = 1;
  end
  pdcvar = bias*(outssq - (outsum.^2)./n)./(n - 1);
else
  pdcvar = [];
end

% this is to achieve the same convention for all connectivity metrics:
% row -> column
for k = 1:prod(siz(4:end))
  pdc(:,:,k)    = transpose(pdc(:,:,k));
  if ~isempty(pdcvar)
    pdcvar(:,:,k) = transpose(pdcvar(:,:,k));
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function datapt = ma_GAAeroTasks(stateLogEntry, subTask, celBodyData)
%ma_GAAeroTasks Summary of this function goes here
%   Detailed explanation goes here

    bodyID = stateLogEntry(8);

    bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
    ut = stateLogEntry(1);
    rVectECI = stateLogEntry(2:4)';
    vVectECI = stateLogEntry(5:7)';

    switch subTask
        case 'dynPress'
            altitude = norm(rVectECI) - bodyInfo.radius;

            if(altitude <= bodyInfo.atmohgt && altitude >= 0)
                [lat, long, ~, ~, ~, ~, ~, vVectECEF] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo, vVectECI);
                density = getAtmoDensityAtAltitude(bodyInfo, altitude, lat, ut, long);
            elseif(altitude <= 0)
                density = 0;
                vVectECEF = [0;0;0];
            else
                density = 0;
                vVectECEF = [0;0;0];
            end

            vVectEcefMag = norm(vVectECEF);
            vVectEcefMagMS = vVectEcefMag*1000;

            dynP = density * (vVectEcefMagMS^2) / 2; %kg/m^3 * m^2 / s^2  = kg/(m*s^2)
            dynP_kPa = dynP/1000;

            datapt = dynP_kPa;
        case 'atmoPress'
            altitude = norm(rVectECI) - bodyInfo.radius;
            pressure = getPressureAtAltitude(bodyInfo, altitude);

            datapt = pressure;

        case 'atmoTemp'
            altitude = norm(rVectECI) - bodyInfo.radius;

            if(altitude <= bodyInfo.atmohgt && altitude >= 0)
                [lat, long, ~, ~, ~, ~, ~, ~] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo, vVectECI);
                temperature = getTemperatureAtAltitude(bodyInfo, altitude, lat, ut, long) ;
            elseif(altitude <= 0)
                temperature = 0;
            else
                temperature = 0;
            end

            datapt = temperature;

        case 'atmoDensity'
            altitude = norm(rVectECI) - bodyInfo.radius;

            if(altitude <= bodyInfo.atmohgt && altitude >= 0)
                [lat, long, ~, ~, ~, ~, ~, ~] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo, vVectECI);
                density = getAtmoDensityAtAltitude(bodyInfo, altitude, lat, ut, long);
            elseif(altitude <= 0)
                density = 0;
            else
                density = 0;
            end

            datapt = density;
        case 'machNumber'
            altitude = norm(rVectECI) - bodyInfo.radius;
            pressure = getPressureAtAltitude(bodyInfo, altitude);

            if(altitude <= bodyInfo.atmohgt && altitude >= 0)
                [lat, long, ~, ~, ~, ~, ~, vVectECEF, ~] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo, vVectECI);
                density = getAtmoDensityAtAltitude(bodyInfo, altitude, lat, ut, long);
            elseif(altitude <= 0)
                density = 0;
            else
                density = 0;
            end

            pressurePa = pressure*1000; %kPa -> Pa

            if(density > 0)
                speedSound = sqrt(1.4 * pressurePa / density);
                vECEFMag = norm(vVectECEF) * 1000; %km/s -> m/s

                datapt = vECEFMag/speedSound;
            else
                datapt = 0;
            end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

%GEOMETRICAL PARAMETERS

l0 = 30;
h0 = 6.5;

h1 = 10;
h2 = 19;

lm = 30;
hm = 23;

lb = l0 - 8.5;

%Minimum radius of the camshaft
v = hm - h2;

%Y coordinate variation between min and max positions of the camshaft
dv = h2 - h1;

%BREATHING CYCLE PARAMETERS

%Duration of the inhale cycle / duration of the whole cycle
lambda1 = 0.50;
lambda2 = 0.08;

%Soft transition between inhale and exhale cycles
dpsi21 = 0.12;
dpsi12 = 0.15;

%Adjust parameters of the inhale curve
ga1 = 4;
gb1 = .85;
ff1 = 50;

%Adjust parametes of the exhale curve
ga2 = 3.5;
gb2 = .70;
ff2 = 55;

%GENERATION OF THE CURVES AND THE CAMSHAFT

%Time increment
dt = 0.01;

%time coordinate during the whole cycle
theta = 0:dt:2*pi;

%Generation of the soft transition between inhale and exhale curves
psi1 = [];
psi2 = [];

for i = 1:numel(theta);

    if theta(i) < (lambda2-dpsi21/2)*2*pi;
        psi1(i) = 0;

    elseif theta(i) < (lambda2+dpsi21/2)*2*pi;
        psi1(i) = 0.5 + 0.5*cos((theta(i)-(lambda2+dpsi21/2)*2*pi)/(2*dpsi21));

    elseif theta(i) < (lambda1-dpsi12/2)*2*pi;
        psi1(i) = 1;

    elseif theta(i) < (lambda1+dpsi12/2)*2*pi;
        psi1(i) = 0.5 - 0.5*cos((theta(i)-(lambda1+dpsi12/2)*2*pi)/(2*dpsi12));

    else
        psi1(i) = 0;

    end;

    psi2(i) = 1 - psi1(i);

end;

%Generation of the complete breathing cycle rho(theta)

rho = [];
rho1 = [];
rho2 = [];
rho2next = [];

rhomin = 1000;
rhomax = 0;

for i = 1:numel(theta);

    %Inhale curve
    rho1(i) = ff1*gampdf(theta(i), ga1, gb1);

    %Exhale curve
    rho2(i) = ff2*gampdf(theta(i), ga2, gb2);
    rho2next(i) = ff2*gampdf(theta(i)+2*pi, ga2, gb2);

    if theta(i) < lambda1*pi
        rho(i) = psi1(i)*rho1(i) + psi2(i)*rho2next(i);

    else
        rho(i) = psi1(i)*rho1(i) + psi2(i)*rho2(i);

    end;

    %Capturing min and max in order to generate the normalized curve
    if rho(i) > rhomax
        rhomax = rho(i);
    end;

    if rho(i) < rhomin
        rhomin = rho(i);
    end;

end;

%Generation of a normalised curve and camshaft

rhonorm = [];
rhocam = [];

for i = 1:numel(theta);

    rhonorm(i) = (rho(i)-rhomin)/(rhomax-rhomin);

    rhocam(i) = v + rhonorm(i)*dv;

end;

%Generation of the first derivate of the camshaft geometry to analize and
%validate the design

drho = [];
drhonorm = [];
drhocam = [];

a = v;
b = 1/dt;

for i = 1 : numel(rho)-1;

    drho(i) = b*(rho(i+1)-rho(i));
    drhonorm(i) = b*(rhonorm(i+1)-rhonorm(i));
    drhocam(i) = b*(rhocam(i+1)-rhocam(i)) + a;

end;

drho(numel(rho)) = b*(rho(1)-rho(numel(rho)));
drhonorm(numel(rhonorm)) = b*(rhonorm(1)-rhonorm(numel(rhonorm)));
drhocam(numel(rhocam)) = b*(rhocam(1)-rhocam(numel(rhocam))) + a;

%X and Y coordinates during two cycles (for 2-cycle camshaft plot)

theta2 = [theta/2, (2*pi+theta)/2];

rhocam2 = [rhocam, rhocam];
drhocam2 = [drho, drho];

%X and Y coordinates during three cycles (for 3-cycle camshaft plot)

theta3 = [theta/3, (2*pi+theta)/3, (4*pi+theta)/3];

rhocam3 = [rhocam, rhocam, rhocam];
drhocam3 = [drho, drho, drho];

%GENERATION OF THE PLOTS

%Trying to print it out in real size (FAIL)
%set(gcf,'PaperUnits','centimeters');
%set(gcf,'PaperSize',[42 29.7]);

%Plot of the breathing cycle interpolation by curves

fcycle = figure('Name', 'Breath Cycle curves', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
plot(theta, psi1, '-k')
plot(theta, psi2, '-k')
plot(theta, rho1, '-g')
plot(theta, rho2, '-g')
plot(theta, rho2next, '-g')
%plot(theta, drho, '-r')
plot(theta, rho, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of the normalized breathing cycle and first derivate

fnorm = figure('Name', 'Normalized Breath Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
plot(theta, drhonorm, '-r')
plot(theta, rhonorm, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of a 1-cycle camshaft

fcam1 = figure('Name', '1-Cycle Camshaft', 'Units', 'centimeters', 'NumberTitle', 'off');

%hold on
%polar(theta, drhocam, '-r')
polar(theta, rhocam, '-b')
%hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of a 2-cycle camshaft

fcam2 = figure('Name', '2-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');

%hold on;
%polar(theta, drhocam2, '-r')
polar(theta2, rhocam2, '-b')
%hold off;

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of a 3-cycle camshaft

fcam3 = figure('Name', '3-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');

%hold on;
%polar(theta,drhocam3)
polar(theta3, rhocam3)
%hold off;

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

author: ProtofyTeam

repo: OxyGEN

save_path: github-repos/MATLAB/ProtofyTeam-OxyGEN


text:
function snf=snirfcreate(varargin)
%
%    snf=snirfcreate
%       or
%    snf=snirfcreate(option)
%    snf=snirfcreate('Format',format,'Param1',value1, 'Param2',value2,...)
%
%    Create a empty SNIRF data structure defined in the SNIRF
%    specification: https://github.com/fNIRS/snirf
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        option (optional): option can be ignored. If it is a string with a
%             value 'snirf', this creates a default SNIRF data structure;
%             otherwise, a JSNIRF data structure is created.
%        format: save as option.
%        param/value:   a list of name/value pairs specify
%             additional subfields to be stored under the /nirs object.
%
%    output:
%        snf: a default SNIRF or JSNIRF data structure.
%
%    example:
%        snf=snirfcreate('data',mydata,'aux',myauxdata,'comment','test');
%
%    this file is part of JSNIRF specification: https://github.com/fangq/snirf
%
%    License: GPLv3 or Apache 2.0, see https://github.com/NeuroJSON/jsnirf for details
%

if(nargin==1)
    snf=jsnirfcreate(varargin{:});
else
    snf=jsnirfcreate('Format','snirf',varargin{:});
end