text:
function reg = modifiedconv(tr,condf,varargin)
% :Usage:
% ::
%
%     model = modifiedconv(tr,condf,heighteq [all opt],delayeq,ttopeakeq,uonseteq)
%
% :Inputs:
%
%   **tr:**
%        repetition time (sampling rate) of scanning, in seconds
%
%   **condf:**
%        condition function
%         an indicator vector of zeros and ones, where ones indicate event
%         onsets
%
% This is a simplified modified convolution that uses nonlinear saturation in height only
% with a guess as to what the decrease in saturation is as a function of
% the time since previous stimulation (exponential model, alpha version)
%
% For the full set of equations modifying height, time-to-peak, and
% dispersion as a function of stumulus history, see
% modifiedconv_wager2005.m
%
%
% :Examples:
% ::
%
%    condf = [1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]';
%    X = modifiedconv(2,condf);
%    % X is convolved predictor
%    plot(X)
%    X2 = conv(condf,spm_hrf(2)./max(spm_hrf(2)));
%    hold on;
%    plot(X2,'r');
%    legend({'Modified' 'Linear'})
%
% :Please see:
% Wager, T. D., Hernandez, L., Vasquez, A., Nichols, T., and Noll, D.
% C. (2005). Accounting for nonlinear BOLD effects in fMRI: Parameter
% estimates and model for accurate prediction in variable-duration blocked
% and rapid event-related studies.  Neuroimage.
%
% ..
%    06/20/01 Tor Wager
% ..

%heighteq = []; delayeq = []; peakeq = []; uonseteq = [];

% ---------------------------------------------------------------------
% * defaults
% ---------------------------------------------------------------------
height = 1; delay = 0; peak = 6; uonset = 16;
dispers = 1; udisp = 1; rtou = 6; klength = 32;
rtou = Inf; % get rid of undershoot

p = [peak uonset dispers udisp rtou delay klength];

hrf = spm_hrf(tr,p)./ max(spm_hrf(tr,p));

%heighteq = inline('1.7141.*(exp(-2.1038.*x)) + 0.4932.*(exp(-0.0770.*x))');
%delayeq = inline('-13.4097.*(exp(-1.0746.*x)) + 4.8733.*(exp(-0.1979.*x))');
%peakeq = inline('37.5445.*(exp(-2.6760.*x)) + -3.2046.*(exp(-0.2120.*x)) + 5.6344');%

% new idea: model nonlinear saturation response as a function of
% a fixed saturation response (that we know is e.g., 1-(.6658 x original) with a
% stimulus 1 s before) x an exponential discount factor for how LONG AGO the
% previous stimulus occurred
% SO THAT: Rnow = Ro - sum(s(t-to)gamm(t-to)
% gam(1) should be 1-.6658, and area under gam should be no more than 1,
% because you can never lose more than 100% of the signal
% BUT we may end up having to take interactions betwn stim into account,
% which is harder: e.g., a stim occurs .5 s before and another occurs .3
% sec before.  The .3 one is going to produce a smaller response and thus
% less saturation, and you can never saturate more than 1 (100%) even w/
% many stimuli occurring in a short time frame (w/i 1 s)
% SO this model, now, can predict negative response values with more than
% 100% sat.

gam = inline('(1-.6658) * (1./exp(-a))*exp(-a*t)','t','a');
a = .7;  % let's assume this is the exp for now. gives 70% sat w/hist 1:5, 74% w/hist 1:10

myzeros = zeros(length(condf),1);
mylen = length(myzeros);
numels = 12 ./ tr;  % number of elements to count

reg = myzeros;

whstim = find(condf);   % which elements contain indicators of stimulation

if any(condf(whstim)>1)
    disp('warning: modifiedconv not valid for onset mag > 1, which you appear to have entered.'),
end

for i = 1:length(whstim)

    j = whstim(i);  % index of which element

    % for each element, get predicted height

			trialdelta = myzeros;
			trialdelta(j) = 1;
			trialp = p;

			% figure out how many of same type came before
            % 1 is "first stim in sequence"

            %myc = condf(j-min(29,j)+1:j)
            %myw = timeweights(end-length(myc)+1:end)

            myc = condf(j-min(numels,j)+1:j);   % recent events, including current stim

            wh = find(myc(1:end-1));            % which elements

            times = length(myc) - wh;  % times at which these occurred in elements
            times = times * tr;        % convert to seconds

            s = myc(wh);               % stimulus intensity at each time, probably 1 for typical indicator vector

            sat = sum(s .* gam(times,a));  % total saturation

            height = max(0,1 - sat);           % magnitude of this event

            trialdelta = myzeros;
            trialdelta(j) = condf(j);

			mytrialpred = conv(height.*hrf,trialdelta);
			mytrialpred = mytrialpred(1:mylen);

			reg = reg + mytrialpred;

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [FV] = mesh_shrink(FV,origin,dist),

% mesh_shrink - implode vertices of mesh by specific distance
%
% FV = mesh_shrink(FV,origin,dist)
%
% FV is a struct with fields:
%
% FV.vertices   - Nx3 matrix of Cartesian vertex coordindates (X,Y,Z)
% FV.faces      - Mx3 matrix of triangulation of FV.vertices
%
% origin        - 1x3 row vector, usually (0,0,0)
%
% dist          - how far to implode the mesh toward the origin;
%                 this distance is relative to current distance from
%                 the origin, not the total distance from the origin.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no implied or express warranties
% History:  10/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    xo = origin(1); yo = origin(2); zo = origin(3);

    Nvert = size(FV.vertices,1);

    fprintf('...mesh implosion...'); tic;

    for v = 1:Nvert,

        x = FV.vertices(v,1);
        y = FV.vertices(v,2);
        z = FV.vertices(v,3);

        % Find direction cosines for line from centre to vertex
        d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );

        l = (x-xo)/d; % cos alpha
        m = (y-yo)/d; % cos beta
        n = (z-zo)/d; % cos gamma

        % now decrease d by dist
        d = d - dist;

        % locate vertex at this new distance
        x = (l * d) + xo;
        y = (m * d) + yo;
        z = (n * d) + zo;

        FV.vertices(v,:) = [ x y z ];
    end

    t = toc; fprintf('...done (%5.2f sec)\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
%SerialLink.cinertia Cartesian inertia matrix
%
% M = R.cinertia(Q) is the NxN Cartesian (operational space) inertia matrix which relates
% Cartesian force/torque to Cartesian acceleration at the joint configuration Q, and N
% is the number of robot joints.
%
% See also SerialLink.inertia, SerialLink.rne.

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, 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 Mx = cinertia(robot, q)
	J = jacob0(robot, q);
	Ji = inv(J);
	M = inertia(robot, q);
	Mx = Ji' * M * Ji;

author: petercorke

repo: robotics-toolbox-matlab


text:
function [BigPopulation,tempPara] = GenerateBigPopulation(PV,groups,Archive)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Huangke Chen

   cDecs  = Archive.decs;
   PVDecs = cDecs(:,PV);

   BigPopulation = cell(1,size(PVDecs,1));	% initialize the big population to contain all the population
   cPopSize      = 16;                      % the population size of each population

   para.lambda = cPopSize;                                      % population size, offspring number
   tempPara    = repmat(para,size(PVDecs,1),length(groups));	% Initialize the parameters for CMA-ES

   for ci = 1: size(PVDecs,1)
       % Construct the population
       BigPopulation{ci} = repmat(cDecs(ci,:),cPopSize,1);

       % Set the CMA-ES parameters
       for g = 1 : length(groups)	% Set the parameters for each variable group
           DVindex = groups{g};     % the variable index of this group
           N = length(DVindex);     % variable number

           % set the CMA-ES parameters for each group
           tempPara(ci,g).N = N;
           mu      = cPopSize/2;                % number of parents/points for recombination
           weights = log(mu+1/2) - log(1:mu)';	% muXone array for weighted recombination
           tempPara(ci,g).mu = floor(mu);

           tempPara(ci,g).weights = weights/sum(weights);	% normalize recombination weights array
           mueff = sum(weights)^2/sum(weights.^2);          % variance-effectiveness of sum w_i x_i
           tempPara(ci,g).mueff = mueff;

           tempPara(ci,g).cc    = (4+mueff/N)/(N+4+2*mueff/N);                                  % time constant for cumulation for C
           tempPara(ci,g).cs    = (mueff+2)/(N+mueff+5);                                        % t-const for cumulation for sigma control
           tempPara(ci,g).c1    = 2/((N+1.3)^2+mueff);                                          % learning rate for rank-one update of C
           tempPara(ci,g).cmu   = min(1-tempPara(ci,g).c1,2*(mueff-2+1/mueff)/((N+2)^2+mueff));	% and for rank-mu update
           tempPara(ci,g).damps = 1 + 2*max(0,sqrt((mueff-1)/(N+1))-1) + tempPara(ci, g).cs;    % damping for sigma
           tempPara(ci,g).chiN  = N^0.5*(1-1/(4*N)+1/(21*N^2));                                 % expectation of ||N(0,I)|| == norm(randn(N,1))

           % variable parameters
           paraDecs = BigPopulation{ci};
           tempPara(ci,g).xmean = mean(paraDecs(:,DVindex))';
           tempPara(ci,g).sigma = 0.1;          % coordinate wise standard deviation (step size)
           tempPara(ci,g).pc    = zeros(N,1);	% evolution paths for C and sigma
           tempPara(ci,g).ps    = zeros(N,1);	% evolution paths for C and sigma
           B = eye(N,N);	% B defines the coordinate system
           tempPara(ci,g).B = B;
           D = ones(N,1);	% diagonal D defines the scaling
           tempPara(ci,g).D = D;
           tempPara(ci,g).C = B*diag(D.^2)*B';	% covariance matrix C
           tempPara(ci,g).eigeneval = 0;
           tempPara(ci,g).counteval = 0;
       end
   end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout = edgesEvalDir( varargin )
% Calculate edge precision/recall results for directory of edge images.
%
% Enhanced replacement for boundaryBench() from BSDS500 code:
%  http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/
% Uses same format for results and is fully compatible with boundaryBench.
% Given default prms results are *identical* to boundaryBench with the
% additional 9th output of R50 (recall at 50% precision).
%
% The odsF/P/R/T are results at the ODS (optimal dataset scale).
% The oisF/P/R are results at the OIS (optimal image scale).
% Naming convention: P=precision, R=recall, F=2/(1/P+1/R), T=threshold.
%
% In addition to the outputs, edgesEvalDir() creates three files:
%  eval_bdry_img.txt - per image OIS results [imgId T R P F]
%  eval_bdry_thr.txt - per threshold ODS results [T R P F]
%  eval_bdry.txt     - complete results (*re-ordered* copy of output)
% These files are identical to the ones created by boundaryBench.
%
% USAGE
%  [odsF,odsP,odsR,odsT,oisF,oisP,oisR,AP,R50] = edgesEvalDir( prms )
%  [ODS,~,~,~,OIS,~,~,AP,R50] = edgesEvalDir( prms )
%
% INPUTS
%  prms       - parameters (struct or name/value pairs)
%   .resDir     - ['REQ'] dir containing edge detection results (.png)
%   .gtDir      - ['REQ'] dir containing ground truth (.mat)
%   .pDistr     - [{'type','parfor'}] parameters for fevalDistr
%   .cleanup    - [0] if true delete temporary files
%   .thrs       - [99] number or vector of thresholds for evaluation
%   .maxDist    - [.0075] maximum tolerance for edge match
%   .thin       - [1] if true thin boundary maps
%   .ids        - list of images to evaluate on
%
% OUTPUTS
%   odsF/P/R/T  - F-measure, precision, recall and threshold at ODS
%   oisF/P/R    - F-measure, precision, and recall at OIS
%   AP          - average precision
%   R50         - recall at 50% precision
%
% EXAMPLE
%
% See also edgesEvalImg, edgesEvalPlot
%
% Structured Edge Detection Toolbox      Version 3.01
% Code written by Piotr Dollar, 2014.
% Licensed under the MSR-LA Full Rights License [see license.txt]

% get additional parameters
dfs={ 'resDir','REQ', 'gtDir','REQ', 'pDistr',{{'type','parfor'}}, ...
  'cleanup',0, 'thrs',99, 'maxDist',.0075, 'thin',1,'ids',[]};
p=getPrmDflt(varargin,dfs,1); resDir=p.resDir; gtDir=p.gtDir;
evalDir=[resDir '-eval/']; if(~exist(evalDir,'dir')), mkdir(evalDir); end

% check if results already exist, if so load and return
fNm = fullfile(evalDir,'eval_bdry.txt');
if(exist(fNm,'file')), R=dlmread(fNm); R=mat2cell2(R,[1 8]);
  varargout=R([4 3 2 1 7 6 5 8]); if(nargout<=8), return; end;
  R=dlmread(fullfile(evalDir,'eval_bdry_thr.txt')); P=R(:,3); R=R(:,2);
  [~,o]=unique(P); R50=interp1(P(o),R(o),max(P(o(1)),.5));
  varargout=[varargout R50]; return;
end

% perform evaluation on each image (this part can be very slow)
assert(exist(resDir,'dir')==7); assert(exist(gtDir,'dir')==7);
if(isempty(p.ids))
  ids=dir(fullfile(gtDir,'*.mat')); ids={ids.name}; n=length(ids);
  for i=1:n, ids{i}=ids{i}(1:end-4); end
else
  ids = p.ids; n = length(p.ids);
end

do=false(1,n); jobs=cell(1,n); res=cell(1,n);
for i=1:n, id = ids{i};
  res{i}=fullfile(evalDir,[id '_ev1.txt']); do(i)=~exist(res{i},'file');
  im1=fullfile(resDir,[id '.png']); gt1=fullfile(gtDir,[id '.mat']);
  jobs{i}={im1,gt1,'out',res{i},'thrs',p.thrs,'maxDist',p.maxDist,...
    'thin',p.thin}; if(0), edgesEvalImg(jobs{i}{:}); end
end
fevalDistr('edgesEvalImg',jobs(do),p.pDistr{:});

% collect evaluation results
I=dlmread(res{1}); T=I(:,1);
Z=zeros(numel(T),1); cntR=Z; sumR=Z; cntP=Z; sumP=Z;
oisCntR=0; oisSumR=0; oisCntP=0; oisSumP=0; scores=zeros(n,5);
for i=1:n
  % load image results and accumulate
  I = dlmread(res{i});
  cntR1=I(:,2); cntR=cntR+cntR1; sumR1=I(:,3); sumR=sumR+sumR1;
  cntP1=I(:,4); cntP=cntP+cntP1; sumP1=I(:,5); sumP=sumP+sumP1;
  % compute OIS scores for image
  [R,P,F] = computeRPF(cntR1,sumR1,cntP1,sumP1); [~,k]=max(F);
  [oisR1,oisP1,oisF1,oisT1] = findBestRPF(T,R,P);
  scores(i,:) = [i oisT1 oisR1 oisP1 oisF1];
  % oisCnt/Sum will be used to compute dataset OIS scores
  oisCntR=oisCntR+cntR1(k); oisSumR=oisSumR+sumR1(k);
  oisCntP=oisCntP+cntP1(k); oisSumP=oisSumP+sumP1(k);
end

% compute ODS R/P/F and OIS R/P/F
[R,P,F] = computeRPF(cntR,sumR,cntP,sumP);
[odsR,odsP,odsF,odsT] = findBestRPF(T,R,P);
[oisR,oisP,oisF] = computeRPF(oisCntR,oisSumR,oisCntP,oisSumP);

% compute AP/R50 (interpolating 100 values, has minor bug: should be /101)
if(0), R=[0; R; 1]; P=[1; P; 0]; F=[0; F; 0]; T=[1; T; 0]; end
[~,k]=unique(R); k=k(end:-1:1); R=R(k); P=P(k); T=T(k); F=F(k); AP=0;
if(numel(R)>1), AP=interp1(R,P,0:.01:1); AP=sum(AP(~isnan(AP)))/100; end
[~,o]=unique(P); R50=interp1(P(o),R(o),max(P(o(1)),.5));

% write results to disk
varargout = {odsF,odsP,odsR,odsT,oisF,oisP,oisR,AP,R50};
writeRes(evalDir,'eval_bdry_img.txt',scores);
writeRes(evalDir,'eval_bdry_thr.txt',[T R P F]);
writeRes(evalDir,'eval_bdry.txt',[varargout{[4 3 2 1 7 6 5 8]}]);

% optionally perform cleanup
if( p.cleanup ), delete([evalDir '/*_ev1.txt']);
  delete([resDir '/*.png']); rmdir(resDir); end

end

function [R,P,F] = computeRPF(cntR,sumR,cntP,sumP)
% compute precision, recall and F measure given cnts and sums
R=cntR./max(eps,sumR); P=cntP./max(eps,sumP); F=2*P.*R./max(eps,P+R);
end

function [bstR,bstP,bstF,bstT] = findBestRPF(T,R,P)
% linearly interpolate to find best thr for optimizing F
if(numel(T)==1), bstT=T; bstR=R; bstP=P;
  bstF=2*P.*R./max(eps,P+R); return; end
A=linspace(0,1,100); B=1-A; bstF=-1;
for j = 2:numel(T)
  Rj=R(j).*A+R(j-1).*B; Pj=P(j).*A+P(j-1).*B; Tj=T(j).*A+T(j-1).*B;
  Fj=2.*Pj.*Rj./max(eps,Pj+Rj); [f,k]=max(Fj);
  if(f>bstF), bstT=Tj(k); bstR=Rj(k); bstP=Pj(k); bstF=f; end
end
end

function writeRes( alg, fNm, vals )
% write results to disk
k=size(vals,2); fNm=fullfile(alg,fNm); fid=fopen(fNm,'w');
if(fid==-1), error('Could not open file %s for writing.',fNm); end
frmt=repmat('%10g ',[1 k]); frmt=[frmt(1:end-1) '\n'];
fprintf(fid,frmt,vals'); fclose(fid);
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function obj = addAttribute(obj, attribute)
% ADDATTRIBUTE Add an attribute to the object.
% ------------------------------------------------------------------------------
% INPUT
% 1 [attribute]
%   Name of attribute which should be added. Actually, possible choices are:
%     * 'slope'
%        Slope (from 0 to 100 gradian) of each point based on its normal vector.
%        This attribute requires the components of the normal vectors nx, ny,
%        and nz for each point. They can be calculated with the method 'normals'
%        (see 'help pointCloud.normals').
%     * 'exposition'
%        Exposition (from 0 to 400 gradian) of each point based on its normal
%        vector. This attribute requires the components of the normal vectors
%        nx, ny, and nz for each point. They can be calculated with the method
%        'normals' (see 'help pointCloud.normals').
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 Calculate and visualize the slope attribute of a point cloud.
%   pc = pointCloud('Lion.xyz', 'Attributes', {'nx' 'ny' 'nz' 'roughness'});
%   pc.addAttribute('slope');
%   pc.plot('Color', 'A.slope', 'MarkerSize', 5);
%
% 2 Calculate and visualize the exposition attribute of a point cloud.
%   pc = pointCloud('Lion.xyz', 'Attributes', {'nx' 'ny' 'nz' 'roughness'});
%   pc.addAttribute('exposition');
%   pc.plot('Color', 'A.exposition', 'MarkerSize', 5);
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

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

validAttribute = {'slope' 'exposition'};

p = inputParser;
p.addRequired('attribute', @(x) any(strcmpi(x, validAttribute)));
p.parse(attribute);
p = p.Results;
% Clear required input to avoid confusion
clear attribute

% Start ------------------------------------------------------------------------

procHierarchy = {'POINTCLOUD' 'ADDATTRIBUTE'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Point cloud label = ''%s''', obj.label));

% Add attribute ----------------------------------------------------------------

% Slope
if strcmpi(p.attribute, 'slope')

    % Points with no normal (i.e. with normal components equal to NaN)
    idxNoNormal = isnan(obj.A.nx) | isnan(obj.A.ny) | isnan(obj.A.nz);

    obj.A.slope = NaN(size(obj.X,1),1); % initialize with NaNs
    obj.A.slope(~idxNoNormal) = acosg(abs(obj.A.nz(~idxNoNormal))); % attention: abs introduced to make sure, slope ranges between 0 and 100, but in TLS slope > 100 may be useful

% Exposition
elseif strcmpi(p.attribute, 'exposition')

    % Points with no normal (i.e. with normal components equal to NaN)
    idxNoNormal = isnan(obj.A.nx) | isnan(obj.A.ny) | isnan(obj.A.nz);

    polar = xyz2polar([obj.A.nx(~idxNoNormal) obj.A.ny(~idxNoNormal) obj.A.nz(~idxNoNormal)]);
    obj.A.exposition = NaN(size(obj.X,1),1); % initialize with NaNs
    obj.A.exposition(~idxNoNormal) = polar(:,2);

    % Set exposition to zero if normal points upwards (otherwise it is equal to NaN)
    obj.A.exposition(obj.A.nx == 0 & obj.A.ny == 0 & obj.A.nz == 1) = 0;

end

% End --------------------------------------------------------------------------

msg('E', procHierarchy);

end

author: pglira


text:
function Im=RGB_YUV(Y,U,V)
% This program transform YUV layers to RGB Layers in ome matrix 'Im'....
%  By   Mohammed Mustafa Siddeq
%  Date 25/7/2010

G=round((Y-(U+V)/4));
R=U+G;
B=V+G;
Im(:,:,1)=R; Im(:,:,2)=G; Im(:,:,3)=B;
%imshow(uint8(Im));
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = tps(varargin)
% VL_TPS  Compute the thin-plate spline basis
%   PHI=VL_TPS(X1,X2,Y) returns the basis PHI of a thin-plate spline
%   (TPS) defined on the domain X1,X2 with control points Y.
%
%   X1 and X2 are MxN matrices specifying the grid vertices.  When
%   warping images, these usually correspond to image pixels.
%
%   Y is a 2xK matrix specifying the control points, one per
%   column. Ofthen Y is a subset of the domain X1,X2, but this is not
%   required.
%
%   PHI is a (K+3)xNxM matrix, with one layer per basis element. Each
%   basis element is a function of the domain X1,X2.
%
%   [PHI,S] = VL_TPS(X1,X2,Y) additionally returns the stiffness matrix S
%   of the TPS.
%
%   See also: VL_WTPS(), VL_HELP().
[varargout{1:nargout}] = vl_tps(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
classdef Plotter < handle
% PLOTTER supports plotting data in 2D with a reference signal.
%
% Properties:
%   r           - Reference signal
%   x           - Time
%   y           - Output signal
%
% Methods:
%   plot_2d_ref - Plots an output and reference signal over time.

    properties
    %% PROPERTIES

        t   % Time
        y   % Output signal
        r   % Reference signal

        h
        g
        a
    end

    methods
        function obj = Plotter()
        % PLOTTER Constructor

            obj.t = 0;
            obj.y = 0;
            obj.r = 0;
            obj.h = -1;
            obj.g = -1;

            figure;
            obj.a = axes;
            set(obj.a, 'NextPlot', 'add');
            hold(obj.a, 'all');
            obj.t = 0;
        end

        function plot_2d_ref(obj, dt, y, r, color)
        %% PLOT_2D_REF Plots an output and reference signal over time
        %   [h,g] = plot_2d_ref(obj, h, g, x, y, r) plots the output signal
        %   (y) and reference signal (r) versus time (t).

            if ~ishandle(obj.h)
                obj.h = plot(obj.a, dt, y, 'b');
                obj.g = plot(obj.a, dt, r, '--');
                set(obj.g, 'Color', color);
                obj.t = obj.t(end);
                obj.y = obj.y(end);
                obj.r = obj.r(end);
            end

            obj.t = [obj.t obj.t(end)+dt];
            obj.y = [obj.y y];
            obj.r = [obj.r r];

            set(obj.h, 'XData', obj.t);
            set(obj.h, 'YData', obj.y);
            set(obj.g, 'XData', obj.t);
            set(obj.g, 'YData', obj.r);

        end

        function switch_2d_ref(obj)
            obj.h = -1;
        end

    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function dtiMrMeshSave(handles)
% Save the dti mesh structure and data required for mrmViewer.
%
%   dtiMrMeshSave(handles)
%
% The mrmViewer code should be:
%
%   load Test
%   msh = dtiGet(handles,'mrMesh');
%   msh.id = 101;   % Set up a number for the window
%   mrmSet(msh,'refresh');
%   msh = dtiInitMrMeshWindow(msh);
%   dtiMrMeshAddROIs(handles,msh);
%   dtiMrMeshAddFGs(handles,msh);
%   dtiMrMeshAddImages(handles,msh,origin,xIm,yIm,zIm);
%
% %Author: Wandell
% (c) Stanford VISTA Team

%   Programming notes.
%   Should we write a dtiLoadMesh, too?  That is not necessary, really,
%   because we should always be able to display from the current
%   data and ....
%

if notDefined('handles'), error('dtiFiberUI handles required.'); end

% Figure out the file
persistent meshPath;
if(isempty(meshPath)), meshPath = fullfile(handles.defaultPath, 'slice'); end
[f, p] = uiputfile({'*.mat';'*.*'}, 'Save current mesh...', meshPath);
if(isnumeric(f)), disp('dtiSaveMesh cancelled.'); return; end
fname = fullfile(p, f);
meshPath = fname;

% Get all of the potential images
anat = dtiGet(handles,'currentanatomydata');
curPosition  = dtiGet(handles,'curposition');
curXform = dtiGet(handles,'curacpcxform');
zIm = dtiGetSlice(curXform, anat, 3, curPosition(3));
yIm = dtiGetSlice(curXform, anat, 2, curPosition(2));
xIm = dtiGetSlice(curXform, anat, 1, curPosition(1));

% Empty out the images not needed for display
[xIm,yIm,zIm] = dtiMrMeshSelectImages(handles,xIm,yIm,zIm); %#ok<NASGU>
origin = dtiGet(handles,'origin'); %#ok<NASGU>

% Make the data size manageable by getting rid of the big variables
handles = rmfield(handles,'dt6');
%handles = rmfield(handles,'vec');
handles = rmfield(handles,'bg');

if isfield(handles,'t1NormParams')
    handles = rmfield(handles,'t1NormParams');
end
handles = rmfield(handles,'brainMask'); %#ok<NASGU>

% Save the file for mrmViewer
save(fname,'handles','origin','xIm','yIm','zIm');

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function    g = lfmGradientH31(preFactor, preFactorGrad, gradThetaGamma, ...
    gradUpsilon1, gradUpsilon2, compUpsilon1, compUpsilon2, mode, term)

% LFMGRADIENTH31 Gradient of the function h_i(z) with respect to some of the
% hyperparameters of the kernel: m_k, C_k, D_k, m_r, C_r or D_r.
% FORMAT
% DESC Computes the gradient of the function h_i(z) with respect to some of
% the parameters of the system (mass, spring or damper).
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG gradThetaGamma : Vector with the gradient of gamma1 and gamma2 with
% respect to the desired parameter.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1)
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN g : Gradient of the function with respect to the desired
% parameter.
%
% COPYRIGHT : David Luengo, 2007, 2008,
%
% COPYRIGHT : Mauricio Alvarez, 2008
%
% SEEALSO : lfmKernGradient, lfmXlfmKernGradient, lfmGradientUpsilon

% KERN

% Gradient evaluation

if nargin<9
    term =[];
end

if ~mode
    if ~term
        g = (preFactor(1)*gradUpsilon1 + preFactorGrad(1)*compUpsilon1)*gradThetaGamma;
    else
        g = (-preFactor(1)*gradUpsilon1 + preFactorGrad(1)*compUpsilon1)*gradThetaGamma(1) + ...
           (preFactor(2)*conj(gradUpsilon1) - preFactorGrad(2)*conj(compUpsilon1))*gradThetaGamma(2);
    end
else
    g = (preFactor(1)*gradUpsilon1 + preFactorGrad(1)*compUpsilon1)*gradThetaGamma(1) +...
        (preFactor(2)*gradUpsilon2 + preFactorGrad(2)*compUpsilon2)*gradThetaGamma(2);
    %g = (preFactor(2)*gradUpsilon2 + preFactorGrad(2)*compUpsilon2)*gradThetaGamma(2);
    %g = (gradUpsilon2)*gradThetaGamma(2);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ A, T, error, alignedShape ] = AlignShapesWithScale( alignFrom, alignTo )
%ALIGNSHAPESWITHSCALE Summary of this function goes here
%   Detailed explanation goes here

    numPoints = size(alignFrom,1);

    meanFrom = mean(alignFrom);
    meanTo = mean(alignTo);

    alignFromMeanNormed = bsxfun(@minus, alignFrom, meanFrom);
    alignToMeanNormed = bsxfun(@minus, alignTo, meanTo);

    % scale now
    sFrom = sqrt(sum(alignFromMeanNormed(:).^2)/numPoints);
    sTo = sqrt(sum(alignToMeanNormed(:).^2)/numPoints);

    s = sTo / sFrom;

    alignFromMeanNormed = alignFromMeanNormed/sFrom;
    alignToMeanNormed = alignToMeanNormed/sTo;

    [R, t] = AlignShapesKabsch(alignFromMeanNormed, alignToMeanNormed);

    A = s * R;
    aligned = (A * alignFrom')';
    T = mean(alignTo - aligned);
    alignedShape = bsxfun(@plus, aligned, T);
    error = sqrt(mean(sum((alignedShape - alignTo).^2,2)));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function y = evaluate_log(this, x)
%EVALUATE_LOG   Evaluate robust function.
% opposite to evaluate, created to be consistent with the class expert

%   EVALUATE(THIS, X) evaluates the robust function THIS at point(s) X.
%
%   This is a member function of the class 'robust_function'.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date:  $
%   $Revision: $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  y = -feval(this.type, x', this.param, 0);

author: kristinbranson

repo: JAABA


text:
% STD_UNIFORMSETINDS - Check uniform channel distribution across datasets
%
% Usage:
%                >> boolval = std_uniformsetinds(STUDY);
% Inputs:
%   STUDY      - EEGLAB STUDY
%
% Outputs:
%   boolval    - [0|1] 1 if uniform
%
% Authors: Arnaud Delorme, SCCN/UCSD, CERCO/CNRS, 2010-

% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function uniformchannels = std_uniformsetinds( STUDY );

uniformchannels = 1;
for c = 1:length(STUDY.changrp(1).setinds(:))
    tmpind = cellfun(@(x)(length(x{c})), { STUDY.changrp(:).setinds });
    if length(unique(tmpind)) ~= 1
        uniformchannels = 0;
    end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%VL_NNBNORM CNN batch normalisation.
%   Y = VL_NNBNORM(X,G,B) applies batch normalization to the input
%   X. Batch normalization is defined as:
%
%      Y(i,j,k,t) = G(k) * (X(i,j,k,t) - mu(k)) / sigma(k) + B(k)
%
%   where:
%
%      mu(k) = mean_ijt X(i,j,k,t),
%      sigma2(k) = mean_ijt (X(i,j,k,t) - mu(k))^2,
%      sigma(k) = sqrt(sigma2(k) + EPSILON)
%
%   are respectively the per-channel mean, variance, and standard
%   deviation of each feature channel in the data X. The parameters
%   G(k) and B(k) are multiplicative and additive constants use to
%   scale each data channel.
%
%   Means and variances are accumulated across all the data items
%   (images) stored in the 4D tensor X (from which the name batch
%   normalization is derived). The constant EPSILON is used to
%   regularize the computation of sigma(k) and to avoid division by
%   zero.
%
%   [DZDX,DZDG,DZDB] = VL_NNBNORM(X,G,B,DZDY) computes the derviatives
%   of the block projected onto DZDY. DZDX, DZDG, DZDB and DZDY have
%   the same dimensions as X, G, B, and Y respectivey.
%
%   Optionally, [Y,MOMENTS] = VL_NNBNORM(...) and
%   [DZDX,DZDG,DZDB,MOMENTS] = VL_NNBNORM(...,DZDY) return the values
%   of the vectors mu and sigma in the formulas above. Here, MOMENTS
%   is a DEPTH x 2 array [MU, SIGMA].
%
%   VL_NNBNROM(..., 'Option', value) takes the following options:
%
%   `Epsilon`:: 1e-4
%       Specifies the constant EPSILON in the formuals above.
%
%   `Moments`:: unspecified
%       Specifies an array MOMENTS with the values of mu and sigma to
%       use instead of computing them according to the equations
%       above. This is useful to disable batch normalization during
%       testing.
%
%   `CuDNN`:: specified
%       If specified, turns on CuDNN. CuDNN is on by default. This
%       option can be useful to undo the effect of a previous
%       `NoCuDNN` option in the argument list.
%
%   `NoCuDNN`:: not specified
%       If specified, turns off CuDNN.
%
%   See also: VL_NNNORMALIZE().

% Copyright (C) 2015 Sébastien Ehrhardt, Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function obj = stat_violin(obj,varargin)
% stat_violin display violin plots of y data for unique x
%
% Options can be given as 'name',value pairs:
% - all options from see stat_density() are present:
% 'bandwidth', 'kernel', 'npoints', 'extra_y'
% - 'normalization': 'area' (default) all violins have the same
% areas. With 'count', the violin areas are proportional to the number of
% points in each category. With 'width', the maximum widths of all violins
% are matched
% - 'half': If set to true, only half violins are drawn, and the color
% index determines if the left half or right hald is drawn. Useful for
% comparing distributions across two colors.
% - 'fill' see stat_bin()
% - 'width' and 'dodge' see stat_summary()
%
% See also stat_summary(), stat_bin(), stat_density()

p=inputParser;
my_addParameter(p,'bandwidth',-1);
my_addParameter(p,'kernel','normal')
my_addParameter(p,'npoints',100)
my_addParameter(p,'extra_y',0)
my_addParameter(p,'normalization','area') %'count' 'width'
my_addParameter(p,'fill','face')
my_addParameter(p,'width',0.6)
my_addParameter(p,'dodge',0.7)
my_addParameter(p,'half',false)
parse(p,varargin{:});

obj.geom=vertcat(obj.geom,{@(dobj,dd)my_violin(dobj,dd,p.Results)});
obj.results.stat_violin={};

end

function hndl=my_violin(obj,draw_data,params)

x=comb(draw_data.x);
y=comb(draw_data.y);

uni_x=unique(x);
uni_x(diff(uni_x)<1e-10)=[];

dens=cell(length(uni_x),1);
dens_pos=cell(length(uni_x),1);

if params.half
    params.dodge=0;
end

if params.dodge>0
    boxw=draw_data.dodge_avl_w*params.width./(draw_data.n_colors);
else
    boxw=draw_data.dodge_avl_w*params.width;
end
boxmid=dodger(uni_x,draw_data,params.dodge);

[face_color , face_alpha , edge_color , ~] = parse_fill (params.fill,draw_data.color);
temp_line_params={'LineStyle',draw_data.line_style,'Color',edge_color,'lineWidth',1};

%Maximum area
max_area=draw_data.dodge_avl_w*(obj.var_lim.maxy-obj.var_lim.miny);

if obj.handle_graphics
    lines=gobjects(1,length(uni_x));
    patches=gobjects(1,length(uni_x));
else
    lines=zeros(1,length(uni_x));
    patches=zeros(1,length(uni_x));
end

for ind_x=1:length(uni_x)
    %And here we have a loose selection also because of
    %potential numerical errors
    ysel=y(abs(x-uni_x(ind_x))<1e-10);

    if ~isempty(ysel)

        extra_y=(max(ysel)-min(ysel))*params.extra_y;
        binranges=linspace(min(ysel)-extra_y,max(ysel)+extra_y,params.npoints);

        if params.bandwidth>0
            [dens{ind_x},dens_pos{ind_x}] = ksdensity(ysel,binranges,'function','pdf',...
                'bandwidth',params.bandwidth,'kernel',params.kernel);
        else
            [dens{ind_x},dens_pos{ind_x}] = ksdensity(ysel,binranges,'function','pdf',...
                'kernel',params.kernel);
        end

        area=sum(dens{ind_x})*2*(binranges(2)-binranges(1));

        switch params.normalization
            case 'area'
                %Normalized area is a third of the max available area
                dens{ind_x}=dens{ind_x}*max_area/(3*area);
            case 'count'
                %Normalized count makes area/max available area correspond to
                %n points/n total points
                dens{ind_x}=dens{ind_x}*length(ysel)*max_area/(obj.data_size*area);
            case 'width'
                %Makes maximum widths equal
                dens{ind_x}=0.5*dens{ind_x}/max(dens{ind_x});
        end

        %Adjust width
        dens{ind_x}=dens{ind_x}*boxw;

        if params.half %If in half mode
            if ~mod(draw_data.color_index,2) %We draw right half if even index
                xpatch=[boxmid(ind_x)-dens{ind_x}(1:end-1) ; ...
                    boxmid(ind_x)+zeros(1,length(dens{ind_x})-1) ; ...
                    boxmid(ind_x)+zeros(1,length(dens{ind_x})-1);  ...
                    boxmid(ind_x)-dens{ind_x}(2:end)];
            else %left half if odd
                xpatch=[boxmid(ind_x)+zeros(1,length(dens{ind_x})-1) ; ...
                    boxmid(ind_x)+dens{ind_x}(1:end-1) ; ...
                    boxmid(ind_x)+dens{ind_x}(2:end) ; ...
                    boxmid(ind_x)+zeros(1,length(dens{ind_x})-1)];

            end

        else %Otherwise draw full violin
            xpatch=[boxmid(ind_x)-dens{ind_x}(1:end-1) ; ...
                boxmid(ind_x)+dens{ind_x}(1:end-1) ; ...
                boxmid(ind_x)+dens{ind_x}(2:end);  ...
                boxmid(ind_x)-dens{ind_x}(2:end)];
        end

        ypatch=[dens_pos{ind_x}(1:end-1) ; ...
            dens_pos{ind_x}(1:end-1) ; ...
            dens_pos{ind_x}(2:end) ; ...
            dens_pos{ind_x}(2:end)];

        %Draw patch
        patches(ind_x)=patch(xpatch,ypatch,[1 1 1],'EdgeAlpha',1,'FaceColor',face_color,'EdgeColor','none','FaceAlpha',face_alpha);

        %Draw lines
        if params.half %If in half mode
            if ~mod(draw_data.color_index,2) %We draw right half if even index
                lines(ind_x)=line(boxmid(ind_x)-[0 dens{ind_x} 0 0] , ...
                    [dens_pos{ind_x}(1) dens_pos{ind_x} dens_pos{ind_x}(end) dens_pos{ind_x}(1)],temp_line_params{:});

            else %left half if odd
                lines(ind_x)=line(boxmid(ind_x)+[0 dens{ind_x} 0 0] , ...
                    [dens_pos{ind_x}(1) dens_pos{ind_x} dens_pos{ind_x}(end) dens_pos{ind_x}(1)],temp_line_params{:});
            end

        else %Otherwise draw full violin
                lines(ind_x)=line([boxmid(ind_x) dens{ind_x}+boxmid(ind_x) boxmid(ind_x) boxmid(ind_x)-fliplr(dens{ind_x}) boxmid(ind_x)] , ...
                 [dens_pos{ind_x}(1) dens_pos{ind_x} dens_pos{ind_x}(end) fliplr(dens_pos{ind_x}) dens_pos{ind_x}(1)],temp_line_params{:});
        end

    end

end

obj.results.stat_violin{obj.result_ind,1}.densities=dens;
obj.results.stat_violin{obj.result_ind,1}.densities_y=dens_pos;
obj.results.stat_violin{obj.result_ind,1}.unique_x=uni_x;
obj.results.stat_violin{obj.result_ind,1}.line_handle=lines;
obj.results.stat_violin{obj.result_ind,1}.fill_handle=patches;

end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function [b0vector b0angle] = getB0direction(vw)
% Return the direction of the b0 field from a scan as (1) a unit vector and
% (2) an angle in degrees
%
%   [b0vector b0angle]  = getB0direction([vw])
%
% Examples:
%   b0vector = getB0direction(vw)
%   [b0vector b0angle] = getB0direction
%
% vw can be INPLANE, GRAY, or VOLUME

mrGlobals

if ~exist('vw', 'var'), vw = getCurView; end

% get the scanner xform (image coords => mm coords)
scannerXform = viewGet(vw, 'scannerXform');

% Define a unit vector in the z-direction in scanner space. This is the B0
% direction (by definition).
origin = [0 0 0]; z = [0 0 1];

% Transform the vector into INPLANE image space
b0ip = scannerXform \ [origin 1; z 1]';
b0ip = b0ip(1:3,2) - b0ip(1:3,1);

% Make it a unit vector
b0ip = 1/norm(b0ip) * b0ip;

% xform the b0 vector to the appropriate view type
switch lower(viewGet(vw, 'viewType'))
    case 'inplane'
        % done
        b0vector = b0ip;
    case {'gray' 'volume'}
        xform = mrSESSION.alignment;
        b0Vol = xform * [[0 0 0 1]' [b0ip; 1]];
        % make it a vector of norm 1
        b0vector = b0Vol(1:3,2) - b0Vol(1:3,1);
        b0vector = 1/norm(b0vector) * b0vector;
        % xform to xyz
        b0vector = [-b0vector(3) b0vector(2) -b0vector(1)];
    otherwise
        error('Can''t get B0 direction in %s vw', viewGet(vw, 'viewType'))
end

% convert b0 vector into angle. This is the angle between z-vector in image
% space and and z-vector in scanner space.
costheta    = dot(b0vector, z)/norm(b0vector)*norm(z);
theta       = acos(costheta);
b0angle     = theta  / pi * 180;


text:
function fn_test ( )

%*****************************************************************************80
%
%% FN_TEST tests the FN library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( )

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FN_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the FN library.\n' );

  acos_test ( );
  acosh_test ( );
  ai_test ( );
  aid_test ( );
  asin_test ( );
  asinh_test ( );
  atan_test ( );
  atan2_test ( );
  atanh_test ( );
  besi0_test ( );
  besi1_test ( );
  besj0_test ( );
  besj1_test ( );
  besk_test ( );
  besk0_test ( );
  besk1_test ( );
  besy0_test ( );
  besy1_test ( );
  beta_test ( );
  betai_test ( );
  bi_test ( );
  bid_test ( );
  binom_test ( );
  cbrt_test ( );
  chi_test ( );
  chu_test ( );
  ci_test ( );
  cin_test ( );
  cinh_test ( );
  cos_test ( );
  cos_deg_test ( );
  cosh_test ( );
  cot_test ( );
  dawson_test ( );
  e1_test ( );
  ei_test ( );
  erf_test ( );
  erfc_test ( );
  exp_test ( );
  fac_test ( );
  gamma_test ( );
  gamma_inc_test ( );
  gamma_inc_tricomi_test ( );
  int_test ( );
  lbeta_test ( );
  li_test ( );
  lngam_test ( );
  log_test ( );
  log10_test ( );
  poch_test ( );
  psi_test ( );
  rand_test ( );
  shi_test ( );
  si_test ( );
  sin_test ( );
  sin_deg_test ( );
  sinh_test ( );
  spence_test ( );
  sqrt_test ( );
  tan_test ( );
  tanh_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FN_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ correlations, rms, mean_correlation, mean_RMSE, long_correlation, long_RMSE, predictions, gts ] = evaluate_CCNF_model( alphas, betas, thetas, x, y, similarityFNs, sparsityFNs, offset, scaling, verbose, PrecalcQ2sFlat)
%evaluate_CCNF_model Evaluate the trained model on test (or training data)

% For visualising time series predictions
num_x_plots = 8;
num_y_plots = 10;

total_plots = num_x_plots * num_y_plots;

if(iscell(x))
    num_seqs = numel(x);
    x = cell2mat(x)';
    % add a bias term
    x =  cat(1, ones(1,size(x,2)), x);
else
    % if not a cell it has already been flattened, and is constant
    % (most likely)
    num_seqs = size(y,2);

end

% if not sure about const assume it is not
const = false;

if(nargin < 11)
    [ ~, ~, PrecalcQ2sFlat, ~ ] = CalculateSimilarities( num_seqs, x, similarityFNs, sparsityFNs, y, const);
end

correlations = zeros(num_seqs, 1);
rms = zeros(num_seqs, 1);

% concatenated data for an alternative correlation
y_predConcat = [];
y_trueConcat = [];

% Predict each sequence
for q=1:num_seqs

    if(iscell(y))
        seq_length = size(y{q},1);
        yq = y{q};
    else
        seq_length = size(y,1);
        yq = y(:,q);
    end

    X = x(:,(q-1)*seq_length+1:q*seq_length);

    h1 = 1./(1 + exp(-thetas * X));
    b = (2 * alphas' * h1)';

    PrecalcQ2flat = PrecalcQ2sFlat{q};

    precalc_eye = eye(seq_length);
    precalc_zeros = zeros(seq_length);

    SigmaInv = CalcSigmaCCNFflat(alphas, betas, seq_length, PrecalcQ2flat, precalc_eye, precalc_zeros);

    y_est = SigmaInv \ b;

    % Can optionally supply the scaling and offset used on the training
    % labels to be applied inversely
    y_est = y_est/scaling + offset;

    if(numel(y_est) > 1)
        R = corrcoef(y_est, yq);
        correlations(q) = R(1,2);
    end

    rms(q) = sqrt( mean((y_est - yq).^2) );

    y_predConcat = cat(1, y_predConcat, y_est);
    y_trueConcat = cat(1, y_trueConcat, yq);

    if(verbose)

        if(mod(q,total_plots) == 1)
            figure;
            remainingPlots = nExamples - q;
            if(remainingPlots < total_plots)
                num_y_plots = ceil(remainingPlots / num_x_plots);
            end
        end

        subplot(num_y_plots,num_x_plots,mod(q-1,total_plots)+1);
        t = 1:nFrames;
        plot(t,y{q},'g',t,y_est,'b');
        title(sprintf('C %.2f, R %.2f', correlations(q), rms(q)));
        set(gca, 'XTick', [], 'YTick', []);

    end

end

% Compute the error metrics
mean_correlation = mean(correlations);
mean_RMSE = mean(rms);
long_correlation = corr(y_predConcat, y_trueConcat).^2;

long_RMSE = sqrt(mean((y_predConcat - y_trueConcat).^2));
predictions = y_predConcat;
gts = y_trueConcat;

if(verbose)
    figure
    plot([1:numel(y_trueConcat)],y_trueConcat,'g',[1:numel(y_trueConcat)],y_predConcat,'b');
    title(sprintf('C %.2f, R %.2f', long_correlation, long_RMSE));
    set(gca, 'XTick', [], 'YTick', []);
end

end


text:
function f=idtwfbreal(c,par,varargin)
%IDTWFBREAL Inverse Dual-tree Filterbank for real-valued signals
%   Usage:  f=idtwfbreal(c,info);
%           f=idtwfbreal(c,dualwt,Ls);
%
%   Input parameters:
%         c           : Input coefficients.
%         info        : Transform params. struct
%         dualwt      : Dual-tree Wavelet Filterbank definition
%         Ls          : Length of the reconstructed signal.
%
%   Output parameters:
%         f     : Reconstructed data.
%
%   `f = idtwfbreal(c,info)` reconstructs real-valued signal *f* from the
%   coefficients *c* using parameters from `info` struct. both returned by
%   |dtwfbreal| function.
%
%   `f = idtwfbreal(c,dualwt,Ls)` reconstructs real-valued signal *f* from the
%   coefficients *c* using dual-tree filterbank defined by `dualwt`. Plese
%   see |dtwfbreal| for supported formats. The *Ls* parameter is mandatory
%   due to the ambiguity of reconstruction lengths introduced by the
%   subsampling operation.
%   Note that the same flag as in the |dtwfbreal| function have to be used,
%   otherwise perfect reconstruction cannot be obtained. Please see help
%   for |dtwfbreal| for description of the flags.
%
%   Examples:
%   ---------
%
%   A simple example showing perfect reconstruction using |idtwfbreal|:::
%
%      f = gspi;
%      J = 7;
%      wtdef = {'qshift3',J};
%      c = dtwfbreal(f,wtdef);
%      fhat = idtwfbreal(c,wtdef,length(f));
%      % The following should give (almost) zero
%      norm(f-fhat)
%
%   See also: dtwfbreal dtwfbinit

complainif_notenoughargs(nargin,2,'IDTWFBREAL');

if(~iscell(c))
   error('%s: Unrecognized coefficient format.',upper(mfilename));
end

if(isstruct(par)&&isfield(par,'fname'))
   complainif_toomanyargs(nargin,2,'IDTWFBREAL');

   if ~strcmpi(par.fname,'dtwfbreal')
      error(['%s: Wrong func name in info struct. ',...
             ' The info parameter was created by %s.'],...
             upper(mfilename),par.fname);
   end

   dtw = dtwfbinit({'dual',par.wt},par.fOrder);
   Ls = par.Ls;
   ext = 'per';
   L = wfbtlength(Ls,dtw,ext);
else
   complainif_notenoughargs(nargin,3,'IDTWFBREAL');

   %% PARSE INPUT
   definput.keyvals.Ls=[];
   definput.keyvals.dim=1;
   definput.import = {'wfbtcommon'};

   [flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);
   complainif_notposint(Ls,'Ls');

   ext = 'per';
   % Initialize the wavelet tree structure
   dtw = dtwfbinit(par,flags.forder);

   [Lc,L]=wfbtclength(Ls,dtw,ext);

   % Do a sanity check
   if ~isequal(Lc,cellfun(@(cEl) size(cEl,1),c))
      error(['%s: The coefficient subband lengths do not comply with the'...
             ' signal length *Ls*.'],upper(mfilename));
   end
end

%% ----- step 3 : Run computation
[nodesBF, rangeLoc, rangeOut] = treeBFranges(dtw,'rev');
outLengths = nodesInLen(nodesBF,L,strcmpi(ext,'per'),dtw);
outLengths(end) = L;

f = comp_idtwfb(c,dtw.nodes(nodesBF),dtw.dualnodes(nodesBF),outLengths,...
                rangeLoc,rangeOut,ext,0);
f = postpad(f,Ls);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function test_suite = test_nndist
%TEST_NNDIST  One-line description here, please.
%
%   output = test_nndist(input)
%
%   Example
%   test_nndist
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2016-07-17,    using Matlab 8.6.0.267246 (R2015b)
% Copyright 2016 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testTwoPoints(testCase) %#ok<*DEFNU>
% compute for two points

p1 = [10 10];
p2 = [20 10];
pts = [p1 ; p2];

[dists, inds] = nndist(pts);

testCase.assertEqual([10; 10], dists, 'AbsTol', .01);
testCase.assertEqual([2;1], inds, 'AbsTol', .01);

function testFourPoints(testCase) %#ok<*DEFNU>
% compute for four points

p1 = [10 10];
p2 = [20 10];
p3 = [20 30];
p4 = [50 30];
pts = [p1 ; p2 ; p3 ; p4];

[dists, inds] = nndist(pts);

testCase.assertEqual([10; 10; 20; 30], dists, 'AbsTol', .01);
testCase.assertEqual([2;1;2;3], inds, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function a = ne( x, y )

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

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function y = nigcdf(x, alpha, beta, mu, delta)
%NIGCDF Normal-Inverse-Gaussian cumulative distribution function (cdf).
%   Y = NIGCDF(X, ALPHA, BETA, MU, DELTA) returns the cdf of the
%   Normal-Inverse-Gaussian distribution with the parameter BETA which
%   determines the skewness, shape parameter ALPHA, location parameter MU
%   and scale parameter DELTA, evaluated at the values in X.
%
%   ALPHA, BETA, MU, DELTA must be scalar values.
%
%   ALPHA and DELTA must be positive values.
%
%   MU must be a value greater than -inf and smaller than inf.
%
%   Default values for ALPHA = 1, BETA = 0; MU = 0; DELTA = 1.
%
%   See also NIGFITC, NIGFITG, NIGFITM, NIGFITP, NIGINV, NIGMM, NIGPAR,
%            NIGPDF, NIGRND, NIGSTAT.

%   References:
%      [1] de Beus, P., Bressers, M. and de Graaf, T. (2003) Alternative
%          Investments and Risk Measurement, Appendix 2.

% -------------------------------------------------------------------------
%
% risklab germany GmbH
% Nymphenburger Strasse 112 - 116
% D-80636 Muenchen
% Germany
% Internet:    www.risklab.de
% email:       info@risklab.de
%
% Implementation Date:  2004 - 10 - 14
% Author:               Dr. Ralf Werner, Michaela Tempes
%                       ralf.werner@risklab.de
% -------------------------------------------------

% Default values
if nargin < 2
    alpha = 1;
end
if nargin < 3
    beta = 0;
end
if nargin < 4
    mu = 0;
end
if nargin < 5
    delta = 1;
end

% Constraints for the parameters
if alpha <= 0
    error('ALPHA must be positive.');
end
if delta <= 0
    error('DELTA must be positive.');
end
if ((mu == -inf) || (mu == inf))
    error('MU muss aus (-inf,inf) sein');
end
if abs(beta)>=alpha
    error('|BETA| muss kleiner gleich ALPHA sein');
end

% approximation of the integral through a summation
N = 1000;           % number of addends
ind = (1:N-1);
ind = 1./ind;
dummy = N*ind-1;

if size(x, 1) <= 1  % make x a column vector
    x = x';
end

M = length(x);

% the following is necessary due to limitations in available memory
% this accelerates computation time

maxSize = 1000;
nSteps = ceil(M / maxSize);
xLarge = x;
y = [];

for iSteps = 1:nSteps

    xSize = min(length(xLarge), maxSize);
    x = xLarge(1:xSize);

    constFact = N*delta/sqrt(2*pi)*exp(delta*sqrt(alpha^2-beta^2));
    firstFact = ind.^2;
    secondFact = dummy.^(-1.5);
    expFact = exp(-1/2*(delta^2./dummy + (alpha^2-beta^2)*dummy));

    X = repmat(x, 1, N-1);
    DUMMY = repmat(dummy, xSize, 1);
    normalFact = normcdf(X, mu + beta*DUMMY, sqrt(DUMMY));

    y1 = constFact*firstFact.*secondFact.*expFact;
    ySmall = repmat(y1, xSize, 1).*normalFact;

    ySmallSum = sum(ySmall, 2);
    y = [y; ySmallSum];

    xLarge = xLarge(xSize + 1:end);

end


text:
function obj = GEVD_method(obj)

%the proposed GEVD-HBF scheme

global H Vn  Nrf Nt Nr W_mopt;
t1 = clock;
i = 0;   %itertion index

%just to achieve W_RF*W_D = W_mopt
W_equal = W_mopt;
w = trace (W_equal' * W_equal);
%random initialization
V_RF = exp( 1i*unifrnd(0,2*pi,Nt,Nrf));
W_RF = exp( 1i*unifrnd(0,2*pi,Nr,Nrf));
%iteration trigger, the normal initialization just for pass into functions
trigger = 1;
m_MSE_new = 100;

%limit the iterations number by i<10
while (trigger > 1e-5 && i<10)

    % precoding
    H1 = H' * W_equal;
    [V_RF, V_U] = gevd_algorithm(V_RF, w, H1);
    V_equal = V_RF *V_U;
    v = trace (V_equal * V_equal');   %beta^(-2)

    %combining
    H2 = H * V_equal;
    [W_RF, W_B] = gevd_algorithm(W_RF, v, H2);
    W_equal = W_RF * W_B;
    w = trace (W_equal' * W_equal);
    %modified MSE
    H_equal = W_equal'*H2;

    m_MSE_old = m_MSE_new;
    m_MSE_new = trace(H_equal * H_equal' - H_equal - H_equal') + Vn * v * w;
    trigger = m_MSE_old - m_MSE_new;

    i = i + 1;
end

V_B = V_U / sqrt(v);

t2 = clock;
runtime = etime(t2,t1);
obj.runtime = obj.runtime + runtime;
obj.V_B = V_B;
obj.W_B = W_B;
obj.V_RF = V_RF;
obj.W_RF = W_RF;
obj = get_metric(obj);

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function hypervolume = hypervolume2D(F,ub)
% Copyright (c) 2011, Johannes
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%
% Method for ND objective function values as described in:
%
% 'M. Fleischer. The measure of Pareto Optima Applications to
%  Multi-objective Metaheuristics. EMO 2003, LNCSS 2632
%  519-533, 2003.'
%
% Author: Johannes W. Kruisselbrink
% Last modified: March 17, 2011
%
% Efficient method for 2D objective function values
    F  = -F' + ones(size(F'));
	L  = sortrows(F',1)';
	l  = length(L(1,:)); ub = ub + ones(1,size(L,1));
	hypervolume = 0;
	for i = 1:l
		hypervolume = hypervolume + ((L(1,i) - ub(1)) * (L(2,i) - ub(2)));
        ub(2)       = L(2,i);
    end
end

author: Eric-Bradford

repo: TS-EMO


text:
function [ t, w ] = u_quadrature_rule ( n )

%*****************************************************************************80
%
%% U_QUADRATURE_RULE: quadrature rule for U(n,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real T(N,1), W(N,1), the points and weights of the rule.
%
  aj = zeros ( n, 1 );

  bj = 0.5 * ones ( n, 1 );

  w = zeros ( n, 1 );
  w(1,1) = sqrt ( pi / 2.0 );

  [ t, w ] = imtqlx ( n, aj, bj, w );

  w(1:n,1) = w(1:n,1).^2;

  return
end


text:
function gt = g01_title ( )

%*****************************************************************************80
%
%% G01_TITLE returns the title for grid 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string GT, the grid title.
%
  gt = 'Franke''s 100 node set';

  return
end


text:
function showargs(funcname)
% Show arguments of an argument-aware function with their descriptions
% Usage: showargs funcname

% query arguments
args = arg_report('lean',str2func(funcname));

% print report
fprintf('Arguments of funcname: \n');
for a=args
    if length(a.names) > 1
        name = a.names{2};
    else
        name = a.names{1};
    end
    fprintf(' %s -- %s (%s,%s)\n',name, a.help{1}, a.type, a.shape);
end


text:
% eeg_multieegplot() - Produce an eegplot() of a the average of an epoched dataset
%                  (with optional pre-labelling of specific trials).
% Usage:
%   >> eeg_multieegplot( data,trialrej, elecrej, ...
%                                'key1', value, 'key2', value ... );
% Inputs:
%   data        - input data (channels x points or channels x points x trials).
%   trialrej    - array of 0s and 1s (depicting rejected trials) (size sweeps)
%   elecrej     - array of 0s and 1s (depicting electrodes rejected in
%                 all trials) (size nbelectrodes x sweeps )
%   oldtrialrej - array of 0s and 1s (depicting rejected trials) (size sweeps)
%   oldelecrej  - array of 0s and 1s (depicting electrodes rejected in
%                 all trials) (size nbelectrodes x sweeps )
%
% Note: 1) {'Key', value } Arguments are passed on to eegplot()
%       2) To ignore previous rejections simply set 'oldtrialrej' and
%          'oldelecrej' arguments to empty ([]).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eegplot(), eegplot2event(), eegplot2trial(), 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

% 01-25-02 reformated help & license -ad
% 03-07-02 corrected help -ad

function eeg_multieegplot( data, rej, rejE, oldrej, oldrejE, varargin);

if nargin < 1
   help eeg_multieegplot;
   return;
end;
if nargin < 4
   oldrej = [];
end;
if nargin < 5
   oldrejE = [];
end;
if ~exist('command')
	command = '';
end;
chans = size(data,1);
pnts  = size(data,2);

colnew = [0.8 0.8 1];
colold = [0.8 1 0.8];

if ndims(data) > 2 % --------------- considering epoched EEG data
   rejeegplot = [];
   % concatenate old rejection if not empty
   % --------------------------------------
   if ~isempty( oldrej )
   		oldrej 	 = find( oldrej > 0);
   		oldrejE  = oldrejE(:, oldrej)';
   		rejeegplot = trial2eegplot( oldrej, oldrejE, pnts, colold);
   end;
   % convert for eegplot
   % -------------------
   if ~isempty(	rej )
   		rej    = find( rej > 0);
   		rejE   = rejE(:,rej)';
   		secondrejeegplot = trial2eegplot( rej, rejE, pnts, colnew); % see bottom of code for this function
   		rejeegplot = [ rejeegplot' secondrejeegplot' ]';

   		% remove duplicates
   		% -----------------
		%[tmp I] = unique( rejeegplot(:,1) );
		%rejeegplot = rejeegplot(I,:);
   end;
else % ---------------------------------------- considering continuous EEG
   % for continuous EEG, electrodes (rejE and oldrejE) are not considered yet
   % because there would be a format problem (these rejection are stored in
   % the event array).

   rejeegplot = [];
   if ~isempty(rej) %assuming nrejection x 2 (2 = begin end)
      s = size(rej, 1);
      rejeegplot = [ rej(3:4) colnew*ones(s,1)  zeros(s, chans) ];
   end;

   % pooling previous rejections
   if ~isempty(oldrej)
       tmp = [ oldrej(:, 3:4) colold*ones(s,1) zeros(size(oldrej,1), chans) ];
       rejeegplot = [rejeegplot; tmp];
   end;
end;

if isempty(varargin)
	eegplot(data, 'winlength', 5, 'position', [100 300 800 500], 'winrej', rejeegplot, 'xgrid', 'off');
else
	eegplot(data, 'winlength', 5, 'position', [100 300 800 500], 'winrej', rejeegplot, 'xgrid', 'off', varargin{:} );
end;
return;

% convert eeglab format to eeplot format of rejection window
% ----------------------------------------------------------
function rejeegplot = trial2eegplot( rej, rejE, pnts, color)
   	rejeegplot = zeros(length(rej), size(rejE,2)+5);
   	rejeegplot(:, 6:end) = rejE;
   	rejeegplot(:, 1) = (rej(:)-1)*pnts;
   	rejeegplot(:, 2) = rej(:)*pnts-1;
   	rejeegplot(:, 3:5) = ones(size(rejeegplot,1),1)*color;
return

author: PatternRecognition

repo: OpenBMI


text:
function [dStates, contactForces] = dynamics_doubleStance(States, Actuators, Parameters)
% function [dStates, Actuators] = dynamics_doubleStance(States, Actuators, Parameters)
%
% Computer Generated File -- DO NOT EDIT
%
% This function was created by the function Write_ContinuousDynamics()
% 10-Dec-2013 19:40:47
%
% Dymanics Model: retractable double pendulum biped
% Motion Phase: Double Stance
%
% Matthew Kelly
% Cornell University
%

x = States(1,:); % (m) Foot One horizontal position
y = States(2,:); % (m) Foot One vertical position
th1 = States(3,:); % (rad) Leg One absolute angle
th2 = States(4,:); % (rad) Leg Two absolute angle
L1 = States(5,:); % (m) Leg One length
L2 = States(6,:); % (m) Leg Two length
dx = States(7,:); % (m/s) Foot One horizontal velocity
dy = States(8,:); % (m/s) Foot One vertical velocity
dth1 = States(9,:); % (rad/s) Leg One absolute angular rate
dth2 = States(10,:); % (rad/s) Leg Two absolute angular rate
dL1 = States(11,:); % (m/s) Leg One extension rate
dL2 = States(12,:); % (m/s) Leg Two extensioin rate

F1 = Actuators(1,:); % (N) Compresive axial force in Leg One
F2 = Actuators(2,:); % (N) Compresive axial force in Leg Two
T1 = Actuators(3,:); % (Nm) External torque applied to Leg One
T2 = Actuators(4,:); % (Nm) External torque applied to Leg Two
Thip = Actuators(5,:); % (Nm) Torque acting on Leg Two from Leg One

m1 = Parameters.m1; % (kg) Foot One mass
m2 = Parameters.m2; % (kg) Foot Two mass
M = Parameters.M; % (kg) Hip mass
g = Parameters.g; % (m/s^2) Gravity

% Constraints for this phase:
ddx = 0; %(m) Foot One horizontal position
ddy = 0; %(m) Foot One vertical position

dStates = zeros(size(States));
dStates(1:6,:) = States((1+6):(6+6),:);
dStates(7,:) = zeros(1,size(States,2));
dStates(8,:) = zeros(1,size(States,2));
dStates(9,:) = -(L2.*Thip - L2.*T1 + L1.*T2.*cos(th1 - th2) + L1.*Thip.*cos(th1 - th2) - F2.*L1.*L2.*sin(th1 - th2) + 2.*L1.*L2.*M.*dL1.*dth1 + L1.*L2.*M.*ddy.*cos(th1) + L1.*L2.*M.*g.*cos(th1) - L1.*L2.*M.*ddx.*sin(th1))./(L1.^2.*L2.*M);
dStates(10,:) = (L2.*Thip.*cos(th1).*cos(th2) - L2.*T1.*cos(th1).*cos(th2) - L2.*T1.*sin(th1).*sin(th2) + L2.*Thip.*sin(th1).*sin(th2) - L1.*L2.*M.*ddy.*cos(th2) + L1.*L2.*M.*ddx.*sin(th2) + L1.*T2.*cos(th1 - th2).*cos(th1).*cos(th2) + L1.*Thip.*cos(th1 - th2).*cos(th1).*cos(th2) + L1.*T2.*cos(th1 - th2).*sin(th1).*sin(th2) - L1.*T2.*sin(th1 - th2).*cos(th1).*sin(th2) + L1.*T2.*sin(th1 - th2).*cos(th2).*sin(th1) + L1.*Thip.*cos(th1 - th2).*sin(th1).*sin(th2) - L1.*Thip.*sin(th1 - th2).*cos(th1).*sin(th2) + L1.*Thip.*sin(th1 - th2).*cos(th2).*sin(th1) + F1.*L1.*L2.*cos(th1).*sin(th2) - F1.*L1.*L2.*cos(th2).*sin(th1) - F2.*L1.*L2.*cos(th1 - th2).*cos(th1).*sin(th2) + F2.*L1.*L2.*cos(th1 - th2).*cos(th2).*sin(th1) - F2.*L1.*L2.*sin(th1 - th2).*cos(th1).*cos(th2) - F2.*L1.*L2.*sin(th1 - th2).*sin(th1).*sin(th2) - 2.*L1.*L2.*M.*dL2.*dth2.*cos(th2).^2 - 2.*L1.*L2.*M.*dL2.*dth2.*sin(th2).^2 + L1.*L2.*M.*ddy.*cos(th1).^2.*cos(th2) + L1.*L2.*M.*g.*cos(th1).^2.*cos(th2) - L1.*L2.*M.*ddx.*cos(th1).^2.*sin(th2) + L1.*L2.*M.*ddy.*cos(th2).*sin(th1).^2 + L1.*L2.*M.*g.*cos(th2).*sin(th1).^2 - L1.*L2.*M.*ddx.*sin(th1).^2.*sin(th2))./(L1.*L2.^2.*M.*(cos(th2).^2 + sin(th2).^2));
dStates(11,:) = -(T2.*sin(th1 - th2) + Thip.*sin(th1 - th2) - F1.*L2 + F2.*L2.*cos(th1 - th2) - L1.*L2.*M.*dth1.^2 + L2.*M.*ddx.*cos(th1) + L2.*M.*ddy.*sin(th1) + L2.*M.*g.*sin(th1))./(L2.*M);
dStates(12,:) = (L2.*T1.*cos(th2).*sin(th1) - L2.*T1.*cos(th1).*sin(th2) + L2.*Thip.*cos(th1).*sin(th2) - L2.*Thip.*cos(th2).*sin(th1) - F1.*L1.*L2.*sin(th1).*sin(th2) - L1.*L2.*M.*ddx.*cos(th2) - L1.*L2.*M.*ddy.*sin(th2) + L1.*L2.^2.*M.*dth2.^2.*cos(th2).^2 + L1.*L2.^2.*M.*dth2.^2.*sin(th2).^2 + L1.*T2.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*T2.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*T2.*sin(th1 - th2).*cos(th1).*cos(th2) + L1.*Thip.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*Thip.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*Thip.*sin(th1 - th2).*cos(th1).*cos(th2) + L1.*T2.*sin(th1 - th2).*sin(th1).*sin(th2) + L1.*Thip.*sin(th1 - th2).*sin(th1).*sin(th2) - F1.*L1.*L2.*cos(th1).*cos(th2) + F2.*L1.*L2.*cos(th1 - th2).*sin(th1).*sin(th2) - F2.*L1.*L2.*sin(th1 - th2).*cos(th1).*sin(th2) + F2.*L1.*L2.*sin(th1 - th2).*cos(th2).*sin(th1) + L1.*L2.*M.*ddx.*cos(th1).^2.*cos(th2) + L1.*L2.*M.*ddx.*cos(th2).*sin(th1).^2 + L1.*L2.*M.*ddy.*cos(th1).^2.*sin(th2) + L1.*L2.*M.*g.*cos(th1).^2.*sin(th2) + L1.*L2.*M.*ddy.*sin(th1).^2.*sin(th2) + L1.*L2.*M.*g.*sin(th1).^2.*sin(th2) + F2.*L1.*L2.*cos(th1 - th2).*cos(th1).*cos(th2))./(L1.*L2.*M.*(cos(th2).^2 + sin(th2).^2));

% contactForces(1,:) == H1 == (N) Foot One, horizontal contact force
% contactForces(2,:) == V1 == (N) Foot One, vertical contact force
% contactForces(3,:) == H2 == (N) Foot Two, horizontal contact force
% contactForces(4,:) == V2 == (N) Foot Two, vertical contact force
contactForces = zeros(4,size(States,2));
contactForces(1,:) = (Thip.*sin(th1) - T1.*sin(th1) + L1.*ddx.*m1 + F1.*L1.*cos(th1))./L1;
contactForces(2,:) = (T1.*cos(th1) - Thip.*cos(th1) + L1.*ddy.*m1 + L1.*g.*m1 + F1.*L1.*sin(th1))./L1;
contactForces(3,:) = (L1.*M.*T2.*sin(th2) + L1.*M.*Thip.*sin(th2) + L2.*T1.*m2.*sin(th1) - L2.*Thip.*m2.*sin(th1) - L1.*L2.*M.*ddx.*m2 - L1.*T2.*m2.*cos(th1 - th2).*sin(th1) + L1.*T2.*m2.*sin(th1 - th2).*cos(th1) - L2.*T1.*m2.*cos(th1 - th2).*sin(th2) - L2.*T1.*m2.*sin(th1 - th2).*cos(th2) - L1.*Thip.*m2.*cos(th1 - th2).*sin(th1) + L1.*Thip.*m2.*sin(th1 - th2).*cos(th1) + L2.*Thip.*m2.*cos(th1 - th2).*sin(th2) + L2.*Thip.*m2.*sin(th1 - th2).*cos(th2) - F2.*L1.*L2.*M.*cos(th2) - F1.*L1.*L2.*m2.*cos(th1) + L1.*T2.*m2.*cos(th1 - th2).^2.*sin(th2) + L1.*Thip.*m2.*cos(th1 - th2).^2.*sin(th2) + L1.*T2.*m2.*sin(th1 - th2).^2.*sin(th2) + L1.*Thip.*m2.*sin(th1 - th2).^2.*sin(th2) + L1.*L2.*M.*ddx.*m2.*cos(th1).^2 + L1.*L2.*M.*ddx.*m2.*cos(th2).^2 + L1.*L2.*M.*ddx.*m2.*sin(th1).^2 + L1.*L2.*M.*ddx.*m2.*sin(th2).^2 + F1.*L1.*L2.*m2.*cos(th1 - th2).*cos(th2) + F2.*L1.*L2.*m2.*cos(th1 - th2).*cos(th1) - F1.*L1.*L2.*m2.*sin(th1 - th2).*sin(th2) + F2.*L1.*L2.*m2.*sin(th1 - th2).*sin(th1) - F2.*L1.*L2.*m2.*cos(th1 - th2).^2.*cos(th2) - F2.*L1.*L2.*m2.*sin(th1 - th2).^2.*cos(th2) + L1.*L2.*M.*ddy.*m2.*sin(th1 - th2).*sin(th1).*sin(th2) + L1.*L2.*M.*g.*m2.*sin(th1 - th2).*sin(th1).*sin(th2) - L1.*L2.*M.*ddx.*m2.*cos(th1 - th2).*cos(th1).*cos(th2) + L1.*L2.*M.*ddy.*m2.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*L2.*M.*ddy.*m2.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*L2.*M.*ddy.*m2.*sin(th1 - th2).*cos(th1).*cos(th2) + L1.*L2.*M.*g.*m2.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*L2.*M.*g.*m2.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*L2.*M.*g.*m2.*sin(th1 - th2).*cos(th1).*cos(th2) - L1.*L2.*M.*ddx.*m2.*cos(th1 - th2).*sin(th1).*sin(th2) + L1.*L2.*M.*ddx.*m2.*sin(th1 - th2).*cos(th1).*sin(th2) - L1.*L2.*M.*ddx.*m2.*sin(th1 - th2).*cos(th2).*sin(th1))./(L1.*L2.*M.*(cos(th2).^2 + sin(th2).^2));
contactForces(4,:) = -(L1.*M.*T2.*cos(th2) + L1.*M.*Thip.*cos(th2) + L2.*T1.*m2.*cos(th1) - L2.*Thip.*m2.*cos(th1) + L1.*L2.*M.*ddy.*m2 - L1.*T2.*m2.*cos(th1 - th2).*cos(th1) - L2.*T1.*m2.*cos(th1 - th2).*cos(th2) - L1.*Thip.*m2.*cos(th1 - th2).*cos(th1) + L2.*Thip.*m2.*cos(th1 - th2).*cos(th2) - L1.*T2.*m2.*sin(th1 - th2).*sin(th1) + L2.*T1.*m2.*sin(th1 - th2).*sin(th2) - L1.*Thip.*m2.*sin(th1 - th2).*sin(th1) - L2.*Thip.*m2.*sin(th1 - th2).*sin(th2) + F2.*L1.*L2.*M.*sin(th2) + F1.*L1.*L2.*m2.*sin(th1) + L1.*T2.*m2.*cos(th1 - th2).^2.*cos(th2) + L1.*Thip.*m2.*cos(th1 - th2).^2.*cos(th2) + L1.*T2.*m2.*sin(th1 - th2).^2.*cos(th2) + L1.*Thip.*m2.*sin(th1 - th2).^2.*cos(th2) + F2.*L1.*L2.*m2.*sin(th1 - th2).^2.*sin(th2) - L1.*L2.*M.*ddy.*m2.*cos(th1).^2 - L1.*L2.*M.*ddy.*m2.*cos(th2).^2 - L1.*L2.*M.*g.*m2.*cos(th1).^2 - L1.*L2.*M.*g.*m2.*cos(th2).^2 - L1.*L2.*M.*ddy.*m2.*sin(th1).^2 - L1.*L2.*M.*ddy.*m2.*sin(th2).^2 - L1.*L2.*M.*g.*m2.*sin(th1).^2 - L1.*L2.*M.*g.*m2.*sin(th2).^2 - F1.*L1.*L2.*m2.*cos(th1 - th2).*sin(th2) - F1.*L1.*L2.*m2.*sin(th1 - th2).*cos(th2) - F2.*L1.*L2.*m2.*cos(th1 - th2).*sin(th1) + F2.*L1.*L2.*m2.*sin(th1 - th2).*cos(th1) + F2.*L1.*L2.*m2.*cos(th1 - th2).^2.*sin(th2) + L1.*L2.*M.*ddx.*m2.*sin(th1 - th2).*sin(th1).*sin(th2) + L1.*L2.*M.*ddy.*m2.*cos(th1 - th2).*cos(th1).*cos(th2) + L1.*L2.*M.*g.*m2.*cos(th1 - th2).*cos(th1).*cos(th2) + L1.*L2.*M.*ddx.*m2.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*L2.*M.*ddx.*m2.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*L2.*M.*ddx.*m2.*sin(th1 - th2).*cos(th1).*cos(th2) + L1.*L2.*M.*ddy.*m2.*cos(th1 - th2).*sin(th1).*sin(th2) - L1.*L2.*M.*ddy.*m2.*sin(th1 - th2).*cos(th1).*sin(th2) + L1.*L2.*M.*ddy.*m2.*sin(th1 - th2).*cos(th2).*sin(th1) + L1.*L2.*M.*g.*m2.*cos(th1 - th2).*sin(th1).*sin(th2) - L1.*L2.*M.*g.*m2.*sin(th1 - th2).*cos(th1).*sin(th2) + L1.*L2.*M.*g.*m2.*sin(th1 - th2).*cos(th2).*sin(th1))./(L1.*L2.*M.*(cos(th2).^2 + sin(th2).^2));

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%ex2 simulates the example of Section 3.4.3
%	_______________________________________________________________________
%  /Program:  ex2.m								\
% / Description:  This program runs the simulation of the Rhors example     \
%|  	found in sections 3.4.3 of the text.  The program allows the user	 |
%|  	to choose between the three algorithms defined in the text and 		 |
%|	between the stable or unstable examples.  The program prompts the    |
%|	prompts the user to enter values for the variables that are adjusted |
%|	within the text.  This program calls and requires the fex2.m file.	 |
% \
%  \_______________________________________________________________________/

%Copyright

% Programmed Summer 1996 by Kenyon Thayer, RPI, Troy, NY

global am bm cm ap bp cp Af Bf Cf Df T Tbar G Qp Qf algor a b u

clear yp ym kpe kpf
temp=exist('Af');
if temp==1
	Af_old=Af;
	Bf_old=Bf;
	Cf_old=Cf;
	Df_old=Df;
	G_old=G;
	Qp_old=Qp;
	Qf_old=Qf;
else
	Af_old=0;
	Bf_old=-10;
	Cf_old=1;
	Df_old=0;
	G_old=-5.714;
	Qp_old=57.14;
	Qf_old=20;
end;

am=-3;bm=1;cm=3;

disp('   **This program simulates the text book example 3.4.3**');

algor=input('Which algorithm do you wish to model? 1,2,3: ');
ch2=input('Is this the stable example? Y/N: ','s');

if ch2=='n' | ch2=='N'
	%example 2 (unstable)
	ap=[-7 -92 100;1 0 0;0 1 0];bp=[1;0;0];cp=[0 0 200];
else
	%example 1 (rohrs)
	ap=[-31 -259 -229;1 0 0;0 1 0];bp=[1;0;0];cp=[0 0 458];
end;
% Initial Conditions
% ************************************************************
xp0=[0 0 0];xm0=0;yf0=0;ki0=[0 0 0 0];x0=[xp0,xm0,yf0,ki0]';
t0=0;tfinal=20;tol=1.e-6;
Af=input(sprintf('Enter a value for Af (default is %5.3f):  ',...
	Af_old));
if isempty(Af)
	Af=Af_old;
end
Bf=input(sprintf('Enter a value for Bf (default is %5.3f):  ',...
	Bf_old));
if isempty(Bf)
	Bf=Bf_old;
end
Cf=input(sprintf('Enter a value for Cf (default is %5.3f):  ',...
	Cf_old));
if isempty(Cf)
	Cf=Cf_old;
end
Df=input(sprintf('Enter a value for Df (default is %5.3f):  ',...
	Df_old));
if isempty(Df)
	Df=Df_old;
end
G=input(sprintf('Enter a value for G (default is %5.3f):  ',...
	G_old));
if isempty(G)
	G=G_old;
end
Qp=input(sprintf('Enter a value for Qp (default is %5.3f):  ',...
	Qp_old));
if isempty(Qp)
	Qp=Qp_old;
end
Qf=input(sprintf('Enter a value for Qf (default is %5.3f):  ',...
	Qf_old));
if isempty(Qf)
	Qf=Qf_old;
end

T=1;Tbar=1;

[tout,xout]=ode45('fex2',t0,tfinal,x0,tol,0);

for i=1:length(tout)
	yp(i)=cp*xout(i,[1:3])';
	ym(i)=cm*xout(i,4) ;
	yf=Cf*xout(i,5);
	eyp(i)=ym(i)-yp(i);
	xm=xout(i,4);
	um=0.3;
	if tout(i)>10
		um=-0.3;
	end;
	[ng,ng]=size(G);
	rvec=[eyp(i);-yf;xm;um];
	v1=Qp*eyp(i)-Qf*yf;
	up=inv(eye(ng)-G*rvec'*Tbar*rvec)*(xout(i,[6:9])+v1*rvec'*Tbar)*rvec;
	v=v1+G*up;
	kp=v*rvec'*Tbar;
	kpe(i)=kp(1,1);
	kpf(i)=kp(1,2);
end;

z=max(length(tout),length(ym(1,:)));
if z>length(tout)
	for i=length(tout):length(ym(1,:))
		tout(i,:)=tout((length(tout)),:);
	end;
end;
figure(1)
plot(tout,yp,tout,ym,'--')
axis([0 20 -.5 .5])
text(.41, 1.05, 'Ym', 'color', [1 0 1], 'FontSize', 12, 'Units', 'normal');
text(.51, 1.05, 'Yp', 'color', [1 1 0], 'FontSize', 12, 'Units', 'normal');

author: Sable

repo: mcbench-benchmarks


text:
function [seriesData, dicomdirpath] = load_multi_dcm_dir(dicomdirpath)
% LOAD_MULTI_DCM_DIR  Read metadata from DICOMDIR file.
%
% LOAD_MULTI_DCM_DIR reads metadata from DICOMDIR file and shows list of
% the different DICOM SERIES in Pop-up menus. After selecting a dicom
% series the program creates the dcmSeries structure containing the
% related file names and path.
%
% DICOMDIR is a standard binary file provided with DICOM data, and has
% information such as paths to files, and patient, study, series IDs, etc.
%
% [SERIESDATA, DICOMDIRPATH] = LOAD_MULTI_DCM_DIR()
% SERIESDATA = LOAD_MULTI_DCM_DIR(DICOMDIRPATH)
%
%   DICOMDIRPATH is a string with the path to the directory containing the
%   DICOMDIR file. If not provided, it opens a file browser to select the
%   directory graphically.
%
%   SERIESDATA is a structure containing image information for series
%   selected by the user in the GUI.

% Original Author: Laszlo Balkay, University of Debrecen, PET Center/LB 2008
%      Available from: http://uk.mathworks.com/matlabcentral/fileexchange/7926-dicomdir-reader
%      under BSD 2-Clause License.
%
% Modified by Christopher Kelly <christopher.kelly28@googlemail.com> and
% Benjamin Villard <b.016434@gmail.com>, University of Oxford.
% Copyright © 2014-2015 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% 10 Dec 2014 - Interface changed to give preview of image.
%             - Also allows the selection of multiple series.
%             - Fewer output arguments, now the only output contains the
%               information of the series selected
%
% This file is distributed as a derivative work of a third-party function
% with project Gerardus.
%
% http://code.google.com/p/gerardus/

% check arguments
narginchk(0, 1);
nargoutchk(0, 2);

hm = [];
dcmPatient = [];
ImgPrcTlbx = ver('images');
ImgPrcTlbxVersion = str2double(ImgPrcTlbx.Version(1));
if ImgPrcTlbxVersion < 4
    disp('The required minimum version of ImgProc. Toolbox is 4!');
    return;
end

if ((nargin == 0) || strcmp(dicomdirpath,''))
    dicomdirpath = uigetdir('','Select DICOMDIR folder');
    if dicomdirpath == 0
        return;
    end
end
filename1mat = fullfile(dicomdirpath,'dicom.mat');
dirres1mat = dir(filename1mat);
filename1 = fullfile(dicomdirpath,'dicomdir');
dirres1 = dir(filename1);
filename2 = fullfile(dicomdirpath,'DICOMDIR');
dirres2 = dir(filename2);
if ~isempty(dirres1mat)
    % if the dicomdir mat structure was saved on the first time
    % open it. It is much faster than reading the DICOMDIR dicom file
    load(filename1mat);
    SeriesListNumInfo = dcmdir.SeriesListNumInfo;
    SeriesList = dcmdir.SeriesList;
    dcmPatient = dcmdir.dcmPatient;
    dcmSeries = dcmdir.dcmSeries;
    clear dcmdir;
    disp(' ');
    disp('load dicomdir.mat');
    disp(' ');
else
    if ~isempty(dirres1)
        dcmdir_path = filename1;
    elseif ~isempty(dirres2)
        dcmdir_path = filename2;
    else
        hm = msgbox('No DICOMDIR in the selected folder !','MIA Info','warn' );
        disp('No DICOMDIR in the selected folder !')
        return;
    end

    % Show up message window and set the cursor type to "watch"
    hm = msgbox('Reading the DICOMDIR structure. It takes time. Please wait!','MIA Info' );
    set(hm, 'tag', 'Msgbox_MIA Info');
    SetData=setptr('watch');set(hm,SetData{:});
    hmc = (get(hm,'children'));
    set(hmc(2),'enable','inactive');
    pause(1);

    % open DICOMDIR file & read DirectoryRecordSequence structure
    dcmhdr = dicominfo(dcmdir_path);
    delete(hm);

    % Read the DirectoryRecordSequence structure until end,
    % and read the all included DirectoryRecordType item sequentially
    DirRecords = dcmhdr.DirectoryRecordSequence;
    DirRecordsFieldNames = fieldnames(DirRecords);
    NumOfField = size(DirRecordsFieldNames,1);
    CurrentPatient = 0;
    SeriesListNum = 0;
    waitbar_h = waitbar(0,'Scanning the DICOM files....');
    for i = 1: NumOfField
        CurrentItem = DirRecords.(DirRecordsFieldNames{i});
        CurrentItemType = CurrentItem.DirectoryRecordType;
        if strcmp(CurrentItemType,'PATIENT')
            CurrentPatient = CurrentPatient + 1;
            CurrentStudy = 0;
            if ImgPrcTlbxVersion > 4
                if isfield(CurrentItem,'PatientName')
                    dcmPatient(CurrentPatient,1).PatientName = ...
                        CurrentItem.PatientName.FamilyName;
                else
                    dcmPatient(CurrentPatient,1).PatientName = 'No value available';
                end
            elseif ImgPrcTlbxVersion == 4
                if isfield(CurrentItem,'PatientsName')
                    dcmPatient(CurrentPatient,1).PatientName = ...
                        CurrentItem.PatientsName.FamilyName;
                else
                    dcmPatient(CurrentPatient,1).PatientName = 'No value available';
                end
            end
        elseif strcmp(CurrentItemType,'STUDY')
            CurrentStudy = CurrentStudy + 1;
            CurrentSeries = 0;
            if isfield(CurrentItem,'StudyDescription')
                dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).StudyDescription = ...
                    CurrentItem.StudyDescription;
            else
                dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).StudyDescription = '';
            end
        elseif strcmp(CurrentItemType,'SERIES')
            if CurrentSeries > 0 % create summary about the previous read series for popupmenu
                SeriesListNum = SeriesListNum + 1;
                Pname = dcmPatient(CurrentPatient,1).PatientName;
                Modality = dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1).Modality;
                SeriesDescription = dcmPatient(CurrentPatient,1).Study(CurrentStudy,1). ...
                    Series(CurrentSeries,1).SeriesDescription;
                StudyDescription = dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).StudyDescription;
                NumOfImages = num2str(CurrentImage);
            end
            CurrentSeries = CurrentSeries + 1;
            CurrentImage = 0;
            dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1).Modality = ...
                CurrentItem.Modality;
            if isfield(CurrentItem,'SeriesDescription')
                dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1). ...
                    SeriesDescription = CurrentItem.SeriesDescription;
            else
                dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1). ...
                    SeriesDescription = '';
            end
            %%%
        elseif strcmp(CurrentItemType,'IMAGE')
            CurrentImage = CurrentImage + 1;
            [pathstr, dcmfname, ext] = fileparts(CurrentItem.ReferencedFileID);
            if CurrentImage == 1
                dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1). ...
                    ImagePath = pathstr;
            end
            if ispc
                dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1). ...
                    ImageNames{CurrentImage,1} = CurrentItem.ReferencedFileID;
            else
                dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1). ...
                    ImageNames{CurrentImage,1} = strrep(CurrentItem.ReferencedFileID,'\','/');
            end
            %ImageNames{CurrentImage,1} = dcmfname;
            dcmPatient(CurrentPatient,1).Study(CurrentStudy,1).Series(CurrentSeries,1). ...
                ReferencedSOPInstanceUIDInFile{CurrentImage,1} = CurrentItem.ReferencedSOPInstanceUIDInFile;
            %%%
        end
        waitbar(i/NumOfField,waitbar_h);
        drawnow;
    end
    close(waitbar_h);

    try
        dcmdir.dcmPatient = dcmPatient;
        save(filename1mat,dcmdir);
    catch
        % Unable to write file 'filename1mat' permission denied.
    end

end

global dcmVal;
dcmVal = 0;

if size(dcmPatient,1) > 1
    selectPatientFigure = figure('Units','normalized','menubar','none',...
        'NumberTitle','off','name','selectPatient',...
        'position',[2/8, 3/8, 4/8, 2/8]);
    uicontrol('Parent',selectPatientFigure,'Units','normalized','Style','text',...
        'String','Select which patient you are analysing?','Position',[0.1 0.7 0.8 0.2]);
    String  = sprintf('%s|', dcmPatient.('PatientName'));
    String(end) = [];
    lbh = uicontrol('Parent',selectPatientFigure,'Units','normalized','Style','popupmenu',...
        'String',String,'Position',[0.1 0.4 0.8 0.2],'tag','patientList');
    okCloseCallback = 'global dcmVal; dcmVal = get(findobj(''tag'',''patientList''),''value''); delete(findobj(''name'',''selectPatient''));';
    uicontrol('Parent',selectPatientFigure,'Units','normalized','Style','pushbutton','String','OK','Position',[0.1 0.1 0.8 0.2],'Callback',okCloseCallback);
    uiwait(selectPatientFigure);
else
    dcmVal = 1;
end

dcmPatient = dcmPatient(dcmVal,:);

clearvars -except dcmPatient dicomdirpath

set(0,'DefaultFigureWindowStyle','normal');

screensize = get(0,'ScreenSize');
dcmdirlistfh = figure('Units','normalized','menubar','none','NumberTitle','off',...
    'name','DICOMDIR info','position',[0.1, 0.1, 0.8, 0.8]);

lbh = uicontrol('Parent',dcmdirlistfh,'Units','normalized','Style','listbox',...
    'Position',[0.05, 0.09, 0.45, 0.88 ],'tag','dcmdirlist_popupmenu');

[seriesDescriptions,order] = sort({dcmPatient.Study.Series.SeriesDescription});
dcmPatient.Study.Series = dcmPatient.Study.Series(order);

set(lbh,'string',seriesDescriptions,'max',1.1);

global dcmdirlistVal;
dcmdirlistVal = 0;

OKCallback = 'global a; a = 1; global dcmdirlistVal; dcmdirlistVal = get(findobj(''tag'',''dcmdirlist_popupmenu''),''value''); delete(findobj(''name'',''DICOMDIR info'')); return;';
CancelCallback = 'global a; a = 1; delete(findobj(''name'',''DICOMDIR info''));';
OK_h = uicontrol('Parent',dcmdirlistfh,'Units','normalized','Style','pushbutton',...
    'String','OK','Position',[0.15, 0.03, 0.2, 0.03],'Callback',OKCallback);
Cancel_h = uicontrol('Parent',dcmdirlistfh,'Units','normalized','Style','pushbutton',...
    'String','Cancel','Position',[0.65, 0.03, 0.2, 0.03],'Callback',CancelCallback);

global a
a = 0;
axes1 = axes('Parent',dcmdirlistfh,'Units','normalized','Position',[0.55, 0.55, 0.4, 0.4]); axis square;
while a == 0
    dcmdirlistVal = get(findobj('tag','dcmdirlist_popupmenu'),'value');
    pause(0.1);
    if length(dcmdirlistVal) >= 1.1
        continue
    end
    try
        imageCurrent = double(uint16(dicomread(fullfile(dicomdirpath,dcmPatient.Study.Series(dcmdirlistVal).ImageNames{1}))));
        imagesc(imageCurrent); axis square; colormap gray; axis off;
        pause(0.1);

        strDetails = sprintf('File Location : %s \n\n Series Description : %s \n\n Number of Frames : %s \n\n ',...
            num2str(fullfile(dicomdirpath,dcmPatient.Study.Series(dcmdirlistVal).ImageNames{1})),...
            num2str(dcmPatient.Study.Series(dcmdirlistVal).SeriesDescription),...
            num2str(numel(dcmPatient.Study.Series(dcmdirlistVal).ImageNames)));

        uicontrol('Parent',dcmdirlistfh,'Units','normalized','Style','text',...
            'Position',[0.55, 0.09, 0.4, 0.4],'String',strDetails,'FontWeight','bold');
    end
end

close all;

if dcmdirlistVal == 0;
    disp('No DICOM Series was selected!');
    return;
end

seriesData = dcmPatient.Study.Series(dcmdirlistVal);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function M = ImaGIN_mexhat(fs, f)
% Generate Mexican hat wavelets
% FORMAT M = spm_eeg_morlet(Rtf, ST, f, ff)
%
% fs  - sample frequency [ms]
% f   - vector of frequencies [Hz]
%
% M   - cell vector, where each element contains the filter for each
%       frequency in f
%__________________________________________________________________________
%
% spm_eeg_morlet generates morlet wavelets for specified frequencies f with
% a specified ratio Rtf, see [1], for sample time ST (ms). One obtains the
% wavelet coefficients by convolution of a data vector with the kernels in
% M. See spm_eeg_tf how one obtains instantaneous power and phase estimates
% from the wavelet coefficients.
%
% [1] C. Tallon-Baudry, O. Bertrand, F. Peronnet and J. Pernier, 1998.
% Induced \gamma-Band Activity during the Delay of a Visual Short-term
% memory Task in Humans. The Journal of Neuroscience (18): 4244-4254.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_eeg_morlet.m 5900 2014-02-27 21:54:51Z karl $

M      = {};
for f0 = f

    nu=f0/fs;
    M{end+1} = mexhat(nu);

end

function  h = mexhat(nu) ;
%MEXHAT	Mexican hat wavelet in time domain.
%	H=MEXHAT(NU) returns the mexican hat wavelet with central
%	frequency NU (NU is a normalized frequency in Hz).
%
%	NU : any real between 0 and 0.5		(default : 0.05).
%	H  : time vector containing the mexhat samples.
%		length(H)=2*ceil(1.5/NU)+1.
%
%	Example :
%	 plot(mexhat);
%
%	See also KLAUDER.

%	P. Goncalves, October 95
%	Copyright (c) 1995 Rice University
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if nargin==0,
 nu=0.05;
elseif nargin>1,
 error('Too many input parameters');
end

if nu>0.5 | nu<0,
 disp('Warning : NU should be between 0 and 0.5');
end

N = 1.5 ;
alpha = pi^2*nu^2 ;
n = ceil(N/nu) ;
t = -n:n ;
h = nu*sqrt(pi)/2*exp(-alpha*t.^2).*(1-2*alpha*t.^2) ;

author: brainstorm-tools


text:
function [xi,w]=secondOrderPrismCubPoints()
%%SECONDORDERPRISMCUBPOINTS Generate second-order cubature points for
%  integration over a standard prism in 3D. This is a triangle that has
%  been extruded upward. The vertices are (1,-1,-1), (-1,-1,-1), (-1,1,-1),
%  (1,-1,1), (-1,-1,1), and (-1,1,1).
%
%INPUTS: None
%
%OUTPUTS: xi This is a 3XnumCubPoints set of points for the standard prism.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the standard prism (4).
%
%This function implements the points given in [1] (5 points).
%
%EXAMPLE:
%We compare a 2nd-order moment computed using these cubature points
%to one computed using monomialIntPrism. The results are the same within
%typical finite precision limits.
% [xi,w]=secondOrderPrismCubPoints();
% alpha=[2;0;0];
% theMoment=findMomentFromSamp(alpha,xi,w);
% intVal=monomialIntPrism(alpha);
% RelErr=(theMoment-intVal)/intVal
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[-0.33333333333333333333333333333333333333,  -0.33333333333333333333333333333333333333,  -0.99999999999999985119303811076522004183,   0.66666666666666686507594918564641756458;
-0.33333333333333333333333333333333333333,  -0.33333333333333333333333333333333333333,   0.99999999999999985119303811076522004183,   0.66666666666666686507594918564641756458;
-0.74158162379719638007464124598498266439,   0.48316324759439276014928249196996532878,                                          0,   0.88888888888888875661603387623572162361;
 0.48316324759439276014928249196996532878,  -0.74158162379719638007464124598498266439,                                          0,   0.88888888888888875661603387623572162361;
-0.74158162379719638007464124598498266439,  -0.74158162379719638007464124598498266439,                                          0,   0.88888888888888875661603387623572162361];

w=M(:,4);
xi=M(:,1:3)';

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function logPr = calcLogPrGaussian( X, mu, invSigma)
% Each row of X is an independent draw from Normal( mu, invSigma )

D = size( invSigma,1);
LOG_2PI = 1.837877066409345;

cholInvSigma = chol( invSigma );
logDetInvSigma = 2*sum( log( diag( cholInvSigma ) ) );

XdiffMu = bsxfun( @minus, X, mu );
U = XdiffMu'*cholInvSigma';
logPr = -0.5*D*LOG_2PI +  0.5*logDetInvSigma  - 0.5*sum(U.^2, 2);

author: michaelchughes

repo: NPBayesHMM


text:
function video2pic( videoFile, outputPath, picFormat )
%VIDEO2PIC converts viedo frames into pictures.
%   This function uses mmreader to read videos. Thus video format supported
%   by system can be converted to pictures. Output file name follows the
%   format: VideoName_FrameNumber.PicSuffix.
%
%   VIDEO2PIC input parameters interactively.
%   VIDEO2PIC( VIDEOFILE, OUTPUTPATH, PICFORMAT ) runs as a
%   function.
%
%   videoFile: full path of video to convert.
%   outputPath: pictures will be saved under outputPath/videoname/.
%   picFormat: check 'help imwrite'
%
%	Shiquan Wang @ CASIA
%	sqwang@nlpr.ia.ac.cn
%   http://www.cbsr.ia.ac.cn/users/sqwang/
%	revised on 2009/10/19

%% Initialization
if nargin < 1
    videoFile = input('Input full video path:', 's');
    outputPath = input('Input video output path[Default is under the video path]:', 's');
    if isempty(outputPath)
        outputPath = videoFile(1, 1:find(videoFile=='\', 1, 'last')-1);
    end
    picFormat = input('Input output pic format[JPG by default]:');
    if isempty(picFormat)
        picFormat = 'JPG';
    end
end

%% Main part
videoObject = mmreader(videoFile);
numFrames = get(videoObject, 'numberOfFrames');
fileName = get(videoObject, 'Name');
% bitsPerPix = get(videoObject, 'BitsPerPixel');
% H = get(videoObject, 'Height');
% W = get(videoObject, 'Width');
% bytesPerFrame = bitsPerPix*H*W/8;
outputPath = fullfile(outputPath, fileName(1, 1:find(fileName=='.', 1, 'last')-1));
numOrder = max(4, size(int2str(numFrames), 2));
mkdir(outputPath);

indexFrame = [1 1];
frameLimit = 100; % This can be set a larger value with enough memory.
while indexFrame(1, 1) <= numFrames
%     memorynow = memory;
%     memoryMax = memorynow.MaxPossibleArrayBytes;
%     frameLimit = floor(memoryMax/bytesPerFrame*0.8);
    indexFrame(1, 2) = min(numFrames, indexFrame(1, 1) + frameLimit);
    frameAll = read(videoObject, indexFrame);

    for i = indexFrame(1, 1):1:indexFrame(1, 2)    %save frames to pic
        imgFrame = frameAll(:,:,:,i - indexFrame(1, 1) + 1);
        saveFormat = strcat('%s\\%s_%0', int2str(numOrder), 'd.%s');
        picName = sprintf(saveFormat, outputPath, fileName, i, picFormat);
        imwrite(imgFrame, picName);
        disp(picName); % This output can be turned off.
    end
    indexFrame(1, 1) = indexFrame(1, 2) + 1;
end

end


text:
function CP = runICPRejection(p, g, CP)

msg('S', {g.procICP{:} 'REJECTION'}, 'LogLevel', 'basic');

for i = 1:size(p.PairList,1)

    CP{i} = CP{i}.rejection('Attribute', ...
                            'AttributeName', 'dAlpha', ...
                            'AttributeMinMax', [p.MaxDeltaAngle*10/9 400]); % conversion from degree to gradian!

    % Rejection based on dp
    med = median(CP{i}.dp);
    sig_mad = 1.4826*mad(CP{i}.dp,1);

    CP{i} = CP{i}.rejection('Attribute', ...
                            'AttributeName', 'dp', ...
                            'AttributeMinMax', [-Inf med-p.MaxSigmaMad*sig_mad]);

    CP{i} = CP{i}.rejection('Attribute', ...
                            'AttributeName', 'dp', ...
                            'AttributeMinMax', [med+p.MaxSigmaMad*sig_mad Inf]);

    % Rejection based on ds
    CP{i} = CP{i}.rejection('Attribute', ...
                            'AttributeName', 'ds', ...
                            'AttributeMinMax', [p.MaxDistance Inf]);

    % Rejection based on roughness
    CP{i} = CP{i}.rejection('Attribute', ...
                            'AttributeName', 'roughness', ...
                            'AttributeMinMax', [p.MaxRoughness Inf]);

end

msg('E', {g.procICP{:} 'REJECTION'}, 'LogLevel', 'basic');

end

author: pglira