text:
function readme()
%
% I. INTRODUCTION
%     The MATLAB(R) program contained within this directory counts the total
% number of loops (cycles) in a network (graph) that consists of nodes and edges.
%     This file describes the setup/installation proceedures for the code. A
% description of the algorithm and how to use the code is contained in the
% 'DETAILS.m' file.
%
%
% II. CONTENTS
%     In the Loops/ directory, you should find the following files:
%         DETAILS.m
%         README.m
%         loops_gui.m
%         run_loops.m
%         nets/**
%         **several sample network files
%
%
% III. SETUP/INSTALLATION
%     1. You must have MATLAB software installed on your computer
%     2. Copy/move the 'Loops' folder to the MATLAB 'work' directory
%     3. Open MATLAB and add the 'Loops' directory to the Path
%         a. Go to: File -> Set Path...
%         b. Click on 'Add with subfolders...'
%         c. Select the 'Loops' directory and click 'Ok'
%         d. Click 'Save'
%         e. Click 'Close'
%     4. To run, open the 'loops_gui.m' file and press F5 or use the
%         command line: >> loops_gui
%     5. Note: this code was written using MATLAB 7R14 through R2006B. It has not
%         been tested on previous versions
%
%
% IV. HELP/REPORT BUGS
%     If you experience difficulties using this program, first make sure that
% the steps in Section III have been completed. Next, make sure your network
% satisfies all the requirements given in the 'DETAILS.m' file.
%
% Please direct questions/comments to:
% Joe Kirk
% jdkirk630@gmail.com
%
%
% V. REVISION NOTES
%     11/2005 Update:
%         1. New file 'loop_gui.m' - GUI file that replaces 'run_loops.m' and
%             displays all of the tools for the user as they are available
%         2. Added ability to save loops in .MAT format
%     10/2005 Update:
%         1. New file 'reduceNet.m' - function which allows networks to be reduced
%             (removes nodes that have only one edge, until no more remain in the net)
%         2. New file 'getStartingNode.m' - function which calculates a (nearly)
%             optimal starting node to make the ILCA more efficient (results in
%             fewer steps to complete the algorithm)
%         3. Removed file 'printNetStats.m' - the basic functionality of this
%             file was separated into two separate files ('calcNumEdges.m'
%             and 'plotHLoops.m')
%         4. New file 'calcNumEdges.m' - function which calculates the number of
%             edges in a network
%         5. New file 'plotHLoops.m' - function which plots the distribution of
%             loops of length 'h'
%         6. Modified file 'generateRandomNet.m' to limit the creation of sparse
%             networks with more nodes than 40 (which are costly to generate with
%             this function)
%     2/2007 Update:
%         1. All of the subfunction files have been deprecated, and their
%             functionality has been combined and added to the end of the
%             'loops_gui.m' and 'run_loops.m' files. The two files 'loops_gui.m'
%             and 'run_loops.m' are now independent, stand-alone m-files
%         2. Added ability to save the network as an edgelist file
%         3. Improved the layout and function of the Loops GUI
%
clc
help readme

author: Sable

repo: mcbench-benchmarks


text:
function [kg] = amu2kg(amu)
% Convert mass from atomic mass units to kilograms.
% Chad Greene 2012
kg = amu*1.6605402e-27;


text:
function im=lpcar2im(ar,np)
%LPCAR2IM Convert AR coefs to impulse response IM=(AR,NP)

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

[nf,p1]=size(ar);
if nargin<2 np=p1-1; end
im=zeros(nf,np+1);
x=[1 zeros(1,np)];
for k=1:nf
  im(k,:)=filter(1,ar(k,:),x);
end

author: decouples


text:
function bnet = mk_hhmm2(varargin)
% MK_HHMM2 Make a 2 level Hierarchical HMM
% bnet = mk_hhmm2(...)
%
% 2-layer hierarchical HMM  (node numbers in parens)
%
%   Q1(1) ---------> Q1(5)
% /  | \            / |
% |  |  v          /  |
% |  |  F2(3) --- /   |
% |  |  ^         \   |
% |  | /           \  |
% |  v              \ v
% |  Q2(2)--------> Q2 (6)
% |  |
% \  |
%  v v
%   O(4)
%
%
% Optional arguments [default]
%
% discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
% obsCPT       - CPT(o,q1,q2) params for O ['rnd']
% mu           - mu(:,q1,q2) params for O [ [] ]
% Sigma        - Sigma(:,q1,q2) params for O [ [] ]
%
% F2toQ1       - 1 if Q2 is an hhmm_CPD, 0 if F2 -> Q2 arc is absent, so level 2 never resets [1]
% Q1args        - arguments to be passed to the constructors for Q1(t=2) [ {} ]
% Q2args        - arguments to be passed to the constructors for Q2(t=2) [ {} ]
%
% F2 only turns on (wp 0.5) when Q2 enters its final state.
% Q1 (slice 1) is clamped to be uniform.
% Q2 (slice 1) is clamped to always start in state 1.

[os nmodels nstates] = size(mu);

ss = 4;
Q1 = 1; Q2 = 2; F2 = 3; obs = 4;
Qnodes = [Q1 Q2];
names = {'Q1', 'Q2', 'F2', 'obs'};
intra = zeros(ss);
intra(Q1, [Q2 F2 obs]) = 1;
intra(Q2, [F2 obs]) = 1;

inter = zeros(ss);
inter(Q1,Q1) = 1;
inter(F2,Q1) = 1;
if F2toQ2
  inter(F2,Q2)=1;
end
inter(Q2,Q2) = 1;

ns = zeros(1,ss);

ns(Q1) = nmodels;
ns(Q2) = nstates;
ns(F2) = 2;
ns(obs) = os;

dnodes = [Q1 Q2 F2];
if discrete_obs
  dnodes = [dnodes obs];
end
onodes = [obs];

bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
eclass = bnet.equiv_class;

% SLICE 1

% We clamp untied nodes in the first slice, since their params can't be estimated
% from just one sequence

% uniform prior on initial model
CPT = normalise(ones(1,nmodels));
bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', CPT, 'adjustable', 0);

% each model always starts in state 1
CPT = zeros(ns(Q1), ns(Q2));
CPT(:, 1) = 1.0;
bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', CPT, 'adjustable', 0);

% Termination probability
CPT = zeros(ns(Q1), ns(Q2), 2);
if 1
  % Each model can only terminate in its final state.
  % 0 params will remain 0 during EM, thus enforcing this constraint.
  CPT(:, :, 1) = 1.0; % all states turn F off ...
  p = 0.5;
  CPT(:, ns(Q2), 2) = p; % except the last one
  CPT(:, ns(Q2), 1) = 1-p;
end
bnet.CPD{eclass(F2,1)}  = tabular_CPD(bnet, F2, 'CPT', CPT);

if discrete_obs
  bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, obs_args{:});
else
  bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, obs_args{:});
end

% SLICE 2

bnet.CPD{eclass(Q1,2)} = hhmm_CPD(bnet, Q1+ss, Qnodes, 1, D, 'args', Q1args);

if F2toQ2
  bnet.CPD{eclass(Q2,2)} = hhmmQD_CPD(bnet, Q2+ss, Qnodes, 2, D, Q2args{:});
else
  bnet.CPD{eclass(Q2,2)} = tabular_CPD(bnet, Q2+ss, Q2args{:});
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% SP_HCURL_ERROR: Evaluate the error in H(curl) norm.
%
%   [errhcurl, errl2, errcurl] = sp_hcurl_error (space, msh, u, uex, curluex);
%
% INPUT:
%
%    space:   object defining the space of discrete functions (see sp_vector)
%    msh:     object defining the domain partition and the quadrature rule (see msh_cartesian)
%    u:       vector of dof weights
%    uex:     function handle to evaluate the exact solution
%    curluex: function handle to evaluate the curl of the exact solution
%
% OUTPUT:
%
%     errhcurl: error in H(curl) norm
%     errl2:    error in L^2 norm
%     errcurl:  error of the curl in L^2 norm
%
% Copyright (C) 2010 Carlo de Falco, Rafael Vazquez
% Copyright (C) 2011, 2015 Rafael Vazquez
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [errhcurl, errl2, errcurl] = sp_hcurl_error (space, msh, u, uex, curluex)

  if (numel(u) ~= space.ndof)
    error ('Wrong size of the vector of degrees of freedom')
  end

  errl2 = 0;
  errcurl = 0;

  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    sp_col  = sp_evaluate_col (space, msh_col, 'value', true, 'curl', true);

    [~, err_l2, err_curl] = sp_hcurl_error (sp_col, msh_col, u, uex, curluex);

    errcurl = errcurl + err_curl.^2;
    errl2 = errl2 + err_l2.^2;
  end

  errhcurl = sqrt (errl2 + errcurl);
  errl2 = sqrt (errl2);
  errcurl = sqrt (errcurl);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function value = p12_exact ( dim_num )

%*****************************************************************************80
%
%% P12_EXACT returns the exact integral for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  value = 1.0;
  for i = 1 : dim_num
    value = value * sin ( i );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% loss_sqp - squared loss function
%
% Copyright(c) 2009 Ryota Tomioka
% This software is distributed under the MIT license. See license.txt
function [floss, gloss]=loss_sqp(zz, bb)

gloss = zz-bb;
floss = 0.5*sum(gloss.^2);

repo: REST


text:
function PA=partialMatrixInverse(A,n)
%%PARTIALMATRIXINVERSE Given an mXm possibly singular matrix A of a rank
%   >=n, where the first n columns are independent, obtain a matrix that is
%   akin to the inverse related to the first n columns. A normal matrix
%   inverse is such that A*inv(A)=eye(m,m). Here, consider A=Q*R, where Q
%   is an orthogonal matrix and R is a lower-triangular matrix (not the
%   standard QR decomposition). The inverse relation is Q*R*inv(A)=eye(m,m)
%   or R*inv(A)=Q'*eye(m,m). This function solves for the nXn upper-left
%   submatrix of inv(A).  This function is useful for taking a singular
%   Fisher informationmatrix, where the first n rows are observable, but
%   the rest are not, and obtaining a Cramer-Rao lower bound for the
%   observable components. Typically, this might be obtaining a CRLB on
%   position components before velocity is fully observable.
%
%INPUTS: A An mXm matrix.
%        n A number n<=m that is <= the rank of A.
%
%OUTPUTS: PA The partial inverse matrix as described above.
%
%This function solves the problem by calling linEqSolveFirstnRows with A
%and vectors [1;0;0;0;...], [0;1;0;0;...], etc. The solution coincides with
%a submatrix of the pseudoinverse.
%
%EXAMPLE:
%This example of the CRLB of a position estimate of a target given a single
%range, direction cosine and range-rate measurement is considered. We show
%that when this function is used, the position accuracy implied by the CRLB
%is not affected by the accuracy of the range rate component (which one
%would assume). However, using the ad-hoc method of taking the inverse of a
%subset of the Fisher information matrix, one gets a different (wrong)
%answer the more accurate the range rate component is. That is due to not
%properly accounting for matrix cross terms.
% %Sensor location near Maui.
% llhRxTx=[deg2rad([20.888645;-156.022474]);0];
% %Convert to Cartesian.
% xRxTx=zeros(6,1);
% %Fill in the position. The reciever will be stationary (0 velocities).
% %Take the first one to be the transmitter.
% xRxTx(1:3,:)=ellips2Cart(llhRxTx);
%
% %Sensor orientation.
% elAboveLevel=deg2rad(15);%The radar point slightly up.
% %Rotation matrices of the radar.
% MRxTx=findRFTransParam(llhRxTx(:,1),pi/2,elAboveLevel);%Facing East
%
% %Construct a target state vector.
% llhTar=[deg2rad([21.14146;-155.058839]);13e3];
% tarLocCart=ellips2Cart(llhTar);
% tarHeading=deg2rad(-140);%Radians East of North.
% angUpFromLevel=0;%Level flight.
% tarSpeed=300;%m/s
% vECEFTar=geogHeading2uVec(llhTar,tarHeading,angUpFromLevel)*tarSpeed;
% xTar=[tarLocCart;vECEFTar];
%
% %Measuremnet and noise parameters of the radar.
% useHalfRange=false;
% sigmaR=10;%Range standard deviation.
% sigmaU=0.01;
% sigmaV=0.01;
% sigmaRR=1;%Range rate standard deviation, meters/ second.
% sigmaRRMassive=1e5;%Effectively uninformative standard deviation.
% %Lower-triangular square root measurement covariance matrix (No cross
% %terms).
% SR=diag([sigmaR;sigmaU;sigmaV;sigmaRR]);
% RInv=inv(SR*SR');
% %The same thing but with the range rate standard deviation so massive it is
% %essentially uninformative.
% SRMassive=diag([sigmaR;sigmaU;sigmaV;sigmaRRMassive]);
% RInvMassive=inv(SRMassive*SRMassive');
%
% %Gradients for the Fisher information matrices.
% H=calcRuvRRJacob(xTar,useHalfRange,xRxTx(:,1),xRxTx(:,1),MRxTx(:,:,1));
% FIM=observedFisherInfo([],RInv,[],H);
% %The position RMSE from the CRLB found correctly using this function.
% RMSE=sqrt(sum(diag(partialMatrixInverse(FIM,3))))
% %An ad-hoc (wrong) approach that one might be tempted to take:
% RMSEWrong=sqrt(sum(diag(inv(FIM(1:3,1:3)))))
% %One can see that the position RMSE is higher when done the correct way.
% %However, a single range rate on its own doesn't inform on the position
% %accuracy. If we do the same thing cranking up the standard deviation of
% %the range rate so it is uninformative, we get:
% FIM=observedFisherInfo([],RInvMassive,[],H);
% RMSEUninfRR=sqrt(sum(diag(partialMatrixInverse(FIM,3))))
% RMSEUninfRRWrong=sqrt(sum(diag(inv(FIM(1:3,1:3)))))
% %In other words, the correct solution doesn't change, but in this case the
% %wrong solution has become correct, because the cross terms that messed it
% %up are essentially gone.
%
%February 2022 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

m=size(A,1);
PA=zeros(n,n);

if(n>m)
    error('n cannot be > the length of A.')
end

e=zeros(m,1);
for k=1:n
    e(k)=1;
    PA(:,k)=linEqSolveFirstnRows(A,e,n);
    e(k)=0;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function filtered = Kuwahara(original,winsize)
%Kuwahara   filters an image using the Kuwahara filter
%   filtered = Kuwahara(original,windowSize) filters the original image with a
%                                            given windowSize and yields the result in filtered
%
% This function is optimised using vectorialisation, convolution and
% the fact that, for every subregion
%     variance = (mean of squares) - (square of mean).
% A nested-for loop approach is still used in the final part as it is more
% readable, a commented-out, fully vectorialised version is provided as
% well.
%
% This function is about 2.3 times faster than KuwaharaFast at
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=13474&objectType=file
% with a 5x5 window, and even faster at higher window sizes (about 4 times on a 13x13 window)
%
% Inputs:
% original      -->    image to be filtered
% windowSize    -->    size of the Kuwahara filter window: legal values are
%                      5, 9, 13, ... = (4*k+1)
%
% Example
% filtered = Kuwahara(original,5);
%
% Filter description:
% The Kuwahara filter works on a window divided into 4 overlapping
% subwindows (for a 5x5 pixels example, see below). In each subwindow, the mean and
% variance are computed. The output value (located at the center of the
% window) is set to the mean of the subwindow with the smallest variance.
%
%    ( a  a  ab   b  b)
%    ( a  a  ab   b  b)
%    (ac ac abcd bd bd)
%    ( c  c  cd   d  d)
%    ( c  c  cd   d  d)
%
% References:
% http://www.ph.tn.tudelft.nl/DIPlib/docs/FIP.pdf
% % http://www.incx.nec.co.jp/imap-vision/library/wouter/kuwahara.html
%
% Copyright Luca Balbi, 2007

%% Incorrect input handling
error(nargchk(2, 2, nargin, 'struct'));

% non-double data will be cast
if ~isa(original, 'double')
    original = double(original);
end % if

% wrong-sized kernel is an error
if mod(winsize,4)~=1
    error([mfilename ':IncorrectWindowSize'],'Incorrect window size: %d',winsize)
end % if

%% Build the subwindows
tmpAvgKerRow = [ones(1,(winsize-1)/2+1) zeros(1,(winsize-1)/2)];
tmpPadder = zeros(1,winsize);
tmpavgker = repmat(tmpAvgKerRow,(winsize-1)/2+1,1);
tmpavgker = [tmpavgker; repmat(tmpPadder,(winsize-1)/2,1)];
tmpavgker = tmpavgker/sum(sum(tmpavgker));

% tmpavgker is a 'north-west' subwindow (marked as 'a' above)
% we build a vector of convolution kernels for computing average and
% variance
avgker(:,:,1) = tmpavgker;                  % North-west (a)
avgker(:,:,2) = fliplr(tmpavgker);      % North-east (b)
avgker(:,:,4) = flipud(tmpavgker);      % South-east (c)
avgker(:,:,3) = fliplr(avgker(:,:,4));      % South-west (d)

% this is the (pixel-by-pixel) square of the original image
squaredImg = original.^2;

% preallocationg these arrays makes it about 15% faster
avgs = zeros([size(original) 4]);
stddevs = zeros([size(original) 4]);

%% Calculation of averages and variances on subwindows
for k=1:4
    avgs(:,:,k) = conv2(original,avgker(:,:,k),'same');      % mean on subwindow
    stddevs(:,:,k) = conv2(squaredImg,avgker(:,:,k),'same'); % mean of squares on subwindow
    stddevs(:,:,k) = stddevs(:,:,k)-(avgs(:,:,k)).^2;        % variance on subwindow
end % for

%% Choice of the index with minimum variance
[minima,indices] = min(stddevs,[],3); %#ok<ASGLU>

%% Building of the filtered image (with nested for loops)
filtered = zeros(size(original));
for k=1:size(original,1)
    for n=1:size(original,2)
        filtered(k,n) = avgs(k,n,indices(k,n));
    end % for
end % for

%% Commented out, completely vectorialised alternative
% [y,x] = meshgrid(1:size(original,2),1:size(original,1));
% lookupIndices = x+size(original,1)*(y-1)+numel(original)*(indices-1);
% filtered = avgs(lookupIndices);

end % function


text:
function [dist, pos] = distancePointLine(point, line)
%DISTANCEPOINTLINE Minimum distance between a point and a line
%
%   D = distancePointLine(POINT, LINE)
%   Return the euclidean distance between line LINE and point POINT.
%
%   LINE has the form: [x0 y0 dx dy], and POINT is [x y].
%
%   If LINE is N-by-4 array, result is N-by-1 array computes for each line.
%
%   If POINT is N-by-2, then result is computed for each point.
%
%   If both POINT and LINE are array, result is computed for each couple of
%   point and line, and is returned in a NP-by-NL array, where NP is the
%   number of points, and NL is the number of lines.
%
%
%   See also:
%   lines2d, points2d, distancePoints, distancePointEdge
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2005-06-24
% Copyright 2016 INRA - BIA-BIBS.

%   HISTORY:
%   2012-10-24 rewrite using bsxfun

% direction vector of each line (row vectors)
vx = line(:, 3)';
vy = line(:, 4)';

% squared norm of direction vectors, with a check of validity
delta = (vx .* vx + vy .* vy);
invalidEdges = delta < eps;
delta(invalidEdges) = 1;

% difference of coordinates between point and line origins
% (NP-by-NE arrays)
dx  = bsxfun(@minus, point(:, 1), line(:, 1)');
dy  = bsxfun(@minus, point(:, 2), line(:, 2)');

% compute position of points projected on the line, by using normalised dot
% product
% (result is a NP-by-NL array)
pos = bsxfun(@rdivide, bsxfun(@times, dx, vx) + bsxfun(@times, dy, vy), delta);

% ensure degenerated lines are correclty processed (consider the line
% origin as closest point)
pos(:, invalidEdges) = 0;

% compute distance between point and its projection on the line
dist = hypot(bsxfun(@times, pos, vx) - dx, bsxfun(@times, pos, vy) - dy);

% if size(line, 1)==1 && size(point, 1)>1
%     line = repmat(line, [size(point, 1) 1]);
% end
%
% if size(point, 1)==1 && size(line, 1)>1
%     point = repmat(point, [size(line, 1) 1]);
% end
%
% dx = line(:, 3);
% dy = line(:, 4);
%
% % compute position of points projected on line
% tp = ((point(:, 2) - line(:, 2)).*dy + (point(:, 1) - line(:, 1)).*dx) ./ (dx.*dx+dy.*dy);
% p0 = line(:, 1:2) + [tp tp].*[dx dy];
%
%
% % compute distances between points and their projections
% dx = point - p0;
% dist  = sqrt(sum(dx.*dx, 2));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ normX, normY, normZ, meanShape, Transform ] = ProcrustesAnalysis3D( x, y, z, tangentSpace, meanShape )
%PROCRUSTESANALYSIS3D Summary of this function goes here
%   Detailed explanation goes here

meanProvided = false;

if(nargin > 4)
    meanProvided = true;
end

% Translate all elements to origin
normX = zeros(size(x));
normY = zeros(size(y));
normZ = zeros(size(z));

for i = 1:size(x,1)

    offsetX = mean(x(i,:));
    offsetY = mean(y(i,:));
    offsetZ = mean(z(i,:));

    Transform.offsetX(i) = offsetX;
    Transform.offsetY(i) = offsetY;
    Transform.offsetZ(i) = offsetZ;

    normX(i,:) = x(i,:) - offsetX;
    normY(i,:) = y(i,:) - offsetY;
    normZ(i,:) = z(i,:) - offsetZ;

end

% Rotate elements untill all of them have the same orientation

% the initial estimate of rotation would be the first element
% if change is less than 1% stop (shouldn't take more than 2 steps)
change = 0.1;

if(~meanProvided)
    meanShape = [ mean(normX); mean(normY); mean(normZ) ]';
end
% scale all the shapes to mean shape

% Get the Frobenius norm, to scale the shapes to mean size (still want to
% retain mm)
meanScale = norm(meanShape, 'fro');

for i = 1:size(x,1)

    scale = norm([normX(i,:) normY(i,:) normZ(i,:)], 'fro')/meanScale;

    normX(i,:) = normX(i,:)/scale;
    normY(i,:) = normY(i,:)/scale;
    normZ(i,:) = normZ(i,:)/scale;

end

Transform.RotationX = zeros(size(x,1),1);
Transform.RotationY = zeros(size(x,1),1);
Transform.RotationZ = zeros(size(x,1),1);

for i = 1:30

    % align all of the shapes to the mean shape

    % remember all orientations to get the mean one (in euler angle form, pitch, yaw roll)
    orientationsX = zeros(size(normX,1),1);
    orientationsY = zeros(size(normX,1),1);
    orientationsZ = zeros(size(normX,1),1);

    for j = 1:size(x,1)

        currentShape = [normX(j,:); normY(j,:); normZ(j,:)]';
        % we want to align the current shape to the mean one
        [ R, T ] = AlignShapesKabsch(currentShape, meanShape);

        eulers = Rot2Euler(R);

        orientationsX(j) = eulers(1);
        orientationsY(j) = eulers(2);
        orientationsZ(j) = eulers(3);

        Transform.RotationX(j) = eulers(1);
        Transform.RotationY(j) = eulers(2);
        Transform.RotationZ(j) = eulers(3);

        currentShape = R * currentShape';

        normX(j,:) = currentShape(1,:);
        normY(j,:) = currentShape(2,:);
        normZ(j,:) = currentShape(3,:);

    end

    % recalculate the mean shape
%     if(~meanProvided)
        oldMean = meanShape;
        meanShape = [mean(normX); mean(normY); mean(normZ)]';
        meanScale = norm(meanShape, 'fro');
%     end

    for j = 1:size(x,1)

        scale = norm([normX(j,:) normY(j,:) normZ(j,:)], 'fro')/meanScale;

        normX(j,:) = normX(j,:)/scale;
        normY(j,:) = normY(j,:)/scale;
        normZ(j,:) = normZ(j,:)/scale;

    end

    if(i==1 && ~meanProvided)

        % rotate the mean shape to mean rotation
        meanOrientationX = mean(orientationsX);
        meanOrientationY = mean(orientationsY);
        meanOrientationZ = mean(orientationsZ);

        R = Euler2Rot([meanOrientationX, meanOrientationY, meanOrientationZ]);
        meanShape = (R * meanShape')';
    end

    % find frobenious norm
    diff = norm(oldMean - meanShape, 'fro');

    if(diff/norm(oldMean,'fro') < change)
        break;
    end

end

% transform to tangent space to preserve linearities

% get the scaling factors for each shape
if(tangentSpace)
    [ normX, normY, normZ] = TangentSpaceTransform(normX, normY, normZ, meanShape);
end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function f=ref_iedgtii_1(c,g,a,M)
%REF_EDGTII_1   Reference Inverse Even DGT type II by DGT
%   Usage  c=ref_edgt(f,g,a,M);
%
%   The input window must be odd-centered of length 2L.
%
%   a must be divisable by 2.

L=size(g,1)/2;
W=size(c,2);

N=L/a;

clong=zeros(M,2*N,W);

cr=reshape(c,M,N,W);
% Copy the first half unchanged
clong(:,1:N,:)=cr;

% Copy the non modulated coefficients.
clong(1,N+1:2*N,:)=cr(1,N:-1:1,:);

% Copy the modulated coefficients.
clong(2:M,N+1:2*N,:)=-cr(M:-1:2,N:-1:1,:);

clong=reshape(clong,2*M*N,W);

fdouble=ref_igdgt(clong,g,a,M,.5,0,floor(a/2));

f=fdouble(1:L,:);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
num = size(AllFeature1,2);
F1 = AllFeature1';
% F1 = sqrt(F1);
F1 = bsxfun(@rdivide, F1, sqrt(sum(F1.^2,2)));
F1 = [F1 F11];
% F1 = bsxfun(@minus,F1,PCAmap.mean);
% F1 = F1 * PCAmap.M;
F2 = AllFeature2';
% F2 = sqrt(F2);
F2 = bsxfun(@rdivide, F2, sqrt(sum(F2.^2,2)));
F2 = [F2 F21];
% F2 = bsxfun(@minus,F2,PCAmap.mean);
% F2 = F2 * PCAmap.M;
% F1 = AllFeature1';
% F2 = AllFeature2';
thresh2 = zeros(num,1);
for i = 1:num
%     thresh2(i) = F1(i,:) * mapping.A * F1(i,:)' + F2(i,:) * mapping.A * F2(i,:)' - 2 * F1(i,:) * mapping.G * F2(i,:)';
    thresh2(i) = pdist2(F1(i,:),F2(i,:));
%     thresh2(i) = F1(i,:) * F2(i,:)';
end;
figure;
hist(thresh2(1:3000),500);
figure;
hist(thresh2(3001:end),500);

accuracies = zeros(10,1);
for i=1:10
    test_idx = [(i-1) * 300 + 1 : i*300, (i-1) * 300 + 3001 : i*300 + 3000];
    train_idx = 1:6000;
    train_idx(test_idx) = [];
    bestc=256;
    same_label = ones(6000,1);
    same_label(3001:6000) = 0;
    % predicted_label = predict(double(lfw_label),sparse(thresh2),model);
    cmd = [' -t 0 -h 0'];
    model = svmtrain(same_label(train_idx),thresh2(train_idx),cmd);
    % model = svmtrain(double(sim_label),thresh,cmd);
    [class, accuracy, deci] = svmpredict(same_label(test_idx),thresh2(test_idx),model);
    accuracies(i) = accuracy(1);
end;
mean(accuracies)
cmd = [' -t 0 -h 0'];
model = svmtrain(same_label,thresh2,cmd);
[class, accuracy, deci] = svmpredict(same_label,thresh2,model);
% mean(thresh2(same_label==1)) / 4 + mean(max(0,1 - thresh2(same_label==0))) / 4
% sum((thresh2<0.22) == same_label) / 6000

author: happynear

repo: FaceVerification


text:
close all
clear

%------------------------------------------------------------------------------------------
nb = 1000;%Total number of initial states
nmpcd=zeros(1,nb);%Store computing time for one time-step for NMPC-DCBF
impcd=zeros(1,nb);%Store computing time for one time-step for iMPC-DCBF
nmpcinf=zeros(1,1);%Store infeasible rate for one time-step for NMPC-DCBF
impcinf=zeros(1,1);%Store infeasible rate for one time-step for iMPC-DCBF
i=24;%i is number of horozon
ts = 1;%Total time steps
load(gamma1_0p6gamma2_0p6\InitialStateData');%Load the generated 1000 initial states
[a, b, c, d]=compare(i,ts,nb,InitialMat);%Involves NMPC-DCBF and iMPC-DCBF
nmpcd(1,:)=a;
impcd(1,:)=b;
nmpcinf(1,1)=c;
impcinf(1,1)=d;
nmpcplot1=nmpcd(1,:);
nmpcplot1(nmpcplot1==-1)=[];%Eliminate infeasible one time-step trajectories
impcplot1=impcd(1,:);
impcplot1(impcplot1==-1)=[];%Eliminate infeasible one time-step trajectories

save('timecom24','nmpcplot1','impcplot1');
save('feasibility24','nmpcinf','impcinf');

distnmpc = fitdist(nmpcplot1','Normal');
distimpc = fitdist(impcplot1','Normal');
mu1=distnmpc.mu;%Get mean of sample of computing time for NMPC-DCBF
mu2=distimpc.mu;%Get mean of sample of computing time for iMPC-DCBF
sigma1=distnmpc.sigma;%Get variance of sample of computing time for NMPC-DCBF
sigma2=distimpc.sigma;%Get variance of sample of computing time for iMPC-DCBF

%Initialize atmosphere parameters
function [tnmpc, timpc, ratiotnmpc, ratiotimpc]=compare(N11,ttsim,samplen,InitialMat)
tnmpc=[];
timpc=[];
for i=1:samplen
N1=N11;
tsim=ttsim;
xini1=InitialMat(1,i);
yini1=InitialMat(2,i);
thetaini1=InitialMat(3,i);
vini1=InitialMat(4,i);
x01=[xini1;yini1;thetaini1;vini1];
x02=[xini1 yini1 thetaini1 vini1];
t1 = nmpcdcbf(x01, N1, tsim);
t2 = impcdcbf(x02, N1, tsim);
tindex=[N11 i];
disp(tindex);
tnmpc=[tnmpc t1];%Computing time for NMPC-DCBF
timpc=[timpc t2];%Computing time for iMPC-DCBF
end
nnmpc1 = length(tnmpc);
nimpc1 = length(timpc);
tnmpcs = tnmpc;
timpcs = timpc;
tnmpcs(tnmpcs==-1)=[];
timpcs(timpcs==-1)=[];
nnmpc2 = length(tnmpcs);
nimpc2 = length(timpcs);
ratiotnmpc = (nnmpc1-nnmpc2)/nnmpc1;%Infeasible rate for NMPC-DCBF
ratiotimpc = (nimpc1-nimpc2)/nimpc1;%Infeasible rate for iMPC-DCBF
end

function t1 = nmpcdcbf(x00, N1, ttsim)
%% General Flags
run_nmpc_dcbf_one = true;
run_nmpc_dcbf_multiple = true;

%% Setup and Parameters
x0 = x00;
dt = 0.1;
time_total = ttsim *dt;
P = zeros(4,4);%Weight matrix P
P(1,1) = 10; P(2,2) = 10; P(3,3) = 10;P(4,4) = 10;
Q = zeros(4,4);%Weight matrix Q
Q(1,1) = 10; Q(2,2) = 10; Q(3,3) = 10;Q(4,4) = 10;
R = zeros(4,4);%Weight matrix R
R(1,1) = 1; R(2,2) = 1;R(3,3) = 1000;R(4,4) = 1000;
%Variables range as below
xmin = [-10; -10; -10;-10];
xmax = [10; 10; 10;10];
umin = [-7; -5;-inf;-inf];
umax = [7; 5;inf;inf];

%% Discrete-time unicycle model
system.dt = dt;
system.A = [1 0 0 0;
    0 1 0 0;
    0 0 1 0;
    0 0 0 1];
system.B = [x0(4,1)*cos(x0(3,1))*dt;
    x0(4,1)*sin(x0(3,1))*dt;
    0;
    0];
system.C =[0 0 0 0;
    0 0 0 0;
    1*dt 0 0 0;
    0 1*dt 0 0];
system.xl = xmin;
system.xu = xmax;
system.ul = umin;
system.uu = umax;

%% NMPC-DCBF parameters
params_nmpc_dcbf.Q = Q;
params_nmpc_dcbf.R = R;
params_nmpc_dcbf.P = P;
params_nmpc_dcbf.gamma1 = 0.6;
params_nmpc_dcbf.gamma2 = 0.6;

%% Obstacle
obs.pos1 = [0; 0];
obs.r1 = 1;

%% Simulate NMPC-DCBF with other N
params_nmpc_dcbf.N = N1;
if run_nmpc_dcbf_one
    fprintf('Run NMPC-DCBF-24\n');
    controller_nmpc_dcbf_multiple24 = NMPCDCBF2(x0, system, params_nmpc_dcbf);%Highest order mcbf=2
    controller_nmpc_dcbf_multiple24.obs = obs;
    controller_nmpc_dcbf_multiple24.sim(time_total);
end
t1=controller_nmpc_dcbf_multiple24.tt;
end

function t2=impcdcbf(x00, N1, ttsim)
t2=0;
xo = 0;
yo = 0;
r = 1;
N = N1;%Number of horizon
K1 = 1000;%Maximum iteration times, jmax
K = 1000;%Maximum iteration times, jmax
dt = 0.1;
gamma1 = 6;
gamma2 = 6;
nsim = ttsim;

tic;

% Constraints

umin = [-7; -5; -inf; -inf;];
umax = [7; 5; Inf; Inf];
xmin = [-10; -10; -10; -10];
xmax = [ 10;  10;  10;  10];

% Objective function
Q = diag([10 10 10 10]);
QN = Q;
R = 1 * eye(4);
R(3,3) = 1000;
R(4,4) = 1000;
% Initial and reference states
x0 = x00;
xr = [3; 0.01; 0; 0];
ur = [0; 0; 1; 1];

% Dynamic system initialization
BB = [0 0 0 0;0 0 0 0;dt 0 0 0;0 dt 0 0];

% Convex MPC frame
[nx, nu] = size(BB);
% Initialize states set (x00,x01,,,x0N）
x_0 = [];
x0 = x0';
u_0 = zeros(nu, N);
u0 = zeros(nu, 1);% we may need to change this for a better warm-up start
impc = zeros(nx, nsim + 1);
impc(:, 1) = x0;
x0new = x0;
for i=1 : (N+1)
    x_0 = [x_0 x0new];
    x0new = [dt*x0new(4)*cos(x0new(3))+x0new(1);dt*x0new(4)*sin(x0new(3))+x0new(2);dt*u0(1)+x0new(3);dt*u0(2)+x0new(4)];
end

abc = tangentline(x_0, r);
Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);%First order CBF constraints
lcbf = -abc(3, :);
lcbf(1) = [];
lcbf = lcbf';
ucbf = inf * ones(N, 1);
A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);%Second order CBF constraints
lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
ucbf2 = inf * ones(N-1, 1);

% Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
% - quadratic objective
P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
% - linear objective
q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
% - linear dynamics
Ax = getax(x_0, dt, N, nx);
Bu = kron([sparse(1, N); speye(N)], BB);
Cx = getcx(x_0, u_0, dt, N, nx);
Aeq = [Ax, Bu];
leq = [-x0; -Cx];
ueq = leq;
% - input and state constraints
Aineq = speye((N+1)*nx + N*nu);
lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

% - OSQP constraints
A = [Aeq; Aineq; Acbf; A2cbf];
l = [leq; lineq; lcbf; lcbf2];
u = [ueq; uineq; ucbf; ucbf2];

% Create an OSQP object
prob = osqp;
% Setup workspace
prob.setup(P, q, A, l, u, 'warm_start', true);
% Solve
res = prob.solve();

% Check solver status
if ~strcmp(res.info.status, 'solved')
%    error('OSQP did not solve the problem!')
        t2=-1;
        return
end
ctrlx = res.x(1:(N+1)*nx);
ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);

for j = 1 : (K1-1)%Iterations for the first time-step
    storagex = ctrlx;%store the state variables and inputs in order to compare them with the next iteration
    storageu = ctrlu;
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
end % Move for the first step
ctrl = ctrlu(1:nu);
x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
ctrlu = rewrite(ctrlu, nu, N);
t2=t2+toc;
return

x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
u_0 = transu(ctrlu, nu, N);
impc(:, 2) = x0;
storagex = ctrlx;
storageu = ctrlu;

for i = 1 : (nsim-1)%Iterations for the left time steps

   for j = 1 : K
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    x0 = ctrlx(1:nx);
    testx0 = [testx0 x0];
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
    storagex = ctrlx;
    storageu = ctrlu;
   end
   ctrl = ctrlu(1:nu);
   x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
   ctrlu = rewrite(ctrlu, nu, N);
   x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
   u_0 = transu(ctrlu, nu, N);
   impc(:, i+2) = x0;
end
end

% Linerize the CBF constraints (get a, b, c for lines)
function abc = tangentline(xy, r)% x and y from initialize states set, abc are coeeficients for linear equation a*x+b*y+c=0
[xx, ~] = size(xy);%xx=2,yy=N+1
xy(xx,:) = []; % this part should be changed for other case
xy((xx-1),:) = []; % this part should be changed other case
[xx, yy] = size(xy);%xx=2,yy=N+1
xyjiao = zeros(xx, yy);%intersection points
for i = 1 : xx
    for j = 1 : yy
        xyjiao(i, j) = r * xy(i, j) * (1 / (xy(:, j)' * xy(:, j)))^(0.5);%calculate coordinates of intersection points
    end
end
cc = -r^2 * ones(1, yy);
abc = [xyjiao; cc];
end

% Get CBF constraints matrix
function Acbf = cbfmat(abc, nx, nu, dt, gamma, x0)
[~, yy] = size(abc);
Acbfx = zeros((yy-1), yy*nx);
Acbfu = zeros((yy-1), (yy-1)*nu);
for i = 1 : (yy-1)
    Acbfx(i, (i*nx)+1) = abc(1, (i+1));
    Acbfx(i, (i*nx)+2) = abc(2, (i+1));
end
for i = 1 : (yy-1)
    Acbfu(i, ((i-1)*nu+3)) = - (1 - dt * gamma)^(i) * (abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
Acbf = [Acbfx Acbfu];
end

% Transfer vector x into matrix
function res = trans(x0, vector, nxx, nyy)%nxx=nx,nyy=N+1
 res = zeros(nxx, nyy);
 res(:,1) = x0;
 for i = 1 : (nyy -1)
     res(:,i+1)= vector(((i)*nxx+1):(i+1)*nxx);
 end
end

% Transfer vector u into matrix
function resu = transu(vector, nxx, nyy)%nxx=nu,nyy=N
 resu = zeros(nxx, nyy);
 for i = 1 : (nyy)
     resu(:,i)= vector(((i-1)*nxx+1):(i)*nxx);
 end
end

% Rewrite u vector
function reu = rewrite(vector, nu, N)
append = vector((N-1)*nu+1:N*nu);
vector(1:nu) = [];
reu = [vector;append];
end

% Get new x_0
function x_0 = newinit(x0, ctrlu, nx, nu, N, dt)
 x_0 = zeros(nx, N+1);
 x_0(:, 1) = x0;
 for i=1 : N
    u0 = ctrlu((i-1)*nu+1:i*nu);
    x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*u0(1)+x0(3);dt*u0(2)+x0(4)];
    x_0(:, i + 1) = x0;
 end
end

% Get AA matrix
function AA = getaa(x0, dt)
 AA = [1 0 -x0(4)*sin(x0(3))*dt cos(x0(3))*dt;0 1 x0(4)*cos(x0(3))*dt sin(x0(3))*dt;0 0 1 0;0 0 0 1];
end
% Get CC matrix
function CC = getcc(x0, x1, u0, dt)
 CC = [x0(4)*sin(x0(3))*x0(3)*dt-x0(4)*cos(x0(3))*dt+x1(1)-x0(1);-x0(4)*cos(x0(3))*x0(3)*dt-x0(4)*sin(x0(3))*dt+x1(2)-x0(2);-u0(1)*dt+x1(3)-x0(3);-u0(2)*dt+x1(4)-x0(4)];
end
% Get Ax matrix
function Ax = getax(x_0, dt, N, nx)
x0 = x_0(:,1);
AA = getaa(x0, dt);
Ax = kron(speye(N+1), -speye(nx)) + kron(sparse(diag(ones(N, 1), -1)), AA);
for i = 1 : (N-1)
   x0 = x_0(:,i+1);
   AA = getaa(x0, dt);
   Ax(nx*(i+1)+1:nx*(i+1)+nx,nx*i+1:nx*i+nx) = AA;
end
end
% Get Cx matrix
function Cx = getcx(x_0, u_0, dt, N, nx)
Cx = zeros(N*nx, 1);
for i = 1 : N
    u0 = u_0(:,i);
    x0 = x_0(:,i);
    x1 = x_0(:,i+1);
    CC = getcc(x0, x1, u0, dt);
    Cx((i-1)*nx+1:(i-1)*nx+nx) = CC;
end
end
% Get A2cbf
function A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0)
[~, yy] = size(abc);
Acbfx2 = zeros((yy-2), yy*nx);
Acbfx22 = zeros((yy-2), yy*nx);
Acbfu2 = zeros((yy-2), (yy-1)*nu);
for i = 1 : (yy-2)
    Acbfx2(i, (i*nx)+1) = (gamma1-1/dt)*abc(1, (i+1));
    Acbfx2(i, (i*nx)+2) = (gamma1-1/dt)*abc(2, (i+1));
    Acbfx2(i, ((i+1)*nx)+1) = (1/dt)*abc(1, (i+2));
    Acbfx2(i, ((i+1)*nx)+2) = (1/dt)*abc(2, (i+2));
end
for i = 1 : (yy-2)
    Acbfx22(i, (1*nx)+1) = -(1-dt*gamma2)^(i)/dt*abc(1, (1+1));
    Acbfx22(i, (1*nx)+2) = -(1-dt*gamma2)^(i)/dt*abc(2, (1+1));
end
Acbfx2 = Acbfx2 + Acbfx22;
for i = 1 : (yy-2)
    Acbfu2(i, ((i-1)*nu+4)) = - (1 - dt * gamma2)^(i)*(gamma1-1/dt)*(abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
A2cbf = [Acbfx2 Acbfu2];
end
% Get lcbf2
function lcbf2 = getlcbf2(abc, dt, gamma1, gamma2)
[~, yy] = size(abc);
lcbf2 = zeros((yy-2),1);
for i = 1 : (yy-2)
    lcbf2(i, 1) = -abc(3, (i+2))/dt-(gamma1-1/dt)*abc(3, (i+1))+(1 - dt * gamma2)^(i)*abc(3, (1+1))/dt;
end
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function loadInputFileCallback(grfBlockHandle)
% This function runs an input file and loads the contents into the block
% masks.

    %% Run the input file selected in the global reference frame
    % Find input filename
    mask = Simulink.Mask.get(grfBlockHandle);
    InputFile = mask.getParameter('InputFile');
    inputFilePath = InputFile.Value;

    %% Set class parameters in block mask
    % Get all simulink blocks
    blocks = find_system(bdroot,'Type','Block');

    % Loop through blocks, writing each to the input file
    for i=1:length(blocks)
        % Variable names and values of a block
        names = get_param(blocks{i},'MaskNames');
        blockHandle = getSimulinkBlockHandle(blocks{i});

        % Check if the block is from the WEC-Sim library
        if any(contains(names,{'simu','waves'})) % Global reference frame
            type = 0;
        elseif any(contains(names,{'body'})) % flexible or rigid body
            type = 1;
        elseif any(contains(names,{'pto'})) % pto
            type = 2;
        elseif any(contains(names,{'constraint'})) % constraint
            type = 3;
        elseif any(contains(names,{'mooring'})) && any(contains(names,{'stiffness'})) % mooring matrix
            type = 4;
        elseif any(contains(names,{'mooring'})) && any(contains(names,{'moorDynLines'})) % moorDyn
            type = 5;
        elseif any(contains(names,{'cable'}))
            type = 6;
        end

        writeBlocksFromInput(blockHandle,type,inputFilePath);
        % write blocks again to account for read only params that are now open
        if type==1 || type==0
            writeBlocksFromInput(blockHandle,type,inputFilePath);
        end
    end
end

author: WEC-Sim

repo: WEC-Sim

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


text:
function [p, dp] = pne_pfcn_pseudo_arc_len(x, xp, step, zp)
%PNE_PFCN_PSEUDO_ARC_LEN  Pseudo arc length parameterization function for PNE
%   [P, DP] = PNE_PFCN_PSEUDO_ARC_LEN(X, XP, STEP, ZP)
%
%   Inputs:
%       X    : solution vector x (last element is parameter lambda)
%       XP   : previous solution vector
%       STEP : continuation parameter step size
%       ZP   : normalized tangent vector at XP
%
%   This function defines a pseudo arc length parameterization for a PNE,
%   where the current point on the solution curve is constrained to lie in
%   the hyperplane running through the predicted solution orthogonal to the
%   tangent line from the previous corrected solution.
%
%   Outputs:
%       P : value of parameterization function
%       DP : Jacobian of paramerization function (transpose of gradient)
%
%   See also PNE_PFCN_NATURAL, PNE_PFCN_ARC_LEN.

%   MP-Opt-Model
%   Copyright (c) 2013-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Shrirang Abhyankar, Argonne National Laboratory
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

dp = zp';                   %% derivative
p = dp * (x - xp) - step;   %% function

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function imdb = getCliveDatabase_train(Dir)

imdb.imgpath = cell(1,1162);
imgpath = fullfile(Dir,'Data','AllImages_release.mat');
img = load(imgpath);
img = img.AllImages_release;

mospath = fullfile(Dir,'Data','AllMOS_release.mat');
mos = load(mospath);
mos = mos.AllMOS_release;
imdb.mos = mos(8:end);

for i = 8:1169
    file_name = img{i,1};
    imdb.imgpath{i-7} = fullfile('Images',file_name);
end

sel = randperm(1162);
train_sel = sel(1:round(0.8*1162));
test_sel = sel(round(0.8*1162)+1:end);

train_path = imdb.imgpath(train_sel);
test_path = imdb.imgpath(test_sel);

train_mos = imdb.mos(train_sel);
test_mos = imdb.mos(test_sel);

imdb.images.id = 1:1162 ;
imdb.images.set = [ones(1,size(train_path,2)),2*ones(1,size(test_path,2))];
imdb.images.label = [train_mos,test_mos];

imdb.classes.description = {'LIVE_CHAN'};
imdb.images.name = [train_path,test_path] ;
imdb.imageDir = Dir ;
% fclose(fileID);

repo: DBCNN


text:
function mind = dijkstra_distance ( nv, ohd )

%*****************************************************************************80
%
%% DIJKSTRA_DISTANCE uses Dijkstra's minimum distance algorithm.
%
%  Discussion:
%
%    We essentially build a tree.  We start with only node 0 connected
%    to the tree, and this is indicated by setting CONNECTED(0) = 1.
%
%    We initialize MIND(I) to the one step distance from node 0 to node I.
%
%    Now we search among the unconnected nodes for the node MV whose minimum
%    distance is smallest, and connect it to the tree.  For each remaining
%    unconnected node I, we check to see whether the distance from 0 to MV
%    to I is less than that recorded in MIND(I), and if so, we can reduce
%    the distance.
%
%    After NV-1 steps, we have connected all the nodes to 0, and computed
%    the correct minimum distances.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2010
%
%  Author:
%
%    Original C version by Norm Matloff, CS Dept, UC Davis.
%    This MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, integer OHD(NV,NV), the distance of the direct link between
%    nodes I and J.
%
%    Output, integer MIND(NV), the minimum distance from node 1 to each node.
%

%
%  Start out with only node 1 connected to the tree.
%
  connected(1) = 1;
  connected(2:nv) = 0;
%
%  Initialize the minimum distance to the one-step distance.
%
  mind(1:nv) = ohd(1,1:nv);
%
%  Attach one more node on each iteration.
%
  for step = 2 : nv
%
%  Find the nearest unconnected node.
%
    [ md, mv ] = find_nearest ( nv, mind, connected );

    if ( mv == - 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'DIJKSTRA_DISTANCE - Warning!\n' );
      fprintf ( 1, '  Search terminated early.\n' );
      fprintf ( 1, '  Graph might not be connected.\n' );
      return
    end
%
%  Mark this node as connected.
%
    connected(mv) = 1;
%
%  Having determined the minimum distance to node MV, see if
%  that reduces the minimum distance to other nodes.
%
    mind = update_mind ( nv, mv, connected, ohd, mind );

  end

  return
end


text:
function [xPred, PPred] = KalmanFilterX_PredictState(x,P,F,Q,u,B,Qu)
% KALMANFILTERX_PREDICTSTATE Perform the discrete-time KF state prediction
% step, under the assumption of additive process noise.
%
% Parameters
% ----------
% x: column vector
%   The (xDim x 1) state estimate at the previous time-step.
% P: matrix
%   The (xDim x xDim) state covariance matrix at the previous
%   time-step.
% F: matrix
%   An (xDim x xDim) state transition matrix.
% Q: matrix
%   The (xDim x xDim) process noise covariance matrix.
% u: column vector, optional
%   An optional (xDim x 1) control input.
%   If omitted, no control input is used.
% B: matrix, optional
%   An optional (xDim x xDim) control gain matrix.
%   If omitted, B is assumed to be 1.
% O: matrix, optional
%   An optional (xDim x xDim) control noise covariance
%   matrix. If omitted, Q is assumed to be 0.
%
% Returns
% -------
% xPred: column vector
%   The (xDim x 1) predicted state estimate.
% PPred: matrix
%   The (xDim x xDim) predicted state covariance matrix.
%
%October 2017 Lyudmil Vladimirov, University of Liverpool.

    switch(nargin)
        case(4)
            u  = 0;
            B  = 0;
            Qu = 0;
        case(5)
            B  = 1;
            Qu = 0;
        case(6)
            Qu = 0;
    end

    % Compute predicted state mean and covariance
    xPred = F*x + B*u;
    PPred =F*P*F' + Q + B*Qu*B';
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [elem2dof,elem2edge,edge,bdDof,freeDof] = dofP3(elem)
%% DOFP3 dof structure for P3 element.
%
%  [elem2dof,edge,bdDof] = DOFP3(elem) constructs the dof structure
%  for the quadratic element based on a triangle. elem2dof(t,i) is the
%  global index of the i-th dof of the t-th element.
%
%  The global indices of the dof is organized  according to the order of
%  nodes, edges and elements, namely, first give index number to the dofs
%  on nodes, then the dofs on edges, last the dofs on elements.
%
%  See also dofP2, dof3P3.
%
%  Created by Jie Zhou.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N = max(max(elem)); NT = size(elem,1);

%% Data structure
totalEdge = uint32(sort([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])],2));
matlabversion = version;
if str2double(matlabversion(end-5:end-2)) > 2012
    [edge, i2, j] = unique(totalEdge,'rows','legacy');
else
    [edge, i2, j] = unique(totalEdge,'rows');
end
NE = size(edge,1);
elem2edge = reshape(j,NT,3);

%% Nodal dof
elem2dof = uint32(zeros(NT,10));
elem2dof(:,1:3) = elem;

%% Two dof on each edge
% edge 1
idx0 = (elem(:,3) > elem(:,2));
elem2dof(idx0,4) = N + 2*(elem2edge(idx0,1))-1;
elem2dof(idx0,5) = N + 2*(elem2edge(idx0,1));
elem2dof(~idx0,4)= N + 2*(elem2edge(~idx0,1));
elem2dof(~idx0,5)= N + 2*(elem2edge(~idx0,1))-1;
% edge 2
idx0 = (elem(:,3) > elem(:,1));
elem2dof(idx0,6) = N + 2*(elem2edge(idx0,2));
elem2dof(idx0,7) = N + 2*(elem2edge(idx0,2))-1;
elem2dof(~idx0,6) = N + 2*(elem2edge(~idx0,2))-1;
elem2dof(~idx0,7) = N + 2*(elem2edge(~idx0,2));
% edge 3
idx0 = (elem(:,2) > elem(:,1));
elem2dof(idx0,8) = N + 2*(elem2edge(idx0,3))-1;
elem2dof(idx0,9) = N + 2*(elem2edge(idx0,3));
elem2dof(~idx0,8) = N + 2*(elem2edge(~idx0,3));
elem2dof(~idx0,9) = N + 2*(elem2edge(~idx0,3))-1;

%% Element dof
elem2dof(:,10) = (N+2*NE+1:N+2*NE+NT)';

%% Boundary dof
i1(j(3*NT:-1:1)) = 3*NT:-1:1;
i1 = i1';
bdEdgeIdx = (i1 == i2);
isBdDof = false(N+2*NE+NT,1);
isBdDof(edge(bdEdgeIdx,:)) = true;   % boundary node
idx = find(bdEdgeIdx);
isBdDof(N+2*idx) = true;      % two dof on boundary edges
isBdDof(N+2*idx-1) = true;
bdDof = find(isBdDof);
 freeDof = find(~isBdDof);

repo: Aero-matlab


text:
function lambda = orth_symm_eigenvalues ( n )

%*****************************************************************************80
%
%% ORTH_SYMM_EIGENVALUES returns eigenvalues of the ORTH_SYMM matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  m = floor ( ( n + 1 ) / 2 );
  lambda(1:m,1)   = +1.0;
  lambda(m+1:n,1) = -1.0;

  return
end


text:
function [vdiv,vcurl]=crldivxyz(vx,vy,vz)
% [vdiv,vcurl]=crldivxyz(vx,vy,vz)
% computes the divergence and curl of a
% vector expressed in x,y,z coordinates
syms x y z real
if ischar(vx), vx=sym(vx); end
if ischar(vy), vy=sym(vy); end
if ischar(vz), vz=sym(vz); end
vdiv=diff(vx,x)+diff(vy,y)+diff(vz,z);
cx=diff(vz,y)-diff(vy,z);
cy=diff(vx,z)-diff(vz,x);
cz=diff(vy,x)-diff(vx,y);
vcurl=[cx;cy;cz];


text:
%% S2FunHarmonicSym
%
% <S2FunHarmonicSym.S2FunHarmonicSym |S2FunHarmonicSym|> is an extension of the <S2FunHarmonic.S2FunHarmonic |S2FunHarmonic|> class which stores
% the symmetry with the function.
%
%
%% Defining a S2FunHarmonic
%
%%
% *Definition via symmetrisation*
%
% The simplest way to define a <S2FunHarmonicSym.S2FunHarmonicSym |S2FunHarmonicSym|> is through the
% symmetrisation of an ordinary <S2FunHarmonic.S2FunHarmonic |S2FunHarmonic|>.

sF = S2Fun.smiley
cs = crystalSymmetry('432');

sFs1 = symmetrise(sF, cs);

%%
% * this symmetrises the function and gives back the result with the
% symmetry attached
%

plot(sFs1)

%%
% * Note that only the important part with respect to the symmetry is plotted
% * you can plot the full sphere using the argument |'complete'|

%%
% *Definition via function handle*
%
% If you have a function handle for the function you could create a
% |S2FunHarmonicSym| via quadrature. At first lets define a symmetry and a
% function handle which takes <vector3d.vector3d.html |vector3d|> as an
% argument and returns double:

f = @(v) 0.1*(v.theta+sin(8*v.x).*sin(8*v.y));
cs = crystalSymmetry('6/m');

%%
% Now you can call the quadrature command to get |sFs2| of type |S2FunHarmonicSym|
sFs2 = S2FunHarmonicSym.quadrature(f, cs)

%% Visualization
% The plot commands for a |S2FunHarmonicSym| by default plot the function only on the fundamenta Sector of the symmetry.
% E.g. the default |plot|-command look as follows

plot(sFs1);

%%
% Another Example is the contour plot

contour(sFs2);

author: mtex-toolbox

repo: mtex

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


text:
function mrtrix_itkclass_wmmask(classfile, fname)

% Convert the white matter segmentation file in itkGray (mrVista) into
% binary white matter mask in mrTrix.
%
% In itkGray (mrVista), white matter is described as 3 and 4 in the matrix
% in nifti file. This code converts itkGray segmentation into binary white
% matter mask (white matter:1, others:0) in order to use it in fiber tractography.
%
% INPUT:
% classfile: A full path to nifti file storing segmentation information in itkGray format
% fname: The name of the output file (does not require filename extension)
%
% EXAMPLE:
% classfile = 't1_class.nii.gz';
% fname = 't1_class_binary';
% mrtrix_itkclass_wmmask(classfile, fname)
% (C) Hiromasa Takemura, CiNet/Stanford VISTA Team, 2015

if notDefined('fname')
    fname = 't1_class_binary';
end

% Set savefile name
niftiname = [fname '.nii.gz'];
mifname = [fname '.mif'];

% Load file
classseg = niftiRead(classfile);

% Make new nifti structure and set file name
nii = classseg;
nii.fname = niftiname;

%% Make binary white matter mask.
% In ITKGray format, white matter is described as 3 and 4. Here, we put
% 1 for all white matter voxels and 0 for all other voxels.

% Set zero for all voxels
nii.data = zeros(size(classseg.data));

% Set one for white matter voxels
nii.data(classseg.data == 3) = 1;
nii.data(classseg.data == 4) = 1;

% Save binary nifti file
niftiWrite(nii);

%% Convert nifti to .mif format
mrtrix_mrconvert(niftiname, mifname);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function res = eq(p,q)
%EQ           Implements  p == q  for polynomials
%
%Result 1 iff p and q are mathematically identical
%

% written  08/28/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  p = removevars(permvars(polynom(p),'lex'));
  q = removevars(permvars(polynom(q),'lex'));
  res = isequal(p.e,q.e) & isequal(p.c,q.c) & isequal(p.v,q.v);

repo: OpenMAS


text:
function nbytes = slstructsize(S, fields, types)
%SLSTRUCTSIZE Compute the size of a struct
%
% $ Syntax $
%   - nbytes = slstructsize(S)
%   - nbytes = slstructsize(S, fields)
%   - nbytes = slstructsize(S, fields, types)
%
% $ Arguments $
%   - nbytes:            the number of bytes in the specified parts of the struct
%   - fields:            the fields that are counted in the computation
%   - types:             the types of the fields when output
%
% $ Description $
%
%   - nbytes = slstructsize(S) computes the bytes of the struct occupied
%     if S is a struct-array, then the resultant nbytes is sum of all sizes
%     of struct entries.
%
%   - nbytes = slstructsize(S, fields) only computes the specified fields.
%     If some specified fields are not in S, an error will be raised.
%
%   - nbytes = slstructsize(S, fields, types) the computation will
%     use the types for computing the size instead of using original types.
%
% $ History $
%   - Created by Dahua Lin on Dec 10th, 2005
%

%% parse and verify input arguments
if ~isstruct(S)
    error('S is not a struct');
end
if nargin < 2 || isempty(fields)
    fields = fieldnames(S);
end
if nargin < 3 || isempty(types)
    use_origin_types = true;
else
    use_origin_types = false;
end

%% compute
n = numel(S);
if n == 1

    nbytes = 0;
    nterms = length(fields);

    for t = 1 : nterms
        cur_term = S.(fields{t});
        if isempty(cur_term)
            continue;
        elseif isstruct(cur_term)
            cur_nbytes = slstructsize(cur_term);
        else
            if use_origin_types || isempty(types{t})
                cur_elem_bytes = sltypesize(class(cur_term));
            else
                cur_elem_bytes = sltypesize(types{t});
            end
            cur_nbytes = cur_elem_bytes * numel(cur_term);
        end

        nbytes = nbytes + cur_nbytes;
    end

else

    nbytes = 0;
    for i = 1 : n
        if use_origin_types
            nbytes = nbytes + slstructsize(S(i), fields);
        else
            nbytes = nbytes + slstructsize(S(i), fields, types);
        end
    end

end

author: lmthang


text:
function val = minc_variable(hdr,var_name,att_name)
% Read a MINC variable inside a header
% VAL = MINC_VARIABLE( HDR , [VAR_NAME] , [ATT_NAME] )
%
% HDR (structure) a minc header read with MINC_READ
% VAR_NAME (string, optional) the name of a variable.
% ATT_NAME (string, optional) the name of an attribute.
% VAL (various types)
%   if VAR_NAME and ATT_NAME are both specified, VAL is the value of the
%      specified variable/attribute.
%   If only VAR_NAME is specified, VAL is the list of attributes of the
%      variable (cell of strings).
%   If VAR_NAME is unspecified, VAL is the list of variables (cell of strings).
%
% Example:
% [hdr,vol] = minc_read('my_vol.mnc');
% dcosinesy = minc_variable(hdr,'yspace','direction_cosines');
%
% Maintainer : pierre.bellec@criugm.qc.ca
% See licensing information in the code.

% Copyright (c) Pierre Bellec, Centre de recherche de l'institut de
% gériatrie de Montréal, Département d'informatique et de recherche
% opérationnelle, Université de Montréal, 2013-2014.
% Keywords : medical imaging, I/O, reader, minc
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.

if (nargin<1)||~isstruct(hdr)||~isfield(hdr,'details')||~isfield(hdr.details,'variables')||~isfield(hdr.details,'globals')
    error('Please specify a valid HDR structure, see MINC_READ')
end

hdr = hdr.details;

list_var = {hdr.variables(:).name};

if nargin == 1
    val = list_var;
    return
end

if nargin >= 2
    ind = find(ismember(list_var,var_name));
    if isempty(ind)
        error('Could not find variable %s in HDR',var_name)
    end
    ind = ind(1);
end

varminc = hdr.variables(ind);
list_att = varminc.attributes;
if nargin == 2
    val = list_att;
    return
end

ind2 = find(ismember(list_att,att_name));
if isempty(ind2)
    error('Could not find attribute %s in variable %s',att_name,var_name)
end
val = varminc.values{ind2};

author: brainstorm-tools


text:
function pointsOut = applyHomography(pointsIn, H)
%% Apply homography transformation to a set of points.

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%

% convert to homogenious coord and transpose
pointsInHomogen = [pointsIn, ones(size(pointsIn,1),1)]';

% apply homography
pointsOutHomogen = H * pointsInHomogen;

% divide by w and convert back to non-homogenious, then transpose
pointsOut = [pointsOutHomogen(1,:)./pointsOutHomogen(3,:); ...
    pointsOutHomogen(2,:)./pointsOutHomogen(3,:)]';
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function [W,fun,iter] = wLassomtl(X,y,XtX,Xty,lambda,W0,tol,maxiter,L,Lflag)
% Lasso Multi-Task Learning
% min_W ||XW - Y||_F^2 + \sum_i lambda_ij*|W^ij|
% By Nesterov's method
% ------------------------------- Input ---------------------------------------------
% X: block diagonal data matrix whose i-th block (n_i*d) is the data matrix of the i-th task
% y: \sum_i{n_i}*1 vector; response vector which is stacked by the reponses of all tasks
% XtX: X'*X
% Xty: X'*y
% lambda: regularized parameter (vector)
% W0: d*m matrix; starting point of W
% tol: stopping tolerance
% maxiter: maximum iterative steps
% Lflag: estimate the upper bound of Lipschitz constant if nonzero, zero otherwise
%
% ------------------------------- Output -------------------------------------------
%
% W: output weight
% fun: function values
% iter: iterative steps
%
% -----------------------------------------------------------------------------------

W = W0;
[d,m] = size(W); % d: dimension, m: the number of tasks

Wn = W;
t_new = 1;
fun = zeros(maxiter+1,1);

% Initial function value
fun(1) = norm(X*W(:) - y)^2 + wL1norm(W,lambda);
count = 0;
for iter = 1:maxiter
    W_old = W;
    t_old = t_new;
    gradvec = 2*(XtX*Wn(:) - Xty);
    gradmat = reshape(gradvec,d,m);
    % If we estimate the upper bound of Lipschitz constant, no line search
    % is needed.
    if Lflag
        W = proximalwL1norm(Wn - gradmat/L, lambda/L);
    else
        % line search
        for inneriter = 1:20
            W = proximalwL1norm(Wn - gradmat/L, lambda/L);
            dW = W - Wn;
            if 2*(dW(:)'*XtX*dW(:)) <= L*sum(sum((dW.*dW)))
                break;
            else
                L = L*2;
            end
        end
    end
    fun(iter+1) = norm(X*W(:) - y)^2 + wL1norm(W,lambda);
    % stopping condition
    if abs(fun(iter) - fun(iter+1))/fun(iter+1) < tol
        count = count + 1;
    else
        count = 0;
    end
    if count >= 1
        break;
    end

    % Update the coefficient
    t_new = (1+sqrt(1+4*t_old^2))/2;
    Wn = W + (t_old-1)/t_new*(W - W_old);

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
% Copyright (c) 2015, Francesco Perrone
% 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.

function [WF,WFtower,t,dy,dz,Y,Z,Zbottom,Ztower] = turbulent_wind_field_generator(U0,I0,Seed,HubHt,Ny,Nz,Ly,Lz,dt,T,xLu,xLv,xLw,Lc,a,shearExp)

%% FREQUENY/TIME PREALLOCATION & INITIALISATION
T = round(T + (Ly/U0));
Nt = Ny*Nz; % Total number of grid points
N = find_steps_number(T,dt); % Number of samples (in frequency and time domain)
Fs = N/T;
df = Fs/N; % Frequency sample
fp = df:df:0.5*Fs; % Positive frequency array [Hz]
dT = T/N; % Wind file time step

%% DEFINITION OF SIMULATION GRID

[Y,Z,dy,dz] = define_grid_coordinates(Ly,Lz,Ny,Nz,HubHt); % Generate (Y,Z) coordinates of grid points and spacing on y-/z-axis
Zunique = unique(Z);
Zbottom = min(Zunique(:));
Ztower = (Zbottom):-dz:0;
% Ytower = zeros(size(Ztower));
Ntower = numel(Ztower);
Uvert = repmat(U0.*((Zunique/HubHt).^shearExp),1,Ny)'; % Exponential vertical wind profile
Uvert_tower = U0.*((Ztower/HubHt).^shearExp);
Uvert1 = Uvert(:);
Distance = abs(sqrt(bsxfun(@plus,sum([Y(:) Z(:)].^2,2),sum([Y(:) Z(:)].^2,2)') - 2*([Y(:) Z(:)]*[Y(:) Z(:)]'))); % Distance matrix

%% POWER SPECTRA

r = 5/3; % Kaimal spectrum exponent
sigma_u = (I0*0.01*U0);
sigma_v = 0.8*sigma_u;
sigma_w = 0.5*sigma_u;
var_u = sigma_u^2;
var_v = sigma_v^2;
var_w = sigma_w^2;

Suu_p = ((4*var_u*(xLu/U0))./((1 + 6*fp*(xLu/U0)).^r))*0.5; % Positive amplitude spectrum (U-component)
Svv_p = ((4*var_v*(xLv/U0))./((1 + 6*fp*(xLv/U0)).^r))*0.5; % Positive amplitude spectrum (V-component)
Sww_p = ((4*var_w*(xLw/U0))./((1 + 6*fp*(xLw/U0)).^r))*0.5; % Positive amplitude spectrum (W-component)

Su = sqrt(Suu_p*N*Fs);
Sv = sqrt(Svv_p*N*Fs);
Sw = sqrt(Sww_p*N*Fs);

%% CALCULATION OF FFT TERMS

rng(Seed,'twister'); % Start Mersenne twister

fft_uu_p = zeros(Nt,N*0.5);
fft_vv_p = zeros(Nt,N*0.5);
fft_ww_p = zeros(Nt,N*0.5);

nn_u = exp(1i*2*pi*rand(Nt,numel(fp)));
nn_v = exp(1i*2*pi*rand(Nt,numel(fp)));
nn_w = exp(1i*2*pi*rand(Nt,numel(fp)));

for idx_chol = 1:numel(fp)
    Coh_uu = chol(exp(-a.*Distance.*sqrt(((fp(idx_chol)./U0).^2) + ((0.12./Lc).^2))),'lower'); %Coherence function U-component (according to IEC 61400-1 Ed.3)
    fft_uu_p(:,idx_chol) = (Coh_uu*Su(idx_chol))*nn_u(:,idx_chol);%chol(Coh_uu,'lower')
    fft_vv_p(:,idx_chol) = (Coh_uu*Sv(idx_chol))*nn_v(:,idx_chol);%chol(Coh_vv_ww,'lower')
    fft_ww_p(:,idx_chol) = (Coh_uu*Sw(idx_chol))*nn_w(:,idx_chol);%chol(Coh_vv_ww,'lower')
end

% Prepare FFT terms for each time series by mirroring the positive side
% about the frequency axis
fft_uu = [zeros(Nt,1) fft_uu_p fliplr(conj(fft_uu_p(:,1:end-1)))];
fft_vv = [zeros(Nt,1) fft_vv_p fliplr(conj(fft_vv_p(:,1:end-1)))];
fft_ww = [zeros(Nt,1) fft_ww_p fliplr(conj(fft_ww_p(:,1:end-1)))];

fft_uu(:,N*0.5 + 1) = real(fft_uu(:,N*0.5 + 1));
fft_vv(:,N*0.5 + 1) = real(fft_vv(:,N*0.5 + 1));
fft_ww(:,N*0.5 + 1) = real(fft_ww(:,N*0.5 + 1));

%% GENERATE TIME SERIES

Ucomp = (ifft((fft_uu),[],2));
Vcomp = (ifft((fft_vv),[],2));
Wcomp = (ifft((fft_ww),[],2));

%% SCALE TIME SERIES

stdU = (std(Ucomp,[],2));
stdV = (std(Vcomp,[],2));
stdW = (std(Wcomp,[],2));

idx_hub = ((Nz+1)*Ny)*0.5 - (Ny-1)*0.5;

SF = [sigma_u./stdU(idx_hub) , sigma_v./stdV(idx_hub) , sigma_w./stdW(idx_hub)];

Ucomp1 = arrayfun(@(i) Ucomp(i,:)*SF(1),1:Nt,'UniformOutput',false);
Ucomp1 = vertcat(Ucomp1{:});
Vcomp1 = arrayfun(@(i) Vcomp(i,:)*SF(2),1:Nt,'UniformOutput',false);
Vcomp1 = vertcat(Vcomp1{:});
Wcomp1 = arrayfun(@(i) Wcomp(i,:)*SF(3),1:Nt,'UniformOutput',false);
Wcomp1 = vertcat(Wcomp1{:});

t = 0:dT:(N-1)*dT;

WF = zeros(N,Nz,Ny,3);
WFtower = zeros(N,numel(Ztower),3);
it = 0;
for iz = 1:Nz
    for iy = 1:Ny
        it = it + 1;
        WF(:,iz,iy,1) = Ucomp1(it,:) + Uvert1(it);
        WF(:,iz,iy,2) = Vcomp1(it,:);
        WF(:,iz,iy,3) = Wcomp1(it,:);
    end
end

for iz = 1:Ntower
        WFtower(:,iz,1) = (WF(:,1,.5*(Ny+1),1) - mean(WF(:,1,.5*(Ny+1),1))) + Uvert_tower(iz);
        WFtower(:,iz,2) = WF(:,1,.5*(Ny+1),2);
        WFtower(:,iz,3) = WF(:,1,.5*(Ny+1),3);
end

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function CropImage(varargin)
% CROPIMAGE Taser wrapper for chipfile.m function
%
% Written by: Tim Cox, NRL
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

p = inputParser; % Extract parameter-value pairs
p.KeepUnmatched=true;
p.addParamValue('filename','');
p.addParamValue('aoi',[0 0 0 0]);
p.addParamValue('segment',1);
p.parse(varargin{:});

filenames = p.Results.filename;
aoi = p.Results.aoi;
segment = p.Results.segment;

%get output folder
%load last path
if ispref('matlab_sar_toolbox','last_used_directory')
    pathstr = getpref('matlab_sar_toolbox','last_used_directory');
    if ~ischar(pathstr)||~exist(pathstr,'dir')
        pathstr = pwd;
    end
else
    pathstr = pwd;
end

outdir = uigetdir(pathstr,'Specify Output Folder');

if iscell(filenames)
    for i=1:length(filenames)
        %construct output filename
        [pathstr, name, ext] = fileparts(filenames{i});
        aoistr = sprintf('_%d_%d_%d_%d_seg_%d',aoi(1),aoi(2),aoi(3),aoi(4),segment);
        outname = [outdir filesep name aoistr '.nitf'];

        %crop image
        chipfile(filenames{i},outname,'azlimits',[aoi(1) aoi(1)+aoi(3)-1],'rnglimits',[aoi(2) aoi(2)+aoi(4)-1],...
                 'framenumber',segment,'file_format','SICD','block_size',2^21);
    end
else
    %construct output filename
    [pathstr, name, ext] = fileparts(filenames);
    aoistr = sprintf('_%d_%d_%d_%d_seg_%d',aoi(1),aoi(2),aoi(3),aoi(4),segment);
    outname = [outdir filesep name aoistr '.nitf'];

    %crop image
    chipfile(filenames,outname,'azlimits',[aoi(1) aoi(1)+aoi(3)-1],'rnglimits',[aoi(2) aoi(2)+aoi(4)-1],...
             'framenumber',segment,'file_format','SICD','block_size',2^21);
end

msgbox('Image Crop Completed');

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function img = read_pfm(filename, bPhotoshopCompatibility)
%
%
%       img = read_pfm(filename, bPhotoshopCompatibility)
%
%        Input:
%           -filename: the name of the file to open
%           -bPhotoshopCompatibility: a flag for enabling compability with
%           Adobe Photoshop. If it is set to 1 the compability is on.
%
%        Output:
%           -img: the opened image
%
%     Copyright (C) 2011-14  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist('bPhotoshopCompatibility', 'var'))
    bPhotoshopCompatibility = 0;
end

%by default PFM files are open in little-endian format
fid = fopen(filename, 'r', 'l');

%reading the header
str = fscanf(fid, '%c', 3);

channels = 3;
iValid = 0;

if(str(1) == 'P')
    iValid = 1;
end

if(str(2) == 'f')
    channels = 1;
    iValid = iValid + 1;
end

if(str(2) == 'F')
    iValid = iValid + 1;
end

if(iValid ~= 2)
    fclose(fid);
    error('This is not a valid PFM file.');
end

%image size
m = fscanf(fid,'%d',1);
fscanf(fid,'%c',1);
n = fscanf(fid,'%d',1);

fscanf(fid,'%c',1);
endian_selector = fscanf(fid,'%f',1);
fscanf(fid,'%c',1);

%we have a big-endian encoded file
if(endian_selector > 0.0)
    fclose(fid);

    %reopen the file in big-endian mode
    fid = fopen(filename, 'r', 'b');
    fscanf(fid,'%c',3);
    m = fscanf(fid,'%d',1);
    fscanf(fid,'%c',1);
    n = fscanf(fid,'%d',1);
    fscanf(fid,'%c',1);
    endian_selector = fscanf(fid,'%f',1);
    fscanf(fid,'%c',1);
end

img = zeros([m, n, channels]);

total = n * m * channels;
tmpImg = fread(fid, total, 'float');

for i=1:channels
    tmpC = i:channels:total;
    img(:,:,i) = reshape(tmpImg(tmpC), m, n);
end

img = imrotate(img, 90, 'nearest');

if(bPhotoshopCompatibility)
    for i=1:channels
        img(:,:,i) = flipud(img(:,:,i));
    end
end

fclose(fid);

end

author: banterle


text:
%% Simulating Pole Figure data
%
%%
% Simulating pole figure data from a given ODF is useful to investigate
% pole figure to ODF reconstruction routines. Let us start with a model ODF
% given as the superposition of 6 components.

cs = crystalSymmetry('orthorhombic');
mod1 = orientation.byAxisAngle(xvector,45*degree,cs);
mod2 = orientation.byAxisAngle(yvector,65*degree,cs);
model_odf = 0.5*uniformODF(cs) + ...
  0.05*fibreODF(Miller(1,0,0,cs),xvector,'halfwidth',10*degree) + ...
  0.05*fibreODF(Miller(0,1,0,cs),yvector,'halfwidth',10*degree) + ...
  0.05*fibreODF(Miller(0,0,1,cs),zvector,'halfwidth',10*degree) + ...
  0.05*unimodalODF(mod1,'halfwidth',15*degree) + ...
  0.3*unimodalODF(mod2,'halfwidth',25*degree);

%%

plot(model_odf,'sections',6,'silent','sigma')

%%
% In order to simulate pole figure data, the following parameters have to be
% specified
%
% # an arbitrary <SO3Fun.SO3Fun.html ODF>
% # a list of <Miller.Miller.html Miller indece>
% # a grid of <S2Grid.S2Grid.html specimen directions>
% # superposition coefficients (optional)
% # the magnitude of error (optional)
%
%%
% The list of <Miller.Miller.html Miller indece>

h = [Miller(1,1,1,cs),Miller(1,1,0,cs),Miller(1,0,1,cs),Miller(0,1,1,cs),...
  Miller(1,0,0,cs),Miller(0,1,0,cs),Miller(0,0,1,cs)];

%%
% The <S2Grid.S2Grid.html grid> of specimen directions

r = regularS2Grid('resolution',5*degree);

%%
% Now the pole figures can be simulated using the command
% <SO3Fun.calcPoleFigure.html calcPoleFigure>.

pf = calcPoleFigure(model_odf,h,r)

%%
% Add some noise to the data. Here we assume that the mean intensity is 1000.

pf = noisepf(pf,1000);

%%
% Plot the simulated pole figures.

plot(pf)

%% ODF Estimation from Pole Figure Data
%
% From these simulated pole figures we can now estimate an ODF,

odf = calcODF(pf)

%%
% which can be plotted,

plot(odf,'sections',6,'silent','sigma')

%%
% and compared to the original model ODF.

calcError(odf,model_odf,'resolution',5*degree)

%% Exploration of the relationship between estimation error and number of pole figures
%
% For a more systematic analysis of the estimation error, we vary the number
% of pole figures used for ODF estimation from 1 to 7 and calculate for any
% number of pole figures the approximation error. Furthermore, we also
% apply ghost correction and compare the approximation error to the
% previous reconstructions.

e = [];
for i = 1:pf.numPF

  odf = calcODF(pf({1:i}),'silent','NoGhostCorrection');
  e(i,1) = calcError(odf,model_odf,'resolution',2.5*degree);
  odf = calcODF(pf({1:i}),'silent');
  e(i,2) = calcError(odf,model_odf,'resolution',2.5*degree);

end

%%
% Plot the error in dependency of the number of single orientations.

close all;
plot(1:pf.numPF,e,'LineWidth',2)
ylim([0.07 0.32])
xlabel('Number of Pole Figures');
ylabel('Reconstruction Error');
legend({'Without Ghost Correction','With Ghost Correction'});


text:
function d = NormGroupL2_dual(groups,x,weights)

if isreal(x)
   d = norm(sqrt(sum(groups * x.^2,2))./weights,inf);
else
   d = norm(sqrt(sum(groups * abs(x).^2,2))./weights,inf);
end

author: mpf


text:
function sphere_triangle_quad_test ( )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_QUAD_TEST tests the SPHERE_TRIANGLE_QUAD library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_TRIANGLE_QUAD_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SPHERE_TRIANGLE_QUAD library.\n' );

  sphere_triangle_quad_test01 ( );
  sphere_triangle_quad_test02 ( );
  sphere_triangle_quad_test03 ( );
  sphere_triangle_quad_test04 ( );
  sphere_triangle_quad_test05 ( );
  sphere_triangle_quad_test06 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_TRIANGLE_QUAD_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 function xs = l1_tv_restore1_fun(yi, A, varargin)
%function xs = l1_tv_restore1_fun(yi, A, varargin)
%|
%| Robust image restoration using a l1 data fit term and a TV-like regularizer.
%|
%| see l1_tv_restore1.m
%| min_x |yi - A(ti) * x|_1 + beta (|Cx x|_1 + |Cy x|_1)
%| where |.|_1 is approximated by a hyperbola
%|
%|
%| in
%|	yi	data
%|	A	system matrix
%|
%| option
%|	several - see code ...
%|
%| out
%|	xs	estimate(s)
%|
%| Copyright 2010-05-24, Jeff Fessler, University of Michigan

if nargin == 1 && streq(yi, 'test')
	l1_tv_restore1
	clear xs
return
end

if nargin < 2, help(mfilename), error(mfilename), end

arg.niter = 20;
arg.delta1 = 0.2; % round corner of |t| in data fit to approximate by hyperbola
arg.delta2 = 0.2; % round corner of |t| in penalty
arg.xinit = [];
arg.mask = [];
arg.l2b = -6;
arg.curvtype = 'oc';
arg.args = {};
arg = vararg_pair(arg, varargin);

if isempty(arg.xinit)
	arg.xinit = yi; % assume restoration
end

if isempty(arg.mask)
	arg.mask = true(size(yi));
end

% R = R_null; % null regularizer
R = Reg1(arg.mask, 'type_denom', 'matlab', ...
	'beta', 2^arg.l2b, 'pot_arg', {'hyper2', arg.delta2});

data = {yi(:), arg.delta1};
xs = pl_pcg_qs_ls(arg.xinit(:), A, data, ...
	@l1_tv_dercurv, R, 'niter', arg.niter, ...
	'curvtype', arg.curvtype, ...
	arg.args{:});

% l1_tv_dercurv()
function [deriv curv] = l1_tv_dercurv(data, yp, curvtype)
yi = data{1};
delta = data{2};
pot = potential_fun('hyper2', delta);
t = yp - yi;
switch curvtype
case 'pc'
	curv = 1;
	deriv = pot.dpot(t);
case 'oc'
	curv = pot.wpot(t);
	deriv = curv .* t;
otherwise
	fail('curvtype %s', curvtype)
end

% R_null()
function R = R_null;
R.dercurv = @R_null_dercurv;
R.C1 = 0;

% R_null_dercurv()
function [a b] = R_null_dercurv(arg1, arg2);
a = 0;
b = 0;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function UID = getIOD_UID(IOD_type)
%getIOD_UID
%   Given a string indicating an IOD type, returns the UID corresponding to
%   that IOD.
%
%   Currently supported UID strings are: CT
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

switch upper(IOD_type)
    case 'CT'
        UID = '1.2.840.10008.5.1.4.1.1.2';
    case 'MR'
        UID = '1.2.840.10008.5.1.4.1.1.4';
    otherwise
        error('Unsupported IOD type.');
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function HaarTransformationMatrix=ConstructHaarMat(WidthOfSquareMatrix)
%--------------------------------------------------------------------------
%Create Haarwavelet transformation matrix H for the matrix vector
%mulplication implimentation of Haar wavelet transformation.
%This function uses the following nice formula to create the Haar
%transformation matrix:
%               H_n=1/sqrt(2)[H_(n/2) kron (1 1)
%                             I_(n/2) kron (1 -1)],
%                              where 'kron' denotes the kronecker product.
%The iteration starts with H_1=[1]. The normalization constant 1/sqrt(2)
%ensure that H_n^T*H_n=I, where I is identity matrix. Haar wavelets are the
%rows of H_n.
%--------------------------------------------------------------------------
%function HaarTransformationMatrix=ConstructHaarWaveletTransformationMatrix(WidthOfSquareMatrix)
% Input:
%       WidthOfSquareMatrix: the width of sqaure Haar wavelet
%                                  transformation matrix.
% Output:
%       HaarTransformationMatrix: Ceated Haar transformation matrix and
%       it's size is the power of 2,i.e., 2, 4,
%       8,16,32,64,128,256,512,1024,2048,4096,etc.
%--------------------------------------------------------------------------
%Author: Jin Qi
%Email: jqichina@hotmail.com
%Date:  11/4/2011
%--------------------------------------------------------------------------
%Test function for one dimentional signal
%
% n=8;
% S=rand(n,1); % for one dimentional signal
% H=ConstructHaarWaveletTransformationMatrix(n);
% C=H*S; %Haar wavelet transformation
% RS=H'*C;% Reconstruct signal
% fprintf('reconstruction error is %f',norm(S-RS));%plot reconstruction error
%--------------------------------------------------------------------------
%Test for two dimentional signal,ie. image
%
% n=8;
% S=rand(n,n); % for one dimentional signal
% H=ConstructHaarWaveletTransformationMatrix(n);
% C=H*S*H'; %Haar wavelet transformation
% RS=H'*C*H;% Reconstruct signal
% fprintf('reconstruction error is %f',norm(S-RS,'fro'));%plot reconstruction error

n=WidthOfSquareMatrix; % copy the parameter

% check input parameter and make sure it's the power of 2
Level=log2(n);
if 2^Level<n, error('please ensure the value of input parameter is the power of 2');end

%Initialization
H=[1];
NC=1/sqrt(2);%normalization constant
LP=[1 1];
HP=[1 -1];

% iteration from H=[1]
for i=1:Level
    H=NC*[kron(H,LP);kron(eye(size(H)),HP)];
end
HaarTransformationMatrix=H;

author: thomaskuestner


function [  ] = DP_save_rows( cutCostsFile, cut_num,  search_k, saveFileAddress )
%DP_SAVE_ROWS Summary of this function goes here
%   Detailed explanation goes here

    file_t = fopen(cutCostsFile);
    cut_costs = fscanf(file_t, '%d');
    fclose(file_t);

    % get the initial averaging-cut

    cut_length = length(cut_costs) / cut_num;

    for i = 1:cut_num-1
       initial_cuts(i,:) = floor(i*cut_length);
    end

    % start Dynamic Programming
    % find the best/min value in the domain near cut_points as initial
    % values

    best_values = [];
    best_cuts = [];
    search_domain_start = floor(cut_length / search_k);
    search_domain_finish = floor(cut_length / search_k);

    sum_best_value = 0;

    for i = 1:length(initial_cuts)

        start = initial_cuts(i,:)-search_domain_start;
        finish = initial_cuts(i,:)+search_domain_finish;

        if start <= 0
            start = 1;
            search_domain_start = initial_cuts(i,:) - start;
        end

        if finish > length(cut_costs)-1
            finish  = length(cut_costs)-1;
            search_domain_finish = finish - initial_cuts(i,:);
        end

        domain = cut_costs(start:finish,:);

        % observe the minimum
%         [min_value, min_cut] = min(domain);
        % obserbe the maximum
        [min_value, min_cut] = max(domain);

        best_values(i,:) = min_value;

        if min_cut < initial_cuts(i,:)
            best_cuts(i,:) = initial_cuts(i,:) + (min_cut - search_domain_start - 1);
        else
            best_cuts(i,:) = initial_cuts(i,:) + (min_cut - search_domain_finish - 1);
        end

        if i > 1 && best_cuts(i,:) == best_cuts(i-1,:)
            best_cuts(i,:) = best_cuts(i,:) + 1;
            best_values(i,:) = domain(min_cut + 1);
        end

        sum_best_value = sum_best_value + best_values(i,:);

    end

    dlmwrite(saveFileAddress, best_cuts);
    sum_best_value

end


text:
function S3G = delete(S3G,ind)
% clear nodes by indece
% Input
%  SOG    - @SO3Grid
%  indece - int32
%
% Output
%   SO3Grid - "not indexed"

% remove orientations
S3G.a(ind) = [];
S3G.b(ind) = [];
S3G.c(ind) = [];
S3G.d(ind) = [];
S3G.i(ind) = [];

S3G.gamma = subGrid(S3G.gamma,~ind);
S3G.alphabeta  = subGrid(S3G.alphabeta,GridLength(S3G.gamma)>0);
S3G.gamma(GridLength(S3G.gamma)==0) = [];


text:
% Function to calculate H2-H1 and creak-F0 using resonator based method
%
% Description
%  Function to calculate H2-H1 and creak-F0 using resonator based method
%
%  This version is a slightly later one than that described in the
%  above published 2013 CSL paper [2]. The algorithm here rather than using binary
%  decision trees using artificial neural networks and combines the
%  features used in the CSL paper with those proposed in Ishi et al.
%  (2008). This updated version has been submitted to CSL for a special
%  issue on glottal source processing on April 14th 2013. It will have
%  reference [1].
%
% Inputs
%  res      : [samples] [Nx1] Linear prediction residual
%  fs       : [Hz]      [1x1] sampling frequency
%  F0mean   : [Hz]      [1x1] fundamental frequency
%
% Outputs
%  H2H1     : [dB]      [Mx2] H2-H1 parameter
%  f0       : [Hz]      [Mx1] Creak-F0 parameter
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% Octave Compatibility
%  As this algorithm relies on the Matlab neural networks toolbox it is not
%  compatible with Octave. Furthermore, this algorithm requires Matlab
%  R2011 or later
%
% References
%  [1] Drugman, T., Kane, J., Gobl, C., `Automatic Analysis of Creaky
%       Excitation Patterns', Submitted to Computer Speech and
%       Language.
%  [2] Kane, J., Drugman, T., Gobl, C., (2013) `Improved automatic
%       detection of creak', Computer Speech and Language 27(4), pp.
%       1028-1047.
%  [3] Drugman, T., Kane, J., Gobl, C., (2012) `Resonator-based creaky
%       voice detection', Interspeech 2012, Portland, Oregon, USA.
%
%
% Copyright (c) 2013 Trinity College Dublin
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  This program is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Authors
%  Thomas Drugman <thomas.drugman@umons.ac.be> & John Kane <kanejo@tcd.ie>
%

function [H2H1,f0] = get_creak_h2h1(res,fs,F0mean)

%% Initial settings
buffer=0;
Start=1;
Stop=50/1000*fs;

Hop=10/1000*fs;
Win=hanning(Stop);
Ind=1;
f0=[];
Delta=[];

moving_average_size=100/1000*fs;
moving_average_frame=moving_average_size/Hop;

% Resonator 2 settings
Phi=2*pi*1*F0mean/fs;
Rho=0.97;
rep=filter([1 0 0],[1 -2*Rho*cos(Phi) Rho^2],res);

rep=rep/max(abs(rep));

% Resonator 1 settings
Phi=2*pi*1*F0mean/fs;
Rho=0.8;
rep2=filter([1 0 0],[1 -2*Rho*cos(Phi) Rho^2],res);

%% Do processing
while Stop<length(rep)
    Sig=rep2(Start:Stop)';
    Sig=Sig.*Win;

    Corrs=xcorr(Sig,Sig);
    Corrs(1:length(Sig))=[];

    for k=1:length(Corrs)
       Corrs(k)=Corrs(k)*(length(Corrs)/(length(Corrs)-(k-1)));
    end

    Corrs2=Corrs;
    k=1;
    while (Corrs2(k)>0)
        Corrs2(k)=0;
        k=k+1;
        if k==length(Corrs2)
            break;
        end
    end
    [maxi,posi]=max(Corrs2);

    if posi<5
        posi=5;
    end

    F0=round(fs/posi);

    f0(Ind)=F0;
    F0=round(f0(Ind));

    Sig=rep(Start:Stop)';
    Sig=Sig.*Win;

    Corrs=xcorr(Sig,Sig);
    Corrs(1:length(Sig))=[]; Sig=Sig.*Win;

    Corrs=xcorr(Sig,Sig);
    Corrs(1:length(Sig))=[];

    Spec=fft(Corrs,fs);
%     Spec=fft(Sig,fs);
    Spec=abs(Spec(1:fs/2));

    Spec=Spec/sqrt(sum(Spec.^2));
    Spec=20*log10(Spec);

    Delta(Ind)=(Spec(2*F0)-Spec(F0))+buffer;

    Start=Start+Hop;
    Stop=Stop+Hop;
    Ind=Ind+1;
end

%% Interpolate to update every sample
if exist('Delta','var')
    if isempty(Delta)==0
        Delta=smooth(Delta,moving_average_frame);
        f0=medfilt1(f0,7);
        f0=smooth(f0,moving_average_frame);
        Delta=interp1(linspace(1,length(res),length(Delta)),Delta,1:length(res));
        %Delta=resample(Delta,length(res),length(Delta));
        Delta=Delta-buffer;

        f0=interp1(linspace(1,length(res),length(f0)),f0,1:length(res));

        %Delta2=smooth(Delta,moving_average_size); % 100ms moving average filter
        %f0=smooth(f0,moving_average_size); % 100ms moving average filter
        Delta2=Delta;

        H2H1=Delta2;
    else H2H1=zeros(1,length(res))-buffer;
        f0=zeros(1,length(res));
    end
else H2H1=zeros(1,length(res))-buffer;
    f0=zeros(1,length(res));
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [ pn, dist ] = circle_arc_point_near_2d ( r, center, theta1, ...
  theta2, p )

%*****************************************************************************80
%
%% CIRCLE_ARC_POINT_NEAR_2D : nearest point on a circular arc.
%
%  Discussion:
%
%    A circular arc is defined by the portion of a circle (R,C)
%    between two angles (THETA1,THETA2).
%
%    Thus, a point (X,Y) on a circular arc satisfies
%
%      ( X - C(1) ) * ( X - C(1) ) + ( Y - C(2) ) * ( Y - C(2) ) = R * R
%
%    and
%
%      Theta1 <= Theta <= Theta2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2,1), the center of the circle.
%
%    Input, real THETA1, THETA2, the angles defining the arc,
%    in radians.  Normally, THETA1 < THETA2.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, real PN(2,1), a point on the circular arc which is
%    nearest to the point.
%
%    Output, real DIST, the distance to the nearest point.
%
  dim_num = 2;
%
%  Special case, the zero circle.
%
  if ( r == 0.0 )
    pn(1:dim_num,1) = center(1:dim_num,1);
    dist = sqrt ( sum ( ( p(1:dim_num,1) - pn(1:dim_num,1) ).^2 ) );
    return
  end
%
%  Determine the angle made by the point.
%
  theta = atan4 ( p(2,1) - center(2,1), p(1,1) - center(1,1) );
%
%  If the angle is between THETA1 and THETA2, then you can
%  simply project the point onto the arc.
%
  if ( r8_modp ( theta  - theta1,  2.0 * pi ) <= ...
       r8_modp ( theta2 - theta1,  2.0 * pi ) )

    r2 = sqrt ( sum ( ( p(1:dim_num,1) - center(1:dim_num,1) ).^2 ) );

    pn(1:dim_num,1) = center(1:dim_num,1) ...
      + ( p(1:dim_num,1) - center(1:dim_num,1) ) * r / r2;
%
%  Otherwise, if the angle is less than the negative of the
%  average of THETA1 and THETA2, it's on the side of the arc
%  where the endpoint associated with THETA2 is closest.
%
  elseif ( r8_modp ( theta - 0.5 * ( theta1 + theta2 ), 2.0 * pi ) <= pi )

    pn(1:dim_num,1) = center(1:dim_num,1) + r * [ cos ( theta2 ), sin ( theta2 ) ];
%
%  Otherwise, the endpoint associated with THETA1 is closest.
  else

    pn(1:dim_num,1) = center(1:dim_num,1) + r * [ cos ( theta1 ), sin ( theta1 ) ];

  end

  dist = sqrt ( sum ( ( p(1:dim_num,1) - pn(1:dim_num,1) ).^2 ) );

  return
end


text:
function [results_table_pos, results_table_neg, r, excluded_region_table, r_excluded, region_list] = table_of_atlas_regions_covered(obj, varargin)
% Make a table of which atlas parcels are covered by a set of regions or map identified in a study
%
% [results_table_pos, results_table_neg, r, excluded_region_table, r_excluded, region_list] = table_of_atlas_regions_covered(obj, [atlas_obj])
%
% fMRI activation maps often include large suprathreshold areas that span
% multiple brain regions.  Traditional tables divide areas by contiguous
% regions ("blobs"), but this type of division is not useful if blobs cannot
% be easily summarized by a single anatomical label. A complementary approach
% is to identify which areas defined in a standard brain atlas ("parcels")
% are covered by the activation map. Parcels can be defined based on anatomy
% (e.g., cytoarchitecture) or function (e.g., resting-state fMRI).
% For example, a large blob may span the insula, claustrum, and putamen.
% This can be described in a table that lists all of these regions, along with
% the percentage of each parcel covered by the activation map.
%
% This function uses the object method subdivide_by_atlas() to create such
% a table, which is dividided into two tables with positive and negative average
% statistic values (results_table_pos, results_table_neg). These are Matlab
% table-class objects.  It also returns a region-class object for reference
% and rendering the subdivided blobs on brain slices or surfaces.
%
% You can enter any atlas-class object as a 2nd argument. If you don't,
% a default atlas object will be used.
%
% Examples:
% -----------------------------------------------------------------
% % Load a thresholded map:
% % (Note: You must have Neuroimaging_Pattern_Masks repo with subfolders on your path)
% % This is a predictive map for drug craving (Koban et al. 2022, Nat Neurosci)
%
% ncsthr = fmri_data(which('NCS_multithr_001k5_005_05_pruned.nii'), 'noverbose');
%
% % Convert to region object:
% r = region(ncsthr);
%
% [results_table_pos, results_table_neg, r, excluded_region_table, table_legend_text] = table_of_atlas_regions_covered(r);
% montage(r, 'regioncenters', 'colormap');
%
% % Also see large blobs that didn't sufficiently cover any one atlas parcel:
% montage(r_excluded(cat(1, r_excluded.numVox) > 20), 'regioncenters', 'colormap');
%
% % You can use the image_vector method to apply this to a single-image
% % thresholded map:
% [results_table_pos, results_table_neg, r, excluded_region_table, r_excluded, region_labels] = table_of_atlas_regions_covered(ncsthr);
%
% Note: For a key to labels when using the Glasser 2016 Nature atlas, or the
% CANlab combined atlas (which uses Glasser), see GlasserTableS1.pdf in the
% original paper or CANlab github, or http://braininfo.rprc.washington.edu/
%
% see also region.table, for autolabeling of regions with an atlas and more
% detailed table output.
%
% Tor Wager, Feb 2023

if (size(obj.dat, 2) > 1)
    error('Use this function with fmri_data or image_vector objects containing only a single image. Use get_wh_image() to select one.')
end

r = region(obj);

[results_table_pos, results_table_neg, r, excluded_region_table, r_excluded, region_list] = table_of_atlas_regions_covered(r, varargin{:});

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore