text:
function [opt, isdefault]= opt_setDefaults(opt, props, check)
%OPT_SETDEFAULTS - Set default values according to a property spec list
%
%Synopsis:
%  [OPT, ISDEFAULT]= opt_setDefaults(OPT, PROPSPEC)
%  [OPT, ISDEFAULT]= opt_setDefaults(OPT, PROPSPEC, CHECK)
%
%Arguments:
%  OPT:      STRUCT of optional properties
%  PROPSPEC: PROPSPECLIST - Property specification list, i.e., CELL of size
%      [N 2] or [N 3], with the first column all being strings.
%  CHECK:    BOOL - Evokes type checking according to the third column of
%      PROPSPEC
%
%Returns:
%  OPT: STRUCT with added properties from the property specification list
%          PROPSPEC
%  ISDEFAULT: STRUCT with the same fields as OPT. Each field has a
%          boolean value indicating whether the default value from PROPSPEC
%          was taken (i.e., that property was not a field of the input OPT).
%
%Description:
%  PROPSPEC is a property specification list that defines default values
%  for a list of properties. The struct OPT is filled with all missing
%  properties from PROPSPEC, but existing fields of OPT are not overwritten
%  by the respective default values.
%  The function makes a case-insensitive match of field of OPT to the
%  property names in PROPLIST. In the case of a match, the variant of
%  PROPSPEC is used.
%
%See also opt_checkProplist, opt_proplistToStruct.
%
%Examples:
%  opt= struct('linewidth',3, 'color','k')
%  props= {'LineWidth', 2;  'Color',[0 0 1];  'LineStyle','--'}
%  [opt, isdefault]= opt_setDefaults(opt, props)
%
%  % If the first input variable (OPT) is not a variable itself, but a field
%  % of a struct, it has to passed as a string:
%  S.opt= opt_setDefaults('S.opt', props);

% 06-2012 Benjamin Blankertz

misc_checkType(opt, 'STRUCT|CHAR');
misc_checkType(props, 'PROPSPEC');
misc_checkTypeIfExists('check', 'BOOL');

if nargin<3,
  check= 0;
end

if ischar(opt),
  structname= opt;
  opt= evalin('caller', opt);
  check= 1;
else
  structname= '';
end

% Perform a case correction in the STRUCT opt, and
% Set 'isdefault' to ones for the field already present in 'opt'
isdefault= [];
if ~isempty(opt),
  for Fld=fieldnames(opt)',
    fld= Fld{1};
    idx= find(strcmpi(fld, props(:,1)));
    if length(idx)>1,
      error('Case ambiguity in propertylist');
    elseif ~isempty(idx) && ~strcmp(fld, props{idx,1}),
      % rename field in opt to match case with props
      oldfld= fld;
      fld= props{idx,1};
      [opt.(fld)]= opt.(oldfld);
      opt= rmfield(opt, oldfld);
    end
    isdefault.(fld)= 0;
  end
end

for k= 1:size(props,1),
  fld= props{k,1};
  if ~isfield(opt, fld),
    %% if opt is a struct *array*, the fields of all elements need to
    %% be set. This is done with the 'deal' function.
    [opt.(fld)]= deal(props{k,2});
    isdefault.(fld)= 1;
  end
end

if check,
  if check==2,
    tcstate= bbci_typechecking('on');
  end
  if isempty(structname),
    structname= inputname(1);
  end
  opt_checkProplist(opt, props, structname);
  if check==2,
    bbci_typechecking(tcstate);
  end
end

author: PatternRecognition

repo: OpenBMI


text:
% VL_DEMO_SLICT  Demo: SLIC superpixels

prefix = fullfile(vl_root,'figures','demo') ;
randn('state',0) ;
rand('state',0) ;
figure(1) ; clf ;

% --------------------------------------------------------------------
%                                                        Load a figure
% --------------------------------------------------------------------

im = imread(fullfile(vl_root,'data','roofs1.jpg')) ;
im = im2single(im) ;
im = im(1:128,end-128+1:end,:) ;

figure(1) ; clf ;
image(im) ;
axis equal off tight ;
vl_demo_print('slic_image') ;

% --------------------------------------------------------------------
%                                    Create various SLIC segmentations
% --------------------------------------------------------------------

regionSizes = [10 30] ;
regularizers = [0.01 0.1 1] ;

figure(2) ; clf ;
for i = 1:numel(regionSizes)
  for j = 1:numel(regularizers)
    regionSize = regionSizes(i) ;
    regularizer = regularizers(j) ;
    segments = vl_slic(im, regionSize, regularizer, 'verbose') ;

    % overaly segmentation
    [sx,sy]=vl_grad(double(segments), 'type', 'forward') ;
    s = find(sx | sy) ;
    imp = im ;
    imp([s s+numel(im(:,:,1)) s+2*numel(im(:,:,1))]) = 0 ;

    vl_tightsubplot(numel(regionSizes),numel(regularizers), (i-1)*numel(regularizers) + j) ;
    imagesc(imp) ; axis image off ; hold on ;
    text(5,5,sprintf('regionSize:%.2g\nregularizer:%.2g', regionSize, regularizer), ...
         'Background', 'white','VerticalAlignment','top')
  end
end

vl_demo_print('slic_segmentation') ;

author: Cloud-CV

repo: object-proposals


text:
function model = mappingOptimise(model, X, Y, varargin)

% MAPPINGOPTIMISE Optimise the given model.

% MLTOOLS

fhandle = str2func([model.type 'Optimise']);
model = fhandle(model, X, Y, varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%tstoolbox/mex/boxcount
%   boxcount is a fast algorithm that partitions a data set of points into
%   equally spaced and sized boxes. The algorithm is based on Robert
%   Sedgewick's Ternary Search Trees which offer a fast and efficient way
%   to create and search a multidimensional histogram. Empty boxes require
%   no storage space, therefore the maximum number of boxes (and memory)
%   used can not exceed the number of points in the data set, regardless
%   of the data set's dimension and the number of partitions per axis.
%
%   During processing, data values are scaled to be within the range
%   [0,1]. All columns of the input matrix are scaled by the same factor,
%   so no skewing is introduced into the point set.
%
%   Syntax:
%
%     * [a,b,c] = boxcount(point_set, partitions)
%
%   Input arguments:
%
%     * pointset - a N by D double matrix containing the coordinates of
%       the point set, organized as N points of dimension D. D is limited
%       to 128.
%     * partitions - number of partitions per axis, limited to 16384. For
%       convenience, if a vector is given, boxcount will iterate over all
%       values of this vector.
%
%   Output arguments:
%
%     * a - vector of size D with: log2(sum(Number of nonempty boxes))
%     * b - vector of size D with: sum(p * log2(p)) , where p is the
%       relative frequency of points falling into a box
%     * c - vector of size D with: log2(sum(p*p)), where p is the relative
%       frequency of points falling into a box
%
%   Example:
%
% p = rand(50000, 4);
% p = p - min(min(p));
% p = p ./ max(max(p));
% [a,b,c] = boxcount(p, 16)
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function conf = voc_config_car_grammar()
% Set up configuration variables

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf.project = 'rel5-dev/car-grammar-3';

conf.training.C                    = 0.006;
conf.training.wlssvm_M             = 1;
conf.training.lbfgs.options.optTol = 0.0001;

conf.eval.max_thresh = -1.4;

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function Xr=pround(X,d)

Xr=round(X.*(10^d))./(10^d);

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function cS = rutile

%
cs_Rt = crystalSymmetry('4/mmm',[1 1 0.644]);
N = Miller({1,0,0},{1,1,0},{1,0,1},cs_Rt);
dist= [0.35, 0.5, 1];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
function ori = exp(m,ori_ref,varargin)
% misorientation vector to misorientation
%
% Description
%
% Syntax
%
%   mori = exp(m)        % misorientation
%
%   ori = exp(m,ori_ref) % orientation update
%
% Input
%  m - @Miller misorientation vector in crystal coordinates
%  ori_ref - @orientation
%
% Output
%  ori - @orientation
%
% See also
% Miller/exp orientation/log

rot = exp@vector3d(m);

ori = orientation(rot,m.CS,m.CS);

if nargin >= 2
  if check_option(varargin,'left')
    ori =  ori * ori_ref;
  else
    ori =  ori_ref .* ori;
  end
end


text:
function out = calcGrainVotes(job,varargin)
% compute votes from grain boundaries
%
% Syntax
%
%   % compute votes from all p2c and c2c boundaries
%   job.calcGBVotes('threshold', 2*degree)
%
%   % compute votes only from p2c boundaries -> growth algorithm
%   job.calcGBVotes('p2c', 'threshold', 3*degree, 'tol', 1.5*degree)
%
% Input
%  job - @parentGrainReconstructor
%
% Output
%  job.votes - table of votes
%
% Options
%  p2c  - consider only parent / child grain boundaries
%  c2c  - consider only child / child grain boundaries
%  threshold - threshold fitting angle between job.p2c and the boundary OR
%  tolerance - range over which the probability increases from 0 to 1 (default 1.5)
%  numFit - number of fits to be computed
%

% noOpt = ~check_option(varargin,{'p2c','c2c'});

if nargin > 1 && isnumeric(varargin{1})
  id = varargin{1};
else
  id = job.grains.id;
end

% the original child orientation
ori = job.grainsPrior('id',id).meanOrientation;

% all variants
oriV = variants(job.p2c, ori);
numV = size(oriV,2);

% fit with neighboring grains
A = job.grains('id',id).neighbors('matrix','maxId',max(job.grains.id));
A = A(id,:);
[grainInd,nId] = find(A);
grainInd = grainInd(:); nId = nId(:);

% some neighbors coorespond to parent and some to child grains
nInd = job.grains.id2ind(nId);
isParent = job.grains.phaseId(nInd) == job.parentPhaseId;
oriP = job.grains(nInd(isParent)).meanOrientation;

isChild = job.grains.phaseId(nInd) == job.childPhaseId;
oriChildV = variants(job.p2c,job.grains(nInd(isChild)).meanOrientation);

% compute fits to all neighbors
fit = nan(length(nId),numV);
for iV = 1:numV

  % parent - parent fit
  fit(isParent,iV) = angle(oriV(grainInd(isParent),iV), oriP);

  % parent - child fit
  %fit(isChild,iV) = min(angle(oriV(grainInd(isChild),iV), oriChildV),[],2);

end

% accumulate votes, i.e. compute a probability for each grain / parentId
% combination
votes = accumVotes(repmat(grainInd,1,numV), repmat(1:numV,length(grainInd),1), fit,...
  max(grainInd), varargin{:});

% output
if nargin > 1 && isnumeric(varargin{1})
  out = votes;
else
  job.votes = votes;
  out = job;
end

end


text:
function [X,obj,err,iter] = groupl1(A,B,G,opts)

% Solve the group l1-minimization problem by ADMM
%
% min_X \sum_{i=1}^n\sum_{g in G} ||(x_i)_g||_2, s.t. AX=B
%
% x_i is the i-th column of X
% ---------------------------------------------
% Input:
%       A       -    d*na matrix
%       B       -    d*nb matrix
%       G       -    a cell indicates a partition of 1:na
%       opts    -    Structure value in Matlab. The fields are
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       X       -    na*nb matrix
%       obj     -    objective function value
%       err     -    residual ||AX-B||_F
%       iter    -    number of iterations
%
% version 1.0 - 18/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

[d,na] = size(A);
[~,nb] = size(B);

X = zeros(na,nb);
Z = X;
Y1 = zeros(d,nb);
Y2 = X;

AtB = A'*B;
I = eye(na);
invAtAI = (A'*A+I)\I;

iter = 0;
for iter = 1 : max_iter
    Xk = X;
    Zk = Z;
    % update X
    for i = 1 : nb
        X(:,i) = prox_gl1(Z(:,i)-Y2(:,i)/mu,G,1/mu);
    end
    % update Z
    Z = invAtAI*(-(A'*Y1-Y2)/mu+AtB+X);
    dY1 = A*Z-B;
    dY2 = X-Z;
    chgX = max(max(abs(Xk-X)));
    chgZ = max(max(abs(Zk-Z)));
    chg = max([chgX chgZ max(abs(dY1(:))) max(abs(dY2(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = compute_obj(X,G);
            err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y1 = Y1 + mu*dY1;
    Y2 = Y2 + mu*dY2;
    mu = min(rho*mu,max_mu);
end
obj = compute_obj(X,G);
err = sqrt(norm(dY1,'fro')^2+norm(dY2,'fro')^2);

function obj = compute_obj(X,G)
obj = 0;
for i = 1 : size(X,2)
    x = X(:,i);
    for j = 1 : length(G)
        obj = obj + norm(x(G{j}));
    end
end

author: canyilu

repo: LibADMM-toolbox


text:
%GETRECTSUBPIX  Retrieves a pixel rectangle from an image with sub-pixel accuracy
%
%     dst = cv.getRectSubPix(src, patchSize, center)
%     dst = cv.getRectSubPix(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Source image, 8-bit integer or 32-bit floating-point, 1- or
%   3-channels.
% * __patchSize__ Size of the extracted patch `[w,h]`.
% * __center__ Floating-point coordinates of the center of the extracted
%   rectangle within the source image. The center `[x,y]` must be inside the
%   image.
%
% ## Output
% * __dst__ Extracted patch that has the size `PatchSize`, the same number of
%   channels as `src`, and the specified type in `PatchType`.
%
% ## Options
% * __PatchType__ Depth of the extracted pixels. By default (-1), they have
%   the same depth as `src`. Supports either `uint8` or `single`.
%
% The function cv.getRectSubPix extracts pixels from `src`:
%
%     dst(x,y) = src(x + center(1) - (size(dst,2)-1)*0.5,
%                    y + center(2) - (size(dst,1)-1)*0.5)
%
% where the values of the pixels at non-integer coordinates are retrieved
% using bilinear interpolation. Every channel of multi-channel images is
% processed independently. Also the image should be a single channel or three
% channel image. While the center of the rectangle must be inside the image,
% parts of the rectangle may be outside. In this case, the pixel values
% outside of the image are extrapolated.
%
% See also: cv.Rect.crop, cv.warpAffine, cv.warpPerspective
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [I,EI] = information(P)
%
%[I,EI] = information(P)
% Obtains the information (pointwise) and its expectation from a distribution of 1
% or two variables, hence:
% - for Pxy a bivariate distribution:
%  [I_Pxy,H_Pxy] = information(Pxy)
% obtains I_Pxy = log Pxy and H_Pxy is joint entropy
%
% - for Px a (row,column) univariate distribution
%  [I_Px,H_Px] = information(Px)
% obtains the self-information and the entropy.

error(nargchk(1,1,nargin));

I = -log2(P);
if nargout == 2
    lp = logical(P);
    if any(size(P) == 1)%row distribution or column distribution, resp
        EI = sum(P(lp) .* I(lp));
    else%bivariate, be very careful
        EI = sum(sum(P(lp) .* I(lp)));
    end
end
return%I,EI

author: Sable

repo: mcbench-benchmarks


text:
% Computes the speed difference between Matlab's binornd and
% Lightspeed's randbinom.

p = 0.13;
n = 123;
tim = [];

nsamples = 1e4;
y = zeros(nsamples,1);
tic
for i = 1:nsamples
  y(i) = randbinom(p,n);
end
tim(1) = toc;
%g = int_hist(y+1,n+1)/nsamples;
tic
for i = 1:nsamples
  y(i) = binornd(n,p);
end
tim(2) = toc;
fprintf('Time for binornd: %g\n', tim(2));
fprintf('Time for randbinom: %g (%g times faster)\n', tim(1), tim(2)/tim(1));

if 0
  % test validity of the sampler (use nsamples = 1e5)
  x = 0:n;
  f = binopdf(x,n,p);
  plot(x,f,x,g)
  legend('true','estimated')
end


text:
function c=xcorr(c,varargin)

% C = XCORR(C)
% This function calculates and fills in the correlation and lag fields in a
% correlation object. The input is a correlation object, presumeably with
% empty correlation and lag fields. c.C is a matrix of maximum correlation
% values normalized so that autocorrelations are 1. c.L is the lag time in
% between the two waveforms required for maximum correlation. To acheive
% maximum alignment, the value in position (i,j) should be added to the
% trigger time of trace j, or subtracted from trigger i. Traces can be
% aligned with the routine ADJUSTTRIG.
%
% By default, peak cross correlation values and lag times are NOT
% interpolated for sub-sample lag time because this requires a 30-40%
% increase in CPU time. For some uses, such as relative earthquake
% locations and coda wave interferometry, such precision is necessary. In
% these cases, consider the INTERP option below.
%
% C = XCORR(C,[PRETRIG POSTTRIG]);
% Perform cross correlation on a cropped portion of the data only. This is
% useful when you wish to keep the entire waveform but align the traces based on
% the correlation of a particular wavelet. PRETRIG and POSTTRIG are the
% time in seconds relative to the trigger time. Note that PRETRIG is
% negative for times before the trigger.
%
% C = XCORR(C,...,'1xr') Use single trace against one row algorithm
% (default).
%
% C = XCORR(C,...,'interp') Perform 2nd-order polynomial fitting to
% estimate sub-sample lag time. Sub-sample alignment requires an addition
% 30-40% CPU time but results in highest precision lag times possible.
%
% C = XCORR(C,...,'row',INDEX) Run correlation only on the traces specified
% by INDEX. Each trace of INDEX is correlated against the entire set of
% waveforms. This is useful if a small number of traces has been added to a
% large correlation matrix. Insead of recomputing the entire correlation
% and lag matrices, the routine allows only the "added" lines to be filled
% in. The syntax is a bit clunky. This routine requires the 'row' algorithm
% and the INDEX list. The two must be used together. Polynomial
% interpolation of lag values is always used with this algorithm.
%
%
% -- DEPRICATED ALGORITHMS ---------------------------------------------
% Because these algorithms seem to have little or no advantages over the
% '1xr' algorithm they will likely not be updated or improved.
%
% C = XCORR(C,...,'dec') Same as 1xr but decomposes complex numbers for
% calculation. Mathworks suggests this approach may be faster in some
% circumstances. However, initial testing found it slower that the 1xr
% algorithm.
%
% C = XCORR(C,...,'1x1') Use single trace against single trace algorithm.
% Conceivably faster when memory is very limited. In practice I have yet to
% encounter a situation where this algorithm benchmarks faster than the
% 1xr.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% GET INPUT PARAMETERS
if ~isa(c,'correlation')
    disp('First input parameter must be a correlation object');
end
algorithm = '1xr';
c1 = correlation;
c1 = set(c1,'WAVEFORM', get(c,'WAVEFORM') );
c1 = set(c1,'TRIG', get(c,'TRIG') );

% CHECK FOR TRACE SUBSET
if length(varargin)>1
    if isa(varargin{end},'double')
        index = varargin{end};
        varargin = varargin(1:end-1);
    end;
end;

    % CHECK ALGORITHM
if ~isempty(varargin)
    if ischar(varargin{end})
        algorithm = lower(varargin{end});
        varargin = varargin(1:end-1);
    end;
end;

    % APPLY CROPPING
if ~isempty(varargin)
    if length(varargin{end})==2       % check for cropping values
        pretrig =  varargin{1}(1);
        posttrig = varargin{1}(2);
        c1 = crop(c,pretrig,posttrig);
    end
end;

% CREATE MATRIX OF DATA FROM WAVEFORM ARRAY
% The correlation object is modified here into a Matlab structure that
% shares similar fields to the object except that the the trace data is
% stored in a matrix instead of in a waveform object. This structure is
% passed to the correlation subroutines as 'd'. The matrix structure
% improves computation speed because fft and other routines are optomized
% for matrices. This structure is based on the correlation object version
% 0.
d.start = get(c1.W,'START_MATLAB');
d.Fs    = get(c1.W(1),'Fs');
d.trig  = c1.trig;
d.w = double(c1.W);
% d.w     = [];
% for i = 1:length(c1.W)
%     d.w(:,i) = get(c1.W(i),'DATA');
% end;
clear c1

% EXECUTE CROSS CORRELATION
if exist('pretrig','var') && exist('posttrig','var')
    disp(['using ' algorithm ' algorithm on the time interval [' num2str(pretrig) ' ' num2str(posttrig) '] ...' ]);
else
    %disp(['using ' algorithm ' algorithm ...']);
end
if strcmp(algorithm,'1x1')==1
	d = xcorr1x1(d);
  elseif strcmp(algorithm,'1xr')
  	d = xcorr1xr(d,0);
  elseif strncmpi(algorithm,'int',3)
  	d = xcorr1xr(d,1);
  elseif strcmp(algorithm,'dec')
  	d = xcorrdec(d);
  elseif strcmp(algorithm,'row')
  	d = xcorrrow(d,c,index);
else
    error('Correlation algorithm not recognized');
end;

% ASSIGN CORRELATION PARAMETERS TO ORIGINAL DATA
c = set( c , 'CORR' , d.C );
c = set( c , 'LAG'  , d.L  );
clear d

author: geoscience-community-codes

repo: GISMO


text:
clear all, close all, clc

% Define domain
dx = 0.001;
L = pi;
x = (-1+dx:dx:1)*L;
n = length(x);   nquart = floor(n/4);

% Define hat function
f = 0*x;
f(nquart:2*nquart) = 4*(1:nquart+1)/n;
f(2*nquart+1:3*nquart) = 1-4*(0:nquart-1)/n;
plot(x,f,'-k','LineWidth',1.5), hold on

% Compute Fourier series
CC = jet(20);
A0 = sum(f.*ones(size(x)))*dx;
fFS = A0/2;
for k=1:20
    A(k) = sum(f.*cos(pi*k*x/L))*dx; % Inner product
    B(k) = sum(f.*sin(pi*k*x/L))*dx;
    fFS = fFS + A(k)*cos(k*pi*x/L) + B(k)*sin(k*pi*x/L);
    plot(x,fFS,'-','Color',CC(k,:),'LineWidth',1.2)
end

%% Plot amplitudes
figure
clear ERR
clear A
fFS = A0/2;
A(1) = A0/2;
ERR(1) = norm(f-fFS);
kmax = 100;
for k=1:kmax
    A(k+1) = sum(f.*cos(pi*k*x/L))*dx;
    B(k+1) = sum(f.*sin(pi*k*x/L))*dx;
%     plot(x,B(k)*sin(2*k*pi*x/L),'k-','LineWidth',1.2);
    fFS = fFS + A(k+1)*cos(k*pi*x/L) + B(k+1)*sin(k*pi*x/L);
    ERR(k+1) = norm(f-fFS)/norm(f);
end
thresh = median(ERR)*sqrt(kmax)*4/sqrt(3);
r = max(find(ERR>thresh));
r = 7;
subplot(2,1,1)
semilogy(0:1:kmax,A,'k','LineWidth',1.5)
hold on
semilogy(r,A(r+1),'bo','LineWidth',1.5)
xlim([0 kmax])
ylim([10^(-7) 1])
subplot(2,1,2)
semilogy(0:1:kmax,ERR,'k','LineWidth',1.5)
hold on
semilogy(r,ERR(r+1),'bo','LineWidth',1.5)

author: dynamicslab


text:
function [data3d] = augmentPascal3Ddata(data,class)
%AUGMENTPASCAL3DDATA Summary of this function goes here
%   Detailed explanation goes here

N = length(data);
globals;
load(fullfile(datadir,'subtypeClusters'));
data3d = data;
cInd = pascalClassIndex(class);
notFound = 0;
for i=1:N
    bbox = data(i).bbox;
    if(data(i).flip)
        bbox(1) = data(i).imsize(2) + 1 - bbox(1) - bbox(3);
    end
    bbox(3:4) = bbox(1:2)+bbox(3:4)-1;
    pascal3Dfile = fullfile(PASCAL3Ddir,'Annotations',[class '_pascal'],[data(i).voc_image_id '.mat']);
    if(exist(pascal3Dfile,'file'))
      record = load(pascal3Dfile);record = record.record;
      bboxes = vertcat(record.objects(:).bbox);
      bboxes(~strcmp(class,{record.objects(:).class}),:) = 0;
      objectInd = (bboxOverlap(bbox,bboxes,0.9));
      if(objectInd)
          viewpoint = record.objects(objectInd).viewpoint;
          [rot,euler] = viewpointToRots(viewpoint);
          data3d(i).rotP3d=rot;
          data3d(i).euler=euler;
          data3d(i).subtype = record.objects(objectInd).cad_index;
          data3d(i).subtype = subtypeClusters{cInd}(data3d(i).subtype);
      end
      data3d(i).objectIndP3d = objectInd;
      if(data3d(i).flip)
          data3d(i).euler = data3d(i).euler.*[-1;1;-1];
          rotX = diag([1,-1,-1]);
          data3d(i).rotP3d = rotX*angle2dcm(data3d(i).euler(1), data3d(i).euler(2)-pi/2, -data3d(i).euler(3),'ZXZ')';
      end
      %disp('blah')
    else
        notFound = notFound+1;
        data3d(i).rotP3d = nan(3,3);
        data3d(i).euler = nan(3,1);
        data3d(i).subtype = nan;
        data3d(i).objectIndP3d = nan;
    end

end

if(notFound~=0)
    fprintf('PASCAL 3D annotations not found - %d\n',notFound);
end
end

function [R,euler] = viewpointToRots(vp)
        euler = [vp.azimuth vp.elevation vp.theta]' .* pi/180;
        rotX = diag([1,-1,-1]);
        R1 = angle2dcm(euler(3), euler(2)-pi/2, -euler(1),'ZXZ'); %took a lot of work to figure this formula out !!
        euler = euler([3 2 1]);
        R = rotX*R1';
end

function isSame = bboxOverlap(bbox1,bboxes,thresh)
    bboxes(:,3:4) = [bboxes(:,3)-bboxes(:,1) bboxes(:,4)-bboxes(:,2)];
    bbox1(3:4) = [bbox1(:,3)-bbox1(:,1) bbox1(:,4)-bbox1(:,2)];
    intersectionArea=rectint(bbox1,bboxes);
    x_g = bbox1(1); y_g = bbox1(2);
    x_p = bboxes(:,1); y_p = bboxes(:,2);
    width_g = bbox1(3); height_g = bbox1(4);
    width_p = bboxes(:,3); height_p = bboxes(:,4);
    unionCoords=[min(x_g,x_p),min(y_g,y_p),max(x_g+width_g-1,x_p+width_p-1),max(y_g+height_g-1,y_p+height_p-1)];
    unionArea=(unionCoords(:,3)-unionCoords(:,1)+1).*(unionCoords(:,4)-unionCoords(:,2)+1);
    overlapArea=intersectionArea(:)./unionArea(:);
    [~,isSame] = max(overlapArea);
end

repo: CategoryShapes


text:
function classifier = oclapsvmp(options,data)
% {oclapsvmp} trains a One-Class Laplacian SVM classifier on the primal.
%
%      classifier = oclapsvmp(options,data)
%
%      Note: This function is only a wrapper to 'lapsvmp', where the
%            contribute of unlabeled data is discarded. See 'lapsvmp'
%            for the description of all required parameters. All parameters
%            related to unlabeled data are obviously not checked.
%
% Author: Stefano Melacci (2012), Salvatore Frandina (2012)
%         mela@dii.unisi.it, salvatore.frandina@gmail.com

options.Hinge=1;
options.UseBias=1;
data.Y(data.Y==-1)=0;
classifier=lapsvmp(options,data);
classifier.name='oclapsvmp';

author: viggin

repo: domain-adaptation-toolbox


text:
function javaudd=dftoolgetudd(uddcmd,varargin);

%   $Revision: 1.1.6.2 $  $Date: 2004/01/24 09:35:57 $
%   Copyright 2003-2004 The MathWorks, Inc.

% unwrap any UDD objects
for i=1:length(varargin)
   if isa(varargin{i}, 'com.mathworks.jmi.bean.UDDObject')
      varargin{i}=handle(varargin{i});
   end
end

% wrap the return UDD object
if nargin == 1
   javaudd=java(eval(uddcmd));
else
   javaudd=java(feval(uddcmd,varargin{:}));
end

author: zouchuhang

repo: LayoutNet


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%

for i = 1: p
  X_poly(:, i) = X' .^i;
end

% =========================================================================

end

author: benoitvallon

repo: coursera-machine-learning


text:
function montage = spm_eeg_montage_ui(montage)
% GUI for EEG montage (rereference EEG data to new reference channel(s))
% FORMAT montage = spm_eeg_montage_ui(montage)
%
% montage     - structure with fields:
%   tra       - MxN matrix
%   labelnew  - Mx1 cell-array - new labels
%   labelorg  - Nx1 cell-array - original labels
%
% Output is empty if the GUI is closed.
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Jean Daunizeau
% $Id: spm_eeg_montage_ui.m 7755 2019-12-16 13:19:28Z spm $

% Create the figure
%--------------------------------------------------------------------------
fig  = figure;
S0   = spm('WinSize','0',1);
pos  = get(fig,'position');
pos2 = [40 70 pos(3)-60 pos(4)-100];
pos  = [S0(1) S0(2) 0 0] + [pos(1) pos(2) 1.8*pos(3) pos(4)];
set(fig,...
    'menubar',     'none',...
    'position',    pos,...
    'numbertitle', 'off',...
    'name',        'Montage edition');

addButtons(fig);

% Display the uitable component
%--------------------------------------------------------------------------
table    = cat(2,montage.labelnew(:),num2cell(montage.tra));
colnames = cat(2,'channel labels',montage.labelorg(:)');

pause(1e-1) % This is weird, but fixes java troubles.
ht       = my_uitable(table,colnames);
set(ht,'position',pos2, 'units','normalized');

% Display the matrix representation of the montage
%--------------------------------------------------------------------------
ax = axes('position',[0.6 0.18 0.4 0.77]);
hi = imagesc(montage.tra,'parent',ax);
axis('image');
colormap('bone');
zoom(fig,'on');

% Store info in figure's userdata and wait for user interaction
%--------------------------------------------------------------------------
ud.hi      = hi;
ud.ht      = ht;
ud.montage = montage;
set(fig,'userdata',ud);
uiwait(fig);

% Get the montage from the GUI
%--------------------------------------------------------------------------
try
    ud = get(fig,'userdata');
    montage = ud.montage;
    close(fig);
catch % GUI was closed without 'OK' button
    montage = [];
end

%==========================================================================
function doAddRow(obj,evd,h)
%==========================================================================
% 'add row' button subfunction
ud = get(h,'userdata');
[M,newLabels] = getM(ud.ht);
M = [M;zeros(1,size(M,2))];
newLabels = cat(1,newLabels(:),num2str(size(M,1)));
set(ud.ht,'units','pixels');
pos = get(ud.ht,'Position');
delete(ud.ht);
table = cat(2,newLabels,num2cell(M));
colnames = cat(2,'channel labels',ud.montage.labelorg(:)');
pause(1) % This is weird, but fixes java troubles.
ht = my_uitable(table,colnames);
set(ht,'position',pos, 'units','normalized');
ud.ht = ht;
set(h,'userdata',ud);
doCheck(obj,evd,h);

%==========================================================================
function doLoad(obj,evd,h)
%==========================================================================
% 'load' button subfunction
[t,sts] = spm_select(1,'mat','Load montage file');
if sts
    montage = load(t);
    if ismember('montage', fieldnames(montage))
        montage    = montage.montage;
        ud         = get(h,'userdata');
        set(ud.ht,'units','pixels');
        pos        = get(ud.ht,'Position');
        delete(ud.ht);
        table      = cat(2,montage.labelnew(:),num2cell(montage.tra));
        colnames   = cat(2,'channel labels',montage.labelorg(:)');
        pause(1) % This is weird, but fixes java troubles.
        ht         = my_uitable(table,colnames);
        set(ht,'position',pos,...
            'units','normalized');
        ud.ht      = ht;
        ud.montage = montage;
        set(h,'userdata',ud);
        pause(1)
        doCheck(obj,evd,h);
    else
        spm('alert!','File did not contain any montage!','Montage edition');
    end
end

%==========================================================================
function doSave(obj,evd,h)
%==========================================================================
% 'save as' button subfunction
doCheck(obj,evd,h);
ud = get(h,'userdata');
[M,newLabels] = getM(ud.ht);
% delete row if empty:
ind              = ~any(M,2);
M(ind,:)         = [];
newLabels(ind)   = [];
montage.tra      = M;
montage.labelorg = ud.montage.labelorg;
montage.labelnew = newLabels;
[filename, pathname] = uiputfile({'*.mat','MAT-files (*.mat)'}, ...
    'Save montage', 'SPMeeg_montage.mat');
if ~isequal(filename, 0)
    save(fullfile(pathname, filename), 'montage', spm_get_defaults('mat.format'));
end

%==========================================================================
function doOK(obj,evd,h)
%==========================================================================
% 'OK' button subfunction
doCheck(obj,evd,h);
ud               = get(h,'userdata');
[M, newLabels]   = getM(ud.ht);
% delete row if empty:
ind              = ~any(M,2);
M(ind,:)         = [];
newLabels(ind)   = [];
montage.tra      = M;
montage.labelorg = ud.montage.labelorg(:);
montage.labelnew = newLabels(:);
ud.montage = montage;
set(h,'userdata',ud);
uiresume(h);

%==========================================================================
function doCheck(obj,evd,h)
%==========================================================================
% Update the montage display
ud = get(h,'userdata');
M  = getM(ud.ht);
set(ud.hi,'cdata',M);
set(gca,'xlim',[0.5 size(M,1)]);
set(gca,'ylim',[0.5 size(M,2)]);
axis('image');
drawnow;

%==========================================================================
function [M,newLabels] = getM(ht)
%==========================================================================
% extracting montage from java object
nnew = get(ht,'NumRows');
nold = get(ht,'NumColumns')-1;
M    = zeros(nnew,nold);
data = get(ht,'data');
for i =1:nnew
    if ~isempty(data(i,1))
        newLabels{i} = data(i,1);
    else
        newLabels{i} = [];
    end
    for j =1:nold
        if ~isempty(data(i,j+1))
            if ~ischar(data(i,j+1))
                M(i,j) = data(i,j+1);
            else
                M0 = str2double(data(i,j+1));
                if ~isempty(M0)
                    M(i,j) = M0;
                else
                    M(i,j) = 0;
                end
            end
        else
            M(i,j) = 0;
        end
    end
end

%==========================================================================
function addButtons(h)
%==========================================================================
% adding buttons to the montage GUI
hAdd = uicontrol('style','pushbutton',...
    'string','Add row','callback',{@doAddRow,h},...
    'position',[60 20 80 20]);
set(hAdd,'units','normalized');
hLoad = uicontrol('style','pushbutton',...
    'string','Load file','callback',{@doLoad,h},...
    'position',[180 20 80 20]);
set(hLoad,'units','normalized');
hSave = uicontrol('style','pushbutton',...
    'string','Save as','callback',{@doSave,h},...
    'position',[280 20 80 20]);
set(hSave,'units','normalized');
hOK = uicontrol('style','pushbutton',...
    'string',' OK ','callback',{@doOK,h},...
    'position',[400 20 80 20]);
set(hOK,'units','normalized');
hCheck = uicontrol('style','pushbutton',...
    'string',' Refresh display ','callback',{@doCheck,h},...
    'position',[760 20 120 20]);
set(hCheck,'units','normalized');

%==========================================================================
function [ht,hc] = my_uitable(varargin)
%==========================================================================
% conversion layer for various MATLAB versions
persistent runOnce
try
    if spm_check_version('matlab','8.4') >= 0
        if isempty(runOnce)
            warning('Consider migrating to the new uitable component.');
            runOnce = true;
        end
        [ht,hc] = uitable('v0',varargin{:});
    else
        [ht,hc] = spm_uitable(varargin{:});
    end
catch
    [ht,hc]     = deal([]);
end

author: spm


text:
function exclude = getLvdGAExcludeList()
    exclude = {};
    exclude{end+1} = 'Distance to Ref. Spacecraft';
    exclude{end+1} = 'Distance to Ref. Station';
    exclude{end+1} = 'Relative Vel. to Ref. Spacecraft';
    exclude{end+1} = 'Relative Pos. of Ref. Spacecraft (In-Track)';
    exclude{end+1} = 'Relative Pos. of Ref. Spacecraft (Cross-Track)';
    exclude{end+1} = 'Relative Pos. of Ref. Spacecraft (Radial)';
    exclude{end+1} = 'Relative Pos. of Ref. Spacecraft (In-Track; Ref. SC-centered)';
    exclude{end+1} = 'Relative Pos. of Ref. Spacecraft (Cross-Track; Ref. SC-centered)';
    exclude{end+1} = 'Relative Pos. of Ref. Spacecraft (Radial; Ref. SC-centered)';
    exclude{end+1} = 'Relative SMA of Ref. Spacecraft';
    exclude{end+1} = 'Relative Eccentricity of Ref. Spacecraft';
    exclude{end+1} = 'Relative Inclination of Ref. Spacecraft';
    exclude{end+1} = 'Relative RAAN of Ref. Spacecraft';
    exclude{end+1} = 'Relative Argument of Periapsis of Ref. Spacecraft';
    exclude{end+1} = 'Elevation Angle w.r.t. Ref. Station';
    exclude{end+1} = 'Line of Sight to Ref. Spacecraft';
    exclude{end+1} = 'Line of Sight to Ref. Station';
    exclude{end+1} = 'Sun-centric Position (X)';
    exclude{end+1} = 'Sun-centric Position (Y)';
    exclude{end+1} = 'Sun-centric Position (Z)';
    exclude{end+1} = 'Body-centric Position (X)';
    exclude{end+1} = 'Body-centric Position (Y)';
    exclude{end+1} = 'Body-centric Position (Z)';
    exclude{end+1} = 'Body-centric Velocity (X)';
    exclude{end+1} = 'Body-centric Velocity (Y)';
    exclude{end+1} = 'Body-centric Velocity (Z)';
    exclude{end+1} = 'Body-Fixed Velocity (X)';
    exclude{end+1} = 'Body-Fixed Velocity (Y)';
    exclude{end+1} = 'Body-Fixed Velocity (Z)';
    exclude{end+1} = 'Body-Fixed Velocity';
    exclude{end+1} = 'Semi-major Axis';
    exclude{end+1} = 'Eccentricity';
    exclude{end+1} = 'Inclination';
    exclude{end+1} = 'Right Asc. of the Asc. Node';
    exclude{end+1} = 'Argument of Periapsis';
    exclude{end+1} = 'True Anomaly';
    exclude{end+1} = 'Mean Anomaly';
    exclude{end+1} = 'Orbital Period';
    exclude{end+1} = 'Radius of Periapsis';
    exclude{end+1} = 'Radius of Apoapsis';
    exclude{end+1} = 'Radius of Spacecraft';
    exclude{end+1} = 'Altitude of Apoapsis';
    exclude{end+1} = 'Altitude of Periapsis';
    exclude{end+1} = 'Speed of Spacecraft';
    exclude{end+1} = 'Flight Path Angle';
    exclude{end+1} = 'Equinoctial H1';
    exclude{end+1} = 'Equinoctial K1';
    exclude{end+1} = 'Equinoctial H2';
    exclude{end+1} = 'Equinoctial K2';
    exclude{end+1} = 'Altitude';
    exclude{end+1} = 'Longitude (East)';
    exclude{end+1} = 'Latitude (North)';
    exclude{end+1} = 'Longitudinal Drift Rate';
    exclude{end+1} = 'Surface Velocity';
    exclude{end+1} = 'Vertical Velocity';
    exclude{end+1} = 'Central Body ID';
    exclude{end+1} = 'C3 Energy';
    exclude{end+1} = 'Hyperbolic Velocity Unit Vector X';
    exclude{end+1} = 'Hyperbolic Velocity Unit Vector Y';
    exclude{end+1} = 'Hyperbolic Velocity Unit Vector Z';
    exclude{end+1} = 'Hyperbolic Velocity Vector Right Ascension';
    exclude{end+1} = 'Hyperbolic Velocity Vector Declination';
    exclude{end+1} = 'Hyperbolic Velocity Magnitude';
    exclude{end+1} = 'Distance to Ref. Celestial Body';
    exclude{end+1} = 'Elevation Angle of Ref. Celestial Body';
    exclude{end+1} = 'Liquid Fuel/Ox Mass';
    exclude{end+1} = 'Monopropellant Mass';
    exclude{end+1} = 'Xenon Mass';
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function I = imreadHistEq( rgbFile, lwirFile )
RGB = imread( rgbFile );
T = rgb2gray( imread( lwirFile ) );
I = cat(3, RGB, T);
% I(:,:,4) = histeq(I(:,:,4));
I(:,:,4) = I(:,:,4) .* (255.0 / 141);
end

author: SoonminHwang

repo: rgbt-ped-detection


text:
function c = mycross(a,b,dim)
%% MYCROSS cross product
%
% c = mycross(a,b) computes the cross product of vectors a and b. It is
% used to replace the built-in cross function which is two times slower.
%
% c = mycross(a,b,dim) the input a, b could be matrix of vectors. The input
% dim specifies the dimension. By default, dim = 2. Namely each row of a
% and b is a vector in R3.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('dim','var'), dim = 2; end

% Check dimensions
if (size(a,dim)~=3) || (size(b,dim)~=3),
  error(message('MATLAB:cross:InvalidDimAorBForCrossProd'))
end

c = zeros(size(a,1),3);
c(:,1) = a(:,2).*b(:,3) - a(:,3).*b(:,2);
c(:,2) = a(:,3).*b(:,1) - a(:,1).*b(:,3);
c(:,3) = a(:,1).*b(:,2) - a(:,2).*b(:,1);

repo: Aero-matlab


text:
%%%% discretization of direct space, to be used for plotting the field
%%%% distribution; X,Y = vectors to be used for FFT; Xi,Yi = vectors to be
%%%% used for interpolation
function [X,Y,Xi,Yi]=prcellgrid(a1,a2,N1,N2)
X=[]; Y=[]; %spatial coordinates after FFT
for l=1:N1
    for m=1:N2
        X(l,m)=(l-1-(N1-1)/2)*a1(1)/N1 + (m-1-(N2-1)/2)*a2(1)/N2;
        Y(l,m)=(l-1-(N1-1)/2)*a1(2)/N1 + (m-1-(N2-1)/2)*a2(2)/N2;
    end
end

M=501;
Xi=[]; Yi=[]; %finer discretization for field interpolation
for l=1:M
    for m=1:M
        Xi(l,m)=(l-1-(M-1)/2)*a1(1)/M + (m-1-(M-1)/2)*a2(1)/M;
        Yi(l,m)=(l-1-(M-1)/2)*a1(2)/M + (m-1-(M-1)/2)*a2(2)/M;
    end
end


text:
function planC = faceOff(planC, scanNum, tmpDirPath, exportDICOMFlag, outputDICOMPath)

%% parse input and set up parameters, variables
if ~exist('exportDICOMFlag', 'var') && ~exist('outputDICOMPath','var')
        exportDICOMFlag = 0;
end

if exportDICOMFlag && (~exist('outputDICOMPath','var') || isempty(outputDICOMPath))
    outputDICOMPath = pwd;
end

if ~exist('tmpDirPath', 'var') || isempty(tmpDirPath)
    tmpDirPath = fullfile(getCERRPath,'ImageRegistration','tmpFiles');
end

%% ANTs setup
% optS = getCERROptions;
% if ~exist(optS.antspath_dir,'dir')
%     error(['ANTSPATH ' optS.antspath_dir ' not found on filesystem. Please review CERROptions.']);
% end
% antspath = fullfile(optS.antspath_dir,'bin');
% setenv('ANTSPATH', antspath);
% antsScriptPath = fullfile(optS.antspath_dir, 'Scripts');
% antsCERRScriptPath = fullfile(getCERRPath,'CERR_core', 'ImageRegistration', 'antsScripts');
% if isunix
%     setenv('PATH',[antspath ':' antsScriptPath ':' antsCERRScriptPath ':' getenv('PATH')]);
% else
%     setenv('PATH',[antspath ';' antsScriptPath ';' antsCERRScriptPath ';' getenv('PATH')]);
% end

%% register scan to template

% set affine transform
inputCmdFile = fullfile(getCERRPath,'ImageRegistration','antsScripts','SyN_a.txt');
[templateImageFile, faceMaskImageFile] = loadTemplatePlanC;
[ ~, planC, ~] = register_scans(templatePlanC, 1, planC, scanNum, 'QuickSyn ANTs', 'ANTs', '', [], [], [], inputCmdFile);

%% inverse-warp the template face mask to scan space
indexS = planC{end};
deformS = planC{indexS.deform}(end);
inverseFlag = 1;
warp_structures(deformS, strCreationScanNum, movStructNumsV, planC, faceMaskImageFIle, inverseFlag);

%% export CERR to DICOM
disp('Exporting to DICOM format...');
if exportDicomFlag
    export_planC_to_DICOM(planC, outputDICOMPath, 0);
%     exportCERRtoDICOM(origCerrPath,allLabelNamesC,outputCERRPath,...
%         outputDicomPath,algorithm,savePlancFlag)
end

function [templateImageFile, faceMaskImageFile] = loadTemplatePlanC
    templateImageFile = fullfile(getCERRPath,'Extras','faceOffBeta','template0.nii.gz');
    faceMaskImageFile = fullfile(getCERRPath,'Extras','faceOffBeta','deface_template_40_M_GD_8.nii.gz');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% FORMATFIGURE: Formats Matlab figures
%
% FormatFigure( FIG, FONTNAME, FONTSIZE, FIGSIZE, BORDER )
%
% Formats a figure for presentations or papers.  The font
% style, font size, figure size and border size can all be
% adjusted.  All parameters are optional.
%
% FIG     : Figure handle.
% FONTNAME: Name of the font as a string.
% FONTSIZE: Size of the font (points) as a vector.
%           Format: [axes_labels tick_labels]
% FIGSIZE : Size of the figure (inches) as a vector.
%           Format: [width height]
% BORDER  : Space around the figure (inches) as a vector.
%           Format: [left bottom right top] or
%                   [left/right  botom/top]
%
% Default settings of the inputs are:
%   FIG     : gcf (current figure)
%   FONTNAME: 'Helvetica'
%   FONTSIZE: [18 14]
%   FIGSIZE : [6 6]
%   BORDER  : [0.8 0.5]
%
% Note: Function has a problem settiing the font size of the
% legend sometimes.  If the legend font is the wrong size,
% call this function again.  The second call corrects the
% problem.  I believe this to be a Matlab problem.
%
% Written By  : Jeffrey DiCarlo
% Last updated: 07-26-00

function FormatFigure( fig, fontname, fontsize, figsize, border )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Format the input parameters.

% Load default paramemters.

if (~exist('fig','var'))
   fig = 0;
end
if (~exist('fontname'))
   fontname = 'Helvetica';
end
if (~exist('fontsize'))
   fontsize = [18 14];
end
if (~exist('figsize'))
   figsize  = [6 6];
end
if (~exist('border'))
   border = [0.8 0.5];
end

% Check the fontsize.

if (length(fontsize) == 1)
   fontsize = [fontsize fontsize];
end

% Check the figure size.

if (length(figsize) == 1)
   figsize = [figsize figsize];
end

% Check the border.

if (length(border) == 2)
   border = [border(1) border(1) border(2) border(2)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the figure axes.

if (exist('fig') ~= 1 | fig == 0)
   fig = gcf;
end
axs = get( fig, 'CurrentAxes' );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the current figure position.

if (strcmp(get(fig, 'Units'),'inches') == 1)
   pos = get(fig, 'Position');
else
   pos = [1 1 1 1];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the figure properties.

set( fig,         'Units', 'inches' );
set( fig,      'Position', [pos(1:2) figsize] );
set( fig, 'PaperPosition', [4.25-figsize(1)/2 5.5-figsize(2)/2 figsize] );
set( fig,         'Color', [1 1 1] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the axes properties.

set( get(axs,  'Title'), 'FontName', fontname, ...
   'FontSize', fontsize(1) );
set( get(axs, 'XLabel'), 'FontName', fontname, ...
	'FontSize', fontsize(1) );
set( get(axs, 'YLabel'), 'FontName', fontname, ...
   'FontSize', fontsize(1) );
set( get(axs, 'ZLabel'), 'FontName', fontname, ...
   'FontSize', fontsize(1) );
set( axs,                'FontName', fontname, ...
   'FontSize', fontsize(2) );

set( axs,    'Units', 'inches' );
set( axs, 'Position', [border(1:2) (figsize-border(1:2)-border(3:4))] );

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function result  = mp(Dict, x, options)
% Matching Pursuit
    if nargin < 3
        options = struct;
    end
    if isa(Dict, 'spx.dict.Operator')
        dict = Dict;
    elseif ismatrix(Dict)
        dict = spx.dict.MatrixOperator(Dict);
    else
        error('Unsupported operator.');
    end
    [m, n] = size(dict);

    % halting criteria
    stop_on_r_norm = false;
    stop_on_iterations = false;
    max_iterations = n;
    max_residual_norm = 0;
    if isfield(options, 'max_iterations')
        stop_on_iterations = true;
        max_iterations = options.max_iterations;
    end
    if isfield(options, 'max_residual_norm')
        stop_on_r_norm = true;
        max_residual_norm = options.max_residual_norm;
    end

    Gamma = [];
    r = x;
    z = zeros(n, 1);
    iter = 0;
    % norm of signal being reconstructed
    x_norm = norm(x);
    while true
        h = dict.apply_ctranspose(r);
        [~, index] = max(abs(h));
        coeff = h(index);
        z(index) = z(index) + coeff;
        r = r - coeff*dict.column(index);
        r_norm = norm(r);
        iter = iter + 1;
        if stop_on_iterations
            if iter >= max_iterations
                break;
            end
        end
        if stop_on_r_norm
            if r_norm < max_residual_norm
                break;
            end
        end
        if r_norm / x_norm < 1e-6
            % we have reached exact exact recovery
            break;
        end
        if iter > 4*n
            % we have reached highest number of iterations.
            break;
        end
    end
    result.z = z;
    result.r = r;
    result.iterations = iter;
    result.support = find(z);
end

author: indigits

repo: sparse-plex


text:
% ---------------------------------------------------------------------
% In this script trajectory optimization otherwise called experiment
% design for dynamic paramters identification is carried out.
%
% First, specify cost function (traj_cost_lgr) and constraints
% (traj_cnstr) for the optimziation. Then choose oprimization algorithm and
% specify trajectory parameters (duration, fundamental frequency, number of
% harmonics, initial (= final) positionin) and max/min positions,
% velocities and accelerations.
%
% Then script runs optimization, plots obtained trajectory and saves its
% parameters into a file.
% ---------------------------------------------------------------------
% get robot description
path_to_urdf = 'ur10e.urdf';
ur10 = parse_urdf(path_to_urdf);

% get mapping from full parameters to base parameters
include_motor_dynamics = 1;
[~, baseQR] = base_params_qr(include_motor_dynamics);

% Choose optimization algorithm: 'patternsearch', 'ga'
optmznAlgorithm = 'patternsearch';

% Trajectory parameters
traj_par.T = 25;          % period of signal
traj_par.wf = 2*pi/traj_par.T;    % fundamental frequency
traj_par.t_smp = 2e-1;   % sampling time
traj_par.t = 0:traj_par.t_smp:traj_par.T;  % time
traj_par.N = 7;          % number of harmonics
traj_par.q0 = deg2rad([0 -90 0 -90 0 0 ]');
% Use different limit for positions for safety
traj_par.q_min = -deg2rad([180  180  100   180  90   90]');
traj_par.q_max =  deg2rad([180  0    100   0    90   90]');
traj_par.qd_max = qd_max;
traj_par.q2d_max = [2 1 1 1 1 2.5]';

%  ----------------------------------------------------------------------
% Otimization
% -----------------------------------------------------------------------
A = []; b = [];
Aeq = []; beq = [];
lb = []; ub = [];

if strcmp(optmznAlgorithm, 'patternsearch')
    x0 = rand(6*2*traj_par.N,1);
    optns_pttrnSrch = optimoptions('patternsearch');
    optns_pttrnSrch.Display = 'iter';
    optns_pttrnSrch.StepTolerance = 1e-1;
    optns_pttrnSrch.FunctionTolerance = 10;
    optns_pttrnSrch.ConstraintTolerance = 1e-6;
    optns_pttrnSrch.MaxTime = inf;
    optns_pttrnSrch.MaxFunctionEvaluations = 1e+6;

    [x,fval] = patternsearch(@(x)traj_cost_lgr(x,traj_par,baseQR), x0, ...
                             A, b, Aeq, beq, lb, ub, ...
                             @(x)traj_cnstr(x,traj_par), optns_pttrnSrch);
elseif strcmp(optmznAlgorithm, 'ga')
    optns_ga = optimoptions('ga');
    optns_ga.Display = 'iter';
    optns_ga.PlotFcn = 'gaplotbestf'; % {'gaplotbestf', 'gaplotscores'}
    optns_ga.MaxGenerations = 50;
    optns_ga.PopulationSize = 1e+3; % in each generation.
    optns_ga.InitialPopulationRange = [-100; 100];
    optns_ga.SelectionFcn = 'selectionroulette';

    [x,fval] = ga(@(x)traj_cost_lgr(x,traj_par,baseQR), 6*2*traj_par.N,...
                  A, b, Aeq, beq, lb, ub, ...
                  @(x)traj_cnstr(x,traj_par), optns_ga);
else
    error('Chosen algorithm is not found among implemented ones');
end

% ------------------------------------------------------------------------
% Plotting obtained trajectory
% ------------------------------------------------------------------------
ab = reshape(x,[12,traj_par.N]);
a = ab(1:6,:); % sin coeffs
b = ab(7:12,:); % cos coeffs
c_pol = getPolCoeffs(traj_par.T, a, b, traj_par.wf, traj_par.N, traj_par.q0);
[q,qd,q2d] = mixed_traj(traj_par.t, c_pol, a, b, traj_par.wf, traj_par.N);

figure
subplot(3,1,1)
    plot(traj_par.t,q)
    ylabel('$q$','interpreter','latex')
    grid on
    legend('q1','q2','q3','q4','q5','q6')
subplot(3,1,2)
    plot(traj_par.t,qd)
    ylabel('$\dot{q}$','interpreter','latex')
    grid on
    legend('qd1','qd2','qd3','qd4','qd5','qd6')
subplot(3,1,3)
    plot(traj_par.t,q2d)
    ylabel('$\ddot{q}$','interpreter','latex')
    grid on
    legend('q2d1','q2d2','q2d3','q2d4','q2d5','q2d6')

% ----------------------------------------------------------------------
% Saving parameters of the optimized trajectory
% ----------------------------------------------------------------------
% %{
pathToFolder = 'trajectory_optmzn/optimal_trjctrs/';
t1 = strcat('N',num2str(traj_par.N),'T',num2str(traj_par.T));
if strcmp(optmznAlgorithm, 'patternsearch')
    filename = strcat(pathToFolder,'ptrnSrch_',t1,'QR.mat');
elseif strcmp(optmznAlgorithm, 'ga')
    filename = strcat(pathToFolder,'ga_',t1,'.mat');
elseif strcmp(optmznAlgorithm, 'fmincon')
    filename = strcat(pathToFolder,'fmncn_',t1,'.mat');
end
save(filename,'a','b','c_pol','traj_par')
%}

author: shamilmamedov


text:
%% t_meshCurvature
%
% Use tricurve_v01 to calculate curvature of a mesh
%
%

%% Load the class file.  A sphere is 1 and everything else 0
fName = fullfile(mrvDataRootPath,'anatomy','harmonic.nii.gz');
niClass = niftiRead(fName);
% showMontage(niClass.data)

%% Use Matlab to build the faces and edges
fv = isosurface(niClass.data,0.9);
smoothMode = 1; nIter = 5;
fv = smoothpatch(fv,smoothMode,nIter);
% This was suggested by Peyre.  It may work in the end.  But it takes time,
% as well, so I may as well stick with the tricurv_v01 for now.
%
% options.method = 'slow';
% fv.faces = perform_faces_reorientation(fv.vertices,fv.faces,options)';

%% One solution.  Slow, but this is the one used in meshFV2msh.m
% curvature1 = tricurv_v01(fv.faces,fv.vertices);

%% A second solution is much faster.  But it fails a lot.

% I will probably remove it from meshFV2msh.m  But there is an interesting
% exchange about it with Peyre.  See email and notes.
%
%   Umin is the direction of minimum curvature
%   Umax is the direction of maximum curvature
%   Cmin is the minimum curvature
%   Cmax is the maximum curvature
%   Cmean=(Cmin+Cmax)/2
%   Cgauss=Cmin*Cmax
%   Normal is the normal to the surface

% BUT, this routine fails a lot.  I wrote the authors.  I put a try/catch
% around it in meshFV2msh.m
% tic
% [Umin,Umax,Cmin,Cmax,Cmean,Cgauss,Normal] = ...
%     compute_curvature(fv.vertices,fv.faces);
% toc

%% Convert the fv mesh to a vistasoft mesh

msh = meshFV2msh(fv);
meshVisualize(msh);

%% End


text:
function [a,n,r,epe,eph,epw,usol_c,vsol_c] = mg_ns_q1cd(xy,ev,flowsol,domain,lncoarse,lnfine,outbnd)
%mg_ns_q1cd   convection-diffusion matrix generator for GMG (Navier-Stokes)
%   [a,n,r,epe,eph,epw,usol_c,vsol_c] = mg_ns_q1cd(xy,ev,flowsol,domain,lncoarse,lnfine,outbnd)
%   input
%           xy         vertex coordinate vector
%           ev         element mapping matrix
%           flowsol    current velocity field
%           domain     domain parameter, 1 for square, 3 for step
%           lncoarse   coarse grid index, log2(nc) for nc x nc square grid
%           lnfine     finest grid index, log2(nf) for nf x nf square grid
%           outbnd     location of outflow boundary
%   output
%           a          discrete diffusion operator
%           n          discrete convection operator
%           r          mass matrix
%           epe        viscosity normalised element peclet numbers
%           eph        flow specific element lengths
%           epw        centroid evaluated wind
%           usol_c     first convection coefficient on coarse grid nodes
%           vsol_c     second convection coefficient on coarse grid nodes
%
%   IFISS function: HCE; 18 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

% Modification of DJS routine mg_q1cd to allow convection coefficients to
% be passed via flowsol.  Analogous to femq1_cd.

% Square domain
if domain==1,
% Size of fine grid, obtained from underlying flow solution flowsol
   nfine  = length(flowsol)/2;
   nf = sqrt(nfine);
   x=xy(:,1); y=xy(:,2);
   nvtx=length(x);   nel=length(ev(:,1));

% Get flow values at coarse grid points
   lshift = lnfine - lncoarse;
   usol=flowsol(1:nfine);
   usol_grid = reshape(usol,nf,nf);
   usol_grid_c = usol_grid(1:2^lshift:nf,1:2^lshift:nf);
   usol_c = reshape(usol_grid_c,nvtx,1);
   vsol=flowsol(nfine+1:2*nfine);
   vsol_grid = reshape(vsol,nf,nf);
   vsol_grid_c = vsol_grid(1:2^lshift:nf,1:2^lshift:nf);
   vsol_c = reshape(vsol_grid_c,nvtx,1);

% Step domain
elseif domain==3,
% Various counts of fine/coarse points in small part (1) and large part (2) of step
% Logic borrowed from AR routine step_mg_prolong
   nelf = 2^lnfine;
   nfine = (((outbnd+1)/2)*nelf+1)*(nelf+1) - (nelf/2)^2;
   nelc = 2^lncoarse;
   ncoarse = (((outbnd+1)/2)*nelc+1)*(nelc+1) - (nelc/2)^2;
   nfh1 = (nelf/2);     nfv1 = (nelf/2+1);  nf1 = nfh1*nfv1;
   nfh2 = (outbnd*nelf/2+1); nfv2 = nelf+1; nf2 = nfh2*nfv2;
   nch1 = (nelc/2);     ncv1 = (nelc/2+1);  nc1 = nch1*ncv1;
   nch2 = (outbnd*nelc/2+1); ncv2 = nelc+1; nc2 = nch2*ncv2;

% Get flow values at coarse grid points
   lshift = lnfine-lncoarse;
   usol=flowsol(1:nfine);
   usol1_grid = reshape(usol(1:nf1),nfh1,nfv1);
   usol1_grid_c = usol1_grid(1:2^lshift:nfh1,1:2^lshift:nfv1);
   usol1_c = reshape(usol1_grid_c,nc1,1);
   usol2_grid = reshape(usol(nf1+1:nfine),nfh2,nfv2);
   usol2_grid_c = usol2_grid(1:2^lshift:nfh2,1:2^lshift:nfv2);
   usol2_c = reshape(usol2_grid_c,nc2,1);
   usol_c = [usol1_c;usol2_c];
   vsol=flowsol(nfine+1:2*nfine);
   vsol1_grid = reshape(vsol(1:nf1),nfh1,nfv1);
   vsol1_grid_c = vsol1_grid(1:2^lshift:nfh1,1:2^lshift:nfv1);
   vsol1_c = reshape(vsol1_grid_c,nc1,1);
   vsol2_grid = reshape(vsol(nf1+1:nfine),nfh2,nfv2);
   vsol2_grid_c = vsol2_grid(1:2^lshift:nfh2,1:2^lshift:nfv2);
   vsol2_c = reshape(vsol2_grid_c,nc2,1);
   vsol_c = [vsol1_c;vsol2_c];
end

%
x=xy(:,1); y=xy(:,2);
nvtx=length(x);
nel=length(ev(:,1));
lx=max(x)-min(x); ly=max(y)-min(y);
hx=max(diff(x)); hy=max(diff(y));
%
% Initialise global matrices
a = sparse(nvtx,nvtx);
n = sparse(nvtx,nvtx);
r = sparse(nvtx,nvtx);
%
% Set up 2x2 Gauss points
gpt=1.0e0/sqrt(3.0e0);
s(1) = -gpt;  t(1) = -gpt;   wt(1) = 1;
s(2) =  gpt;  t(2) = -gpt;   wt(2) = 1;
s(3) =  gpt;  t(3) =  gpt;   wt(3) = 1;
s(4) = -gpt;  t(4) =  gpt;   wt(4) = 1;

%
% Inner loop over elements
for ivtx = 1:4
   xl_v(:,ivtx) = x(ev(:,ivtx));
   yl_v(:,ivtx) = y(ev(:,ivtx));
   xsl(:,ivtx) = usol_c(ev(:,ivtx));
   ysl(:,ivtx) = vsol_c(ev(:,ivtx));
end
ae = zeros(nel,4,4);
ne = zeros(nel,4,4);
re = zeros(nel,4,4);
% Loop over 2x2 Gauss points
for igpt = 1:4
   sigpt=s(igpt);
   tigpt=t(igpt);
   wght=wt(igpt);
%  Evaluate derivatives etc
   [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
   u_x = zeros(nel,1); u_y=zeros(nel,1);
   for k=1:4
      u_x(:) = u_x(:) + xsl(:,k) .* phi(:,k);
      u_y(:) = u_y(:) + ysl(:,k) .* phi(:,k);
   end
   for j = 1:4
      for i = 1:4
          ae(:,i,j) = ae(:,i,j)  + dphidx(:,i).*dphidx(:,j) .* invjac(:);
          ae(:,i,j) = ae(:,i,j)  + dphidy(:,i).*dphidy(:,j) .* invjac(:);
          re(:,i,j) = re(:,i,j)  + phi(:,i).*phi(:,j) .* jac(:);
%%%       ne(:,i,j) = ne(:,i,j) + flowx(:) .* phi(:,i) .* dphidx(:,j);
%%%       ne(:,i,j) = ne(:,i,j) + flowy(:) .* phi(:,i) .* dphidy(:,j);
          ne(:,i,j)  = ne(:,i,j)  + wght*u_x(:).*phi(:,i).*dphidx(:,j);
          ne(:,i,j)  = ne(:,i,j)  + wght*u_y(:).*phi(:,i).*dphidy(:,j);
       end
	end
% End of Gauss point loop
end
%
% Perform assembly of global matrix  and source vector
for krow=1:4
   nrow=ev(:,krow);
   for kcol=1:4
      ncol=ev(:,kcol);
      a = a + sparse(nrow,ncol,ae(:,krow,kcol),nvtx,nvtx);
      r = r + sparse(nrow,ncol,re(:,krow,kcol),nvtx,nvtx);
      n = n + sparse(nrow,ncol,ne(:,krow,kcol),nvtx,nvtx);
   end
end
%
% Computation of element Peclet number (at the centroid)
% Rectangle specific calculation here
hx=abs(xl_v(:,2)-xl_v(:,1)); hy=abs(yl_v(:,3)-yl_v(:,2));
[jac,invjac,phi,dphidx,dphidy] = deriv(0,0,xl_v,yl_v);
u_x = zeros(nel,1); u_y=zeros(nel,1);
for k=1:4
   u_x(:) = u_x(:) + xsl(:,k) .* phi(:,k);
   u_y(:) = u_y(:) + ysl(:,k) .* phi(:,k);
end
flowx = u_x(:);   flowy = u_y(:);
flow_l2 = sqrt(flowx(:) .* flowx(:) + flowy(:) .* flowy(:));
if all(flowx==0), flow_h=hy;
elseif all(flowy==0), flow_h=hx;
else
   angle = atan(abs(flowy./flowx));
   flow_h = min([hx./cos(angle),hy./sin(angle)],[],2);
end
eph = flow_h;
epe = flow_h.*flow_l2/2;
epw = flow_l2;

author: johannesgerer

repo: jburkardt-m


text:
function [kSpace, evalMask] = extractKSpace_main(dData, iLC, iEvalInfoMask, measPara)
% extract kSpace data
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

% set repetitions
if(measPara.LCall(4) ~= 1)
    iRepetitions = 1:measPara.LCall(4);
else
    iRepetitions = 1;
end

% set averages
if(measPara.LCall(2) ~= 1)
    iAverages = 1:measPara.LCall(2);
else
    iAverages = 1;
end

nSeg = double(max(iLC(:,11)) - min(iLC(:,11)) + 1);

fprintf(1, 'Preprocessing Channels\n');
if(strcmp(measPara.dimension,'2D'))
    % cell: slice - cha
    % in cell: k_y - k_x - (t)
    if(measPara.LCall(1) == 1)
        % one 2D slice was acquired
        iSlices = 1;
    else
        % several 2D slices were acquired
        iSlices = 1:measPara.LCall(1);
    end

    kSpace = cell(length(iSlices),measPara.dim(5),length(iRepetitions),length(iAverages));
    evalMask = cell(length(iSlices),measPara.dim(5),length(iRepetitions),length(iAverages));

    % extract phase correction data
    if(~isempty(iEvalInfoMask))
        lPhaseCorr = false(size(iLC,1),1);
        for iPhaseCorr=1:length(lPhaseCorr)
            lPhaseCorr(iPhaseCorr) = bitget(iEvalInfoMask(iPhaseCorr,1),22);
        end
        iLC = iLC(~lPhaseCorr,:);
        dData = dData(~lPhaseCorr,:);
        if(~isempty(iEvalInfoMask))
            iEvalInfoMask = iEvalInfoMask(~lPhaseCorr,:);
        end
    end

    totalCount = length(iRepetitions)*length(iAverages)*measPara.dim(5)*length(iSlices)*measPara.dim(4)*nSeg;
    dispProgress('Extracting K-Space', 0, totalCount);
    for iRep=1:length(iRepetitions)
        dDataRep = dData(iLC(:,9) == iRepetitions(iRep)-1, :);
        iLCRep = iLC(iLC(:,9) == iRepetitions(iRep)-1, :);
        if(~isempty(iEvalInfoMask))
            iEvalInfoMaskRepLC = iEvalInfoMask(iLC(:,9) == iRepetitions(iRep)-1, :);
        end
        for iAvg=1:length(iAverages)
            dDataAvg = dDataRep(iLCRep(:,4) == iAverages(iAvg)-1, :);
            iLCAvg = iLCRep(iLCRep(:,4) == iAverages(iAvg)-1, :);
            if(~isempty(iEvalInfoMask))
                iEvalInfoMaskAvgLC = iEvalInfoMaskRepLC(iLCRep(:,4) == iAverages(iAvg)-1, :);
            end
            for iSli=1:length(iSlices)
                dDataSli = dDataAvg(iLCAvg(:,5) == iSlices(iSli)-1, :);
                iLCSli = iLCAvg(iLCAvg(:,5) == iSlices(iSli)-1, :);
                if(~isempty(iEvalInfoMask))
                    iEvalInfoMaskSliLC = iEvalInfoMaskAvgLC(iLCAvg(:,5) == iSlices(iSli)-1, :);
                else
                    iEvalInfoMaskSliLC = [];
                end

                [kSpace(iSli,:,iRep,iAvg), evalMask(iSli,:,iRep,iAvg)] = extractKSpace2D(dDataSli, iLCSli, iEvalInfoMaskSliLC, measPara.dim, [iRep, iAvg, iSli, totalCount, length(iAverages), length(iSlices)]);
                dispProgress('Extracting K-Space', ((iRep-1)*length(iAverages)*measPara.dim(5)*length(iSlices)*measPara.dim(4)*nSeg + (iAvg-1)*measPara.dim(5)*length(iSlices)*measPara.dim(4)*nSeg + iSli*measPara.dim(5)*measPara.dim(4)*nSeg)/totalCount);
            end
            dispProgress('Extracting K-Space', ((iRep-1)*length(iAverages)*measPara.dim(5)*length(iSlices)*measPara.dim(4)*nSeg + iAvg*measPara.dim(5)*length(iSlices)*measPara.dim(4)*nSeg)/totalCount);
        end
        dispProgress('Extracting K-Space', iRep/length(iRepetitions));
    end
    dispProgress('Extracting K-Space', 'Close');

elseif(strcmp(measPara.dimension,'3D'))
    % cell: 1 - cha
    % in cell: k_y - k_x - k_z

    kSpace = cell(1,measPara.dim(5),length(iRepetitions),length(iAverages));
    evalMask = cell(1,measPara.dim(5),length(iRepetitions),length(iAverages));
    totalCount = length(iRepetitions)*measPara.dim(5)*measPara.dim(3)*nSeg;
    dispProgress('Extracting K-Space', 0, totalCount);
    for iRep=1:length(iRepetitions)
        dDataRep = dData(iLC(:,9) == iRepetitions(iRep)-1, :);
        iLCRep = iLC(iLC(:,9) == iRepetitions(iRep)-1, :);
        if(~isempty(iEvalInfoMask))
            iEvalInfoMaskRepLC = iEvalInfoMask(iLC(:,9) == iRepetitions(iRep)-1, :);
        end
        for iAvg=1:length(iAverages)
            dDataAvg = dDataRep(iLCRep(:,4) == iAverages(iAvg)-1, :);
            iLCAvg = iLCRep(iLCRep(:,4) == iAverages(iAvg)-1, :);
            if(~isempty(iEvalInfoMask))
                iEvalInfoMaskAvgLC = iEvalInfoMaskRepLC(iLCRep(:,4) == iAverages(iAvg)-1, :);
            else
                iEvalInfoMaskAvgLC = [];
            end
            [kSpace(1,:,iRep,iAvg), evalMask(1,:,iRep,iAvg)] = extractKSpace3D(dDataAvg, iLCAvg, iEvalInfoMaskAvgLC, measPara.dim, [iRep,iAvg,totalCount]);
            dispProgress('Extracting K-Space', ((iRep-1)*length(iAverages) + iAvg)/(length(iRepetitions)*length(iAverages)));
        end
        dispProgress('Extracting K-Space', iRep/length(iRepetitions));
    end
    dispProgress('Extracting K-Space', 'Close');

elseif(strcmp(measPara.dimension,'4D'))
    % cell: 1 - cha
    % in cell: k_y - k_x - k_z - t

    % extract prescans/Navi for CS_Trufi
    if(max(iLC(:,1)) ~= min(iLC(:,1)))
        dData = dData(iLC(:,1) == measPara.dim(2),:);
        iLC = iLC(iLC(:,1) == measPara.dim(2),:);
        if(~isempty(iEvalInfoMask))
            iEvalInfoMask = iEvalInfoMask(iLC(:,1) == measPara.dim(2),:);
        end
    end
    if(strcmp(measPara.sequenceName,'CS_Trufi'))
        phaseLCcol = 10;
    else
        phaseLCcol = 8;
    end
    kSpace = cell(1,measPara.dim(5),length(iRepetitions),length(iAverages));
    evalMask = cell(1,measPara.dim(5),length(iRepetitions),length(iAverages));
    totalCount = length(iRepetitions)*measPara.dim(5)*measPara.dim(3)*nSeg*measPara.dim(4);
    dispProgress('Extracting K-Space', 0, totalCount);
    for iRep=1:length(iRepetitions)
        dDataRep = dData(iLC(:,9) == iRepetitions(iRep)-1, :);
        iLCRep = iLC(iLC(:,9) == iRepetitions(iRep)-1, :);
        if(~isempty(iEvalInfoMask))
            iEvalInfoMaskRepLC = iEvalInfoMask(iLC(:,9) == iRepetitions(iRep)-1, :);
        end
        for iAvg=1:length(iAverages)
            dDataAvg = dDataRep(iLCRep(:,4) == iAverages(iAvg)-1, :);
            iLCAvg = iLCRep(iLCRep(:,4) == iAverages(iAvg)-1, :);
            if(~isempty(iEvalInfoMask))
                iEvalInfoMaskAvgLC = iEvalInfoMaskRepLC(iLCRep(:,4) == iAverages(iAvg)-1, :);
            else
                iEvalInfoMaskAvgLC = [];
            end
            [kSpace(1,:,iRep,iAvg), evalMask(1,:,iRep,iAvg)] = extractKSpace4D(dDataAvg, iLCAvg, iEvalInfoMaskAvgLC, measPara.dim, [iRep,iAvg,totalCount], phaseLCcol);
            dispProgress('Extracting K-Space', ((iRep-1)*length(iAverages) + iAvg)/(length(iRepetitions)*length(iAverages)));
        end
        dispProgress('Extracting K-Space', iRep/length(iRepetitions));
    end
    dispProgress('Extracting K-Space', 'Close');

end

end

author: thomaskuestner


text:
function points = duplicateGerms(lx, ly, points)
% Duplicate germs for voronoi periodic boundary condition.
%
%   GERMS2 = duplicateGerms(LX, LY, GERMS);
%   LX and LY are vectors containing positions of pixels in an image
%   GERMS are points within thin image
%   The result is the set of germes repeated in each of the 8 directions
%   aroundthe central image. The number of returned points is 9 times the
%   number of input points.
%
%   Example
%   duplicateGerms([1 100], [1 100], rand(30, 2)*100);
%
%   See also
%     imvoronoi2d, imAWVoronoi, imPowerDiagram
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2009-05-29,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

N = size(points, 1);

% width of window in each dimension
if length(lx)>1
    width = [lx(end)-2*lx(1)+lx(2) ly(end)-2*ly(1)+ly(2)];
else
    width = [lx ly];
end

% duplicate the array of points with same coordinates
points = repmat(points, [9 1]);

% add x-shift for left points
for i=[1 2 3]
    points((i-1)*N+1:i*N, 1) = points((i-1)*N+1:i*N, 1) - width(1);
end

% add x-shift for right points
for i=[7 8 9]
    points((i-1)*N+1:i*N, 1) = points((i-1)*N+1:i*N, 1) + width(1);
end

% add y-shift for bottom points
for i=[1 4 7]
    points((i-1)*N+1:i*N, 2) = points((i-1)*N+1:i*N, 2) - width(2);
end

% add y-shift for top points
for i=[3 6 9]
    points((i-1)*N+1:i*N, 2) = points((i-1)*N+1:i*N, 2) + width(2);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function KLD = prtRvUtilWishartKld(q,Q,p,P)
% WISHARTKLD  Kulback Liebler Divergence between two Wishart densities
%       KLD(Q||P)
%
%   KL-Divergence of Normal, Gamma, Dirichlet and Wishart densities
%      Penny, 2001
%
% Syntax: KLD = wishartKLD(aQ,BQ,aP,BP)
%
% Inputs:
%   aQ - The strength parameter of the Q distribution
%   BQ - The mean parameter of the Q distribution
%   aP - The strength parameter of the P distribution
%   BP - The mean parameter of the P distribution
%
% Outputs:
%   KLD - The KLD for the Wishart distributions

d = size(Q,1);

dims = 1:d;

KLD = p/2*(prtUtilLogDet(Q) - prtUtilLogDet(P)) + q/2*(trace(inv(Q)*P) - d) + prtRvUtilGeneralizedGammaLn(p/2,d) - prtRvUtilGeneralizedGammaLn(q/2,d) + (q/2 - p/2)*sum(psi((q+1-dims)/2));

if KLD < 0
    KLD = 0; % This only happens in the range of 0;
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function A = meanDiscrete(s, hyp, x, i)

% Mean function for discrete inputs x. Given a function defined on the
% integers 1,2,3,..,s, the mean function is parametrized as:
%
% m(x) = mu_x,
%
% where mu is a fixed vector of length s.
%
% This implementation assumes that the inputs x are given as integers
% between 1 and s, which simply index the provided vector.
%
% The hyperparameters are:
%
% hyp = [ mu_1
%         mu_2
%         ..
%         mu_s ]
%
% Copyright (c) by Roman Garnett, 2014-08-14.
%
% See also COVDISCRETE.M, MEANFUNCTIONS.M.

if nargin==0, error('s must be specified.'), end           % check for dimension
if nargin<=2, A = num2str(s); return; end     % report number of hyperparameters
mu = hyp(:);
if nargin==3
  A = mu(x(:));                                                  % evaluate mean
else
  A = zeros(numel(x),1);                                            % derivative
  A(x==i) = 1;
end

author: benfulcher

repo: hctsa


text:
function diff_fab = usim(fa, f)
diff_fab = 0;
N = length(f);
for i = 1 : N
    fb = f{i};
    diff_ab = usimole(fa, fb);
    diff_fab = diff_fab + diff_ab;
end

author: decouples


text:
function bands = SpectrogramBands(spectrogram,frequencies,varargin)

%SpectrogramBands - Determine running power in physiological bands.
%
%  USAGE
%
%    bands = SpectrogramBands(spectrogram,frequencies,<options>)
%
%    spectrogram    spectrogram obtained using MTSpectrogram
%    frequencies    frequency bins for spectrogram
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'smooth'      smoothing for ratio (0 = no smoothing) (default = 2)
%     'delta'       set delta band (default = [0 4])
%     'theta'       set theta band (default = [7 10])
%     'spindles'    set spindle band (default = [10 20])
%     'lowGamma'    set low gamma band (default = [30 80])
%     'highGamma'   set high gamma band (default = [80 120])
%     'ripples'     set ripple band (default = [100 250])
%     'broadLow'    set broad low frequency band (default = [1 12])
%     'amyGamma'    set amygdala gamma band (default = [45 65])
%    =========================================================================
%
%  OUTPUT
%
%    bands.delta      delta power
%    bands.theta      theta power
%    bands.spindles   spindle power
%    bands.lowGamma   low gamma power
%    bands.highGamma  high gamma power
%    bands.ripples    ripple power
%    bands.broadLow   broad low frequency power
%    bands.amyGamma   amygdala gamma power
%
%    bands.ratios.hippocampus      theta/delta ratio
%    bands.ratios.cortex           [0.5 4.5]/[0.5 9] and [0.5 20]/[0.5 55]
%    bands.ratios.amygdala         gamma/broad low ratio
%
%    Heuristic ratios for cortex are from Gervasoni et al. (2004).
%
%  SEE
%
%    See also MTSpectrogram.

% Copyright (C) 2004-2014 by Michaël Zugaro, Gabrielle Girardeau
%
% 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.

% Defaults
smooth = 2;
delta = [0 4];
theta = [7 10];
spindles = [10 20];
lowGamma = [30 80];
highGamma = [80 120];
ripples = [100 250];
broadLow = [1 12];
amyGamma = [45 65];

% Check number of parameters
if nargin < 2,
  error('Incorrect number of parameters (type ''help SpectrogramBands'' for details).');
end

% Check parameter sizes
if ~isdmatrix(spectrogram),
	error('Parameter ''spectrogram'' is not a matrix (type ''help SpectrogramBands'' for details).');
end
if ~isdvector(frequencies) | length(frequencies) ~= size(spectrogram,1),
	error('Parameter ''frequencies'' is not a vector or does not match spectrogram size (type ''help SpectrogramBands'' for details).');
end

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i) ' is not a property (type ''help SpectrogramBands'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'smooth',
			smooth = varargin{i+1};
			if ~isdscalar(smooth,'>=0'),
				error('Incorrect value for property ''smooth'' (type ''help SpectrogramBands'' for details).');
			end
		case 'delta',
			delta = varargin{i+1};
			if ~isdvector(delta,'>=0','<','#2'),
				error('Incorrect value for property ''delta'' (type ''help SpectrogramBands'' for details).');
			end
		case 'theta',
			theta = varargin{i+1};
			if ~isdvector(theta,'>=0','<','#2'),
				error('Incorrect value for property ''theta'' (type ''help SpectrogramBands'' for details).');
			end
		case 'spindles',
			spindles = varargin{i+1};
			if ~isdvector(spindles,'>=0','<','#2'),
				error('Incorrect value for property ''spindles'' (type ''help SpectrogramBands'' for details).');
			end
		case 'lowgamma',
			lowGamma = varargin{i+1};
			if ~isdvector(lowGamma,'>=0','<','#2'),
				error('Incorrect value for property ''lowGamma'' (type ''help SpectrogramBands'' for details).');
			end
		case 'highgamma',
			highGamma = varargin{i+1};
			if ~isdvector(highGamma,'>=0','<','#2'),
				error('Incorrect value for property ''highGamma'' (type ''help SpectrogramBands'' for details).');
			end
		case 'ripples',
			ripples = varargin{i+1};
			if ~isdvector(ripples,'>=0','<','#2'),
				error('Incorrect value for property ''ripples'' (type ''help SpectrogramBands'' for details).');
			end
		case 'broadlow',
			broadLow = varargin{i+1};
			if ~isdvector(broadLow,'>=0','<','#2'),
				error('Incorrect value for property ''broadLow'' (type ''help SpectrogramBands'' for details).');
			end
		case 'amygamma',
			amyGamma = varargin{i+1};
			if ~isdvector(amyGamma,'>=0','<','#2'),
				error('Incorrect value for property ''amyGamma'' (type ''help SpectrogramBands'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help SpectrogramBands'' for details).']);
	end
end

n = size(spectrogram,1);

% Select relevant frequency bins
thetaBins = frequencies >= theta(1) & frequencies <= theta(2);
deltaBins = frequencies >= delta(1) & frequencies <= delta(2);
spindleBins = frequencies >= spindles(1) & frequencies <= spindles(2);
lowGammaBins = frequencies >= lowGamma(1) & frequencies <= lowGamma(2);
highGammaBins = frequencies >= highGamma(1) & frequencies <= highGamma(2);
rippleBins = frequencies >= ripples(1) & frequencies <= ripples(2);
broadLowBins = frequencies >= broadLow(1) & frequencies <= broadLow(2);
amyGammaBins = frequencies >= amyGamma(1) & frequencies <= amyGamma(2);

% Compute physiological bands
bands.theta = mean(spectrogram(thetaBins,:))';
bands.delta = mean(spectrogram(deltaBins,:))';
bands.spindles = mean(spectrogram(spindleBins,:))';
bands.lowGamma = mean(spectrogram(lowGammaBins,:))';
bands.highGamma = mean(spectrogram(highGammaBins,:))';
bands.ripples = mean(spectrogram(rippleBins,:))';
bands.broadLow = mean(spectrogram(broadLowBins,:))';
bands.amyGamma = mean(spectrogram(amyGammaBins,:))';

% Theta/delta ratio
bands.ratios.hippocampus = Smooth(bands.theta./(bands.delta+eps),smooth);
bands.ratio = bands.ratios.hippocampus;

% Heuristic ratios from Gervasoni et al. (2004)
range1 = [0.5 4.5];
bins1 = frequencies >= range1(1) & frequencies <= range1(2);
band1 =  mean(spectrogram(bins1,:))';
range2 = [0.5 9];
bins2 = frequencies >= range2(1) & frequencies <= range2(2);
band2 =  mean(spectrogram(bins2,:))';
bands.ratios.cortex(:,1) = Smooth(band1./(band2+eps),smooth);
bands.ratio1 = bands.ratios.cortex;

range1 = [0.5 20];
bins1 = frequencies >= range1(1) & frequencies <= range1(2);
band1 =  mean(spectrogram(bins1,:))';
range2 = [0.5 55];
bins2 = frequencies >= range2(1) & frequencies <= range2(2);
band2 =  mean(spectrogram(bins2,:))';
bands.ratios.cortex(:,2) = Smooth(band1./(band2+eps),smooth);
bands.ratio2 = bands.ratios.cortex(:,2);

% Amygdala ratio
bands.ratios.amygdala = Smooth(bands.amyGamma./(bands.broadLow+eps),smooth);
bands.ratio3 = bands.ratios.amygdala;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% eeg_eventformat() - Convert the event information of a dataset from struct
%                 to array or vice versa.
%
% Usage: >> [eventout fields] = eeg_eventformat( event, 'format', fields );
%
% Inputs:
%   event  - event array or structure
%   format - ['struct'|'array'] see below
%   fields - [optional] cell array of strings containing the names of
%            the event struct fields. If this field is empty, it uses
%            the following list for
%            the names of the fields { 'type' 'latency' 'var1' ...
%            'var2' ... }.
% Output:
%   eventout  - output event array or structure
%   fields    - output cell array with the name of the fields
%
% Event formats:
%   struct - Events are organised as an array of structs with at
%            least two fields ('type' and 'latency')
%            (Ex: reaction_time may be type 1).
%   array  - events are organized as an array, the first column
%            representing the type, the second the latency and the
%            other ones user-defined variables.
%
% Note: 1) The event structure is defined only for continuous data
%          or epoched data derived from continuous data.
%       2) The event 'struct' format is more comprehensible.
%          For instance, to see all the properties of event 7,
%          type >> EEG.event(7)
%          Unfortunately, structures are awkward for expert users to deal
%          with from the command line (Ex: To get an array of latencies,
%           >> [ EEG.event(:).latency ] )
%          In array format, the same information is obtained by typing
%           >> EEG.event(:,2)
%       3) This function automatically updates the 'eventfield'
%          cell array depending on the format.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 27 Jan 2002
%
% See also: eeglab(), pop_selectevent(), pop_importevent()

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

% 2/06/02 modifed header - sm & ad
% 2/08/02 add field input - ad
% 2/12/02 reprogrammed function using epochformat.m - ad

function [event, eventfield] = eeg_eventformat(event, format, fields);

if nargin < 2
   help eeg_eventformat;
   return;
end;

if exist('fields') ~= 1, fields = { 'type', 'latency' }; end;

[event eventfield] = eeg_epochformat( event, format, fields);

author: PatternRecognition

repo: OpenBMI


text:
function value = i4_is_prime ( n )

%*****************************************************************************80
%
%% I4_IS_PRIME reports whether an integer is prime.
%
%  Discussion:
%
%    A simple, unoptimized sieve of Erasthosthenes is used to
%    check whether N can be divided by any integer between 2
%    and SQRT(N).
%
%    Note that negative numbers, 0 and 1 are not considered prime.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer to be tested.
%
%    Output, logical VALUE, is TRUE if N is prime, and FALSE
%    otherwise.
%
  if ( n <= 0 )
    value = 0;
    return
  end

  if ( n == 1 )
    value = 0;
    return
  end

  if ( n <= 3 )
    value = 1;
    return
  end

  nhi = floor ( sqrt ( n ) );

  for i = 2 : nhi
    if ( mod ( n, i ) == 0 )
      value = 0;
      return
    end
  end

  value = 1;

  return
end


text:
function  checkout()
%Commit Checkout the currect directory to TortoiseSVN
% See also: svn
svn('checkout',pwd)

author: Sable

repo: mcbench-benchmarks


text:
function [fResult] = calc_log10poisspdf2(nX, fLambda)
    % Calculates the logarithm to the basis of 10 of the Poisson probability density function.
    % [fResult] = calc_log10poisspdf2(nX, fLambda)
    % --------------------------------------------------
    % Calculates the logarithm to the basis of 10 of the Poisson probability density function.
    % !!! Does not use zero values as input for fLambda!!!!!!
    %
    % Input parameters:
    %   nX          Parameter x (see help for 'poisspdf')
    %   fLambda     Parameter lambda (see help for 'poisspdf')
    %
    % Output parameters:
    %   fResult     Logarithm to the basis of 10 of the Poisson probability density
    %
    % J. Woessner, woessner@seismo.ifg.ethz.ch
    % updated: 03.12.2002

    % Create emtpy matrix for results
    fResult = zeros(size(nX));
    if isempty(fResult)
        return;
    end
    fResult(fLambda < 0) = NaN;

    % Select all computable elements
    vSel = (nX >= 0 & nX == round(nX) & fLambda > 0);

    % Adding of realmin to 0 cases is to get the effect of 0^0 = 1.
    if (any(vSel))
        fResult(vSel) = 1/log(10) * (-fLambda(vSel) + nX(vSel) .* log(fLambda(vSel) + realmin*(fLambda(vSel)==0)) ...
            - gammaln(nX(vSel) + 1));
        %fResult(vSel) = log10(poisspdf(nX(vSel),fLambda(vSel)));
    end
end

author: CelsoReyes


text:
report_this_filefun(mfilename('fullpath'));

str = [];
[newmatfile, newpath] = uiputfile([ hodo '*.dat'], 'Save As'); %Syntax change Matlab Version 7, no window positioning on macs

s = [a.Longitude   a.Latitude  a.Date  a.Date.Month  a.Date.Day  a.Magnitude  a.Depth a.Date.Hour a.Date.Minute a(:,10) a(:,11) a(:,12) a(:,13) a(:,14) ];
fid = fopen([newpath newmatfile],'w') ;
fprintf(fid,'%7.3f  %7.3f %6.2f %6.0f  %6.0f %6.1f %6.2f  %6.0f  %6.0f  %7.2f  %7.2f  %7.2f %7.2f  %7.2f\n',s');
fclose(fid);
clear s
return


text:
function save_mat(qual,outfilename)
% Saves a Seg_Quality object to a mat file.
% Inputs:
%   qual: The Seg_Quality object to be saved.
%   outfilename: The name for the mat output file.

assert(nargin==2)
assert(isa(qual,'Seg_Quality'))
assert(isa(outfilename,'char'))
assert(qual.validate())

warning('off','MATLAB:structOnObject')
qual = struct(qual);
warning('on','MATLAB:structOnObject')
save(outfilename,'qual');

author: nesl


text:
function [ cx, cy ] = zdrot ( n, cx, incx, cy, incy, c, s )

%*****************************************************************************80
%
%% ZDROT applies a complex plane rotation.
%
%  Discussion:
%
%    The cosine C and sine S are real and the vectors CX and CY are complex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, complex CX(*), one of the vectors to be rotated.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Input, complex CY(*), one of the vectors to be rotated.
%
%    Input, integer INCY, the increment between successive elements of CY.
%
%    Input, real C, S, parameters (presumably the cosine and sine of
%    some angle) that define a plane rotation.
%
%    Output, complex CX(*), one of the vectors to be rotated.
%
%    Output, complex CY(*), one of the vectors to be rotated.
%
  temp(1:n)               =   c * cx(1:incx:1+(n-1)*incx) ...
                            + s * cy(1:incy:1+(n-1)*incy);

  cy(1:incy:1+(n-1)*incy) = - s * cx(1:incx:1+(n-1)*incx) ...
                            + c * cy(1:incy:1+(n-1)*incy);

  cx(1:incx:1+(n-1)*incx) = temp(1:n);

  return
end


text:
% Implementation of a simple version of the Dijkstra shortest path algorithm
% Returns the distances from a single vertex to all others, doesn't save the path
%
% INPUTS: adjacency matrix, adj (nxn), start node s (index between 1 and n)
% OUTPUTS: shortest path length from the start node to all other nodes, 1xn
%
% Note: Works for a weighted/directed graph.
% GB: last updated, September 28, 2012

function d = simpleDijkstra(adj,s)

n=length(adj);
d = inf*ones(1,n); % distance s-all nodes
d(s) = 0;    % s-s distance
T = 1:n;    % node set with shortest paths not found yet

while not(isempty(T))
    [dmin,ind] = min(d(T));
    for j=1:length(T)
        if adj(T(ind),T(j))>0 && d(T(j))>d(T(ind))+adj(T(ind),T(j))
            d(T(j))=d(T(ind))+adj(T(ind),T(j));
        end
    end
    T = setdiff(T,T(ind));
end

author: aeolianine

repo: octave-networks-toolbox


text:
% Add a prior on the beta parameter to avoid low Signal to Noise Ratio
% problems.
function model = hsvargplvmControlSNR(model, meanSNR, layer, view, priorInfo, priorScale)

if nargin < 1, error('At least one argument needed!'); end
if nargin < 6, priorScale = 25; end
if nargin < 5, priorInfo = []; end
if nargin < 4 || isempty(view), view = 1; end
if nargin < 3 || isempty(layer), layer = model.H; end
% Where I want the expected value of my inv gamma if it was on SNR
if nargin < 2 || isempty(meanSNR), meanSNR = 150; end

if isempty(priorInfo)
    priorInfo.name = 'invgamma'; % What type of prior
    varData = var(model.layer{layer}.comp{view}.mOrig(:));
    meanB = meanSNR./varData;
    a=0.08;%1.0001; % Relatively large right-tail
    b=meanB*(a+1); % Because mode = b/(a-1)
    priorInfo.params = [a b];
end

model = hsvargplvmAddParamPrior(model, layer, 1, 'beta', priorInfo.name, priorInfo.params);
if ~isempty('priorScale')
    model.layer{layer}.comp{view}.paramPriors{1}.prior.scale = priorScale;
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function [X,Y,Z] = EquinNdes3D(N)

% function [X,Y,Z] = EquinNdes3D(N)
% Purpose: compute the equidistributed nodes on the reference tetrahedron

% total number of nodes
Np = (N+1)*(N+2)*(N+3)/6;

% 2) create equidistributed nodes on equilateral triangle
X = zeros(Np,1); Y = zeros(Np,1); Z = zeros(Np,1);

sk = 1;
for n=1:N+1
  for m=1:N+2-n
    for q=1:N+3-n-m
      X(sk) = -1 + (q-1)*2/N; Y(sk) = -1 + (m-1)*2/N; Z(sk) = -1 + (n-1)*2/N;
      sk = sk+1;
    end
  end
end
return;


text:
function x = cs_qrsol (A,b,order)                                           %#ok
%CS_QRSOL solve a sparse least-squares problem.
%   x = cs_qrsol(A,b) solves the over-determined least squares problem to
%   find x that minimizes norm(A*x-b), where b is a full vector and
%   A is m-by-n with m >= n.  If m < n, it solves the underdetermined system
%   Ax=b.  A 3rd input argument specifies the ordering method to use
%   (0: natural, 3: amd(A'*A)).  The default ordering is 3.
%
%   Example:
%       Prob = UFget ('HB/well1033') ; A = Prob.A ; [m n] = size (A) ;
%       b = rand (m,1) ;
%       x1 = cs_qrsol (A,b) ;
%       x2 = A\b ;
%       norm (x1-x2)
%
%   For this example, cs_qrsol is about 3 times faster than A\b in MATLAB 7.3.
%
%   See also CS_QR, CS_AMD, CS_LUSOL, CS_CHOLSOL, MLDIVIDE.

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

error ('cs_qrsol mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% I think this script was written before understanding that rudder derivatives
% are readily available from Tornado (?).
% Compute the first order derivative of all coefficients with respect to
% elevator (horizontal stabilizer rudder) deflection at different alphas. Yields
% derivative of coefficients per radian.
% Tornado\T135_export has to be in the path
% See also solverloop5.m line 312 on how to to a sweep over rudder deflections.

% @
% Copyright (C) 2017 Jonas Ruesch
%
% 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/>.
%
% @

currDir = pwd;
cd(fileparts(which(mfilename)));
scriptDir = pwd;

% Load the geometry and the initial state
load('../aircraft/ExperimentalCarrier.mat')
load('../aircraft/State_1.5alpha_12.5ms.mat')

% Must run in tornado root directory for loading airfoil profile later
cd('F:\svn\dev\matlab\tornado\T135_export')

% Alpha range
alphaStart = -10/180*pi;
alphaEnd = 20/180*pi;
numAlphas = 31; % must be > 1

resultsPos = [];
resultsNeg = [];
settings=config('startup');

%lattictype=1; %Standard VLM
lattictype=0;%Tornado freestream following wake VLM

alphaStep = (alphaEnd - alphaStart)/(numAlphas-1);
alpha = zeros(numAlphas,1);
CXdr = zeros(numAlphas, 1);
CYdr = zeros(numAlphas, 1);
CZdr = zeros(numAlphas, 1);
Cldr = zeros(numAlphas, 1);
Cmdr = zeros(numAlphas, 1);
Cndr = zeros(numAlphas, 1);

[n,m]=find(geo.flapped');
resetdelta=geo.flap_vector;
positiveFlap = 1*pi/180; % in radian
negativeFlap = -1*pi/180; % in radian
flapDiff = positiveFlap - negativeFlap;
rudder = 1; % the first rudder is on the horizontal stabilizer in ExperimentalCarrier geometry.

%% loop over alphas
for alphaIndex=1:numAlphas

    % Modify state
    alphaCurrent = alphaStart + (alphaIndex-1)*alphaStep;
    state.alpha = alphaCurrent;
    alpha(alphaIndex) = alphaCurrent;
    state.betha = 0.0;

    %% Positive Flap
    % Modify geometry
    geo.flap_vector(m(rudder),n(rudder))=positiveFlap;
    % Regenerate lattice.
    [lattice,ref]=fLattice_setup2(geo,state,lattictype);
    % Compute the solution
    [resultsPos]=solver9(resultsPos,state,geo,lattice,ref);
    [resultsPos]=coeff_create3(resultsPos,lattice,state,ref,geo);

    %% Negative Flap
    % Modify geometry
    geo.flap_vector(m(rudder),n(rudder))=negativeFlap;
    % Regenerate lattice.
    [lattice,ref]=fLattice_setup2(geo,state,lattictype);
    % Compute the solution
    [resultsNeg]=solver9(resultsNeg,state,geo,lattice,ref);
    [resultsNeg]=coeff_create3(resultsNeg,lattice,state,ref,geo);

    %% Derivative

    CXdr(alphaIndex, 1) = (resultsPos.CX - resultsNeg.CX)/flapDiff;
    CYdr(alphaIndex, 1) = (resultsPos.CY - resultsNeg.CY)/flapDiff;
    CZdr(alphaIndex, 1) = (resultsPos.CZ - resultsNeg.CZ)/flapDiff;
    Cldr(alphaIndex, 1) = (resultsPos.Cl - resultsNeg.Cl)/flapDiff;
    Cmdr(alphaIndex, 1) = (resultsPos.Cm - resultsNeg.Cm)/flapDiff;
    Cndr(alphaIndex, 1) = (resultsPos.Cn - resultsNeg.Cn)/flapDiff;

    %% Output
    disp(['alpha ' num2str(alphaCurrent*180/pi) ' deg']);
end

geo.flap_vector=resetdelta;

cd(currDir);

save

alphaDegrees = alpha*180/pi;

figure(23);
subplot(2,3,1);
plot(alphaDegrees,CXdr);
title('CXdr (longitudinal)');
xlabel('alpha');
zlabel('longitudinal coefficient derivative');

subplot(2,3,2);
plot(alphaDegrees,CYdr);
title('CYdr (lateral)');
xlabel('alpha');
zlabel('lateral coefficient derivative');

subplot(2,3,3);
plot(alphaDegrees, CZdr);
title('CZdr (lift)');
xlabel('alpha');
zlabel('lift force coefficient derivative');

subplot(2,3,4);
plot(alphaDegrees, Cldr);
title('Cldr (roll)');
xlabel('alpha');
zlabel('roll moment coefficient derivative');

subplot(2,3,5);
plot(alphaDegrees, Cmdr);
title('Cmdr (pitch)');
xlabel('alpha');
zlabel('pitch moment coefficient derivative');

subplot(2,3,6);
plot(alphaDegrees, Cndr);
title('Cndr (yaw)');
xlabel('alpha');
zlabel('yaw moment coefficient derivative');

author: jrgenerative

repo: fixed-wing-sim


text:
% LFMapRectifiedToMeasured - Applies a calibrated camera model to map desired samples to measured samples
%
% Usage:
%
%     InterpIdx = LFMapRectifiedToMeasured( InterpIdx, CalInfo, RectOptions )
%
% Helper function used by LFCalRectifyLF.  Based on a calibrated camera model, including distortion parameters and a
% desired intrinsic matrix, the indices of a set of desired sample is mapped to the indices of corresponding measured
% samples.
%
% Inputs :
%
%     InterpIdx : Set of desired indices, in homogeneous coordinates
%     CalInfo : Calibration info as returned by LFFindCalInfo
%     RectOptions : struct controlling the rectification process, see LFCalRectify
%
% Outputs:
%
%      InterpIdx : continuous-domain indices for interpolating from the measured light field
%
% User guide: LFToolbox.pdf
% See also: LFCalRectifyLF

% Copyright (c) 2013-2020 Donald G. Dansereau

function InterpIdx = LFMapRectifiedToMeasured( InterpIdx, CalInfo, RectOptions )

RectOptions = LFDefaultField( 'RectOptions', 'Precision', 'single' );

%---Cast the to the required precision---
InterpIdx = cast(InterpIdx, RectOptions.Precision);

%---Convert the index of the desired ray to a ray representation using ideal intrinsics---
InterpIdx = RectOptions.RectCamIntrinsicsH * InterpIdx;

%---Apply inverse lens distortion to yield the undistorted ray---
k1 = CalInfo.EstCamDistortionV(1); % r^2
k2 = CalInfo.EstCamDistortionV(2); % r^4
k3 = CalInfo.EstCamDistortionV(3); % r^6
b1 = CalInfo.EstCamDistortionV(4); % decentering of lens distortion
b2 = CalInfo.EstCamDistortionV(5); % decentering of lens distortion

InterpIdx(3:4,:) = bsxfun(@minus, InterpIdx(3:4,:), [b1; b2]); % decentering of lens distortion

%---Iteratively estimate the undistorted direction----
DesiredDirection = InterpIdx(3:4,:);
for( InverseIters = 1:RectOptions.NInverse_Distortion_Iters )
    R2 = sum(InterpIdx(3:4,:).^2); % compute radius^2 for the current estimate
    % update estimate based on inverse of distortion model
    InterpIdx(3:4,:) = DesiredDirection ./ repmat((1 + k1.*R2 + k2.*R2.^2 + k3.*R2.^3),2,1);
end
clear R2 DesiredDirection

InterpIdx(3:4,:) = bsxfun(@plus, InterpIdx(3:4,:), [b1; b2]); % decentering of lens distortion

%---Convert the undistorted ray to the corresponding index using the calibrated intrinsics---
% todo[optimization]: The variable InterpIdx could be precomputed and saved with the calibration
InterpIdx = CalInfo.EstCamIntrinsicsH^-1 * InterpIdx;

%---Interpolate the required values---
InterpIdx = InterpIdx(1:4,:); % drop homogeneous coordinates

repo: LFToolbox


text:
function R = InvIllco(A,iter)
%INVILLCO     Inverse of extremely ill-conditioned matrices
%
%   R = InvIllco(A,iter)
%
%On return, R is a cell array such that sum(R{i}) is an approximate inverse
%  of A and  I-R*A  is convergent.
%The input matrix A itself may be a cell array. This is convenient to store
%  not exactly representable input data in higher precision. Then R is an
%  approximate inverse of sum(A{i}).
%
%The parameter iter is optional. If specified, an extra iteration is
%  executed to produce  I-R*A  of order eps.
%
%Example:
%
%  n = 20;              % dimension os matrix
%  C = 1e50;            % anticipated condition number
%  A = randmat(n,C);    % ill-conditioned matrix
%  R = InvIllco(A);     % approximate inverse, stored in cell array R
%  norm(AccDot(R,A,-1,eye(n)),'fro')
%
%Implements algorithm InvIllco from
%  S.M. Rump: Inversion of extremely ill-conditioned matrices in floating-point,
%    Japan J. Indust. Appl. Math. (JJIAM), 26:249-277, 2009.
%
%Reference implementation! Slow due to interpretation!
%

% written  06/23/08     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest, warning
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end
  wng = warning;
  warning off

  if iscell(A)
    n = size(A{1},1);
    R = inv(A{1});
  else
    n = size(A,1);
    R = inv(A);
  end
  k = 1;
  iter = ( nargin==2 );
  while 1
    k = k+1;
    P = ProdKL(R,A,k,1);
    X = inv(P);
    while any(any(isinf(X)))
      X = inv(P.*(1+eps*randn(n)));
    end
    R = ProdKL(X,R,k,k);
    if norm(X,'fro')*norm(P,'fro')<.01/eps
      if iter, iter = 0; else break, end
    end
  end

  if rndold
    setround(rndold)
  end
  warning(wng)

repo: OpenMAS


text:
function SO3F = times(SO3F1,SO3F2)
% overloads |SO3F1 .* SO3F2|
%
% Syntax
%   sF = SO3F1 .* SO3F2
%   sF = a .* SO3F2
%   sF = SO3F1 .* a
%
% Input
%  SO3F1, SO3F2 - @SO3FunBingham
%  a - double
%
% Output
%  SO3F - @SO3Fun
%

if isnumeric(SO3F1)
  SO3F = SO3F2;
  SO3F.weight = SO3F.weight .* SO3F1;
  return
end

SO3F = times@SO3Fun(SO3F1,SO3F2);

end

author: mtex-toolbox

repo: mtex


text:
function [yd2] = mm22yd2(mm2)
% Convert area from square millimeters to square yards.
% Chad A. Greene 2012
yd2 = mm2*0.000001195990046301;