text:
function [vor, dist] = spm_voronoi(img, seeds, distance)
% Geodesic Discrete Voronoi Diagram - a compiled routine
% FORMAT [vor, dist] = spm_voronoi(img, seeds, distance)
%
% img       - binary image:  > 0  : inside
%                            <= 0 : outside
% seeds     - {n x 3} array of the n seeds positions [in voxels]
% distance  - type of chamfer distance to use ('d4', 'd8', 'd34' or 'd5711')
%             (default is 'd34')
%
% vor       - Geodesic Discrete Voronoi diagram
%             (label is equal to the index of the seed in 'seeds')
% dist      - Geodesic Distance map of img with seeds as objects
%
% Compute the geodesic discrete Voronoi Diagram of an image of labelled
% objects using front propagation. The distance map is also available
% on output.
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_voronoi.m 6079 2014-06-30 18:25:37Z spm $

%-This is merely the help file for the compiled routine
error('spm_voronoi.c not compiled - see Makefile')

author: spm


text:
% VL_KDTREEQUERY Query KD-tree
%   [INDEX, DIST] = VL_KDTREEQUERY(KDTREE, X, Y) computes the nearest
%   column of X to each column of Y (in Euclidean distance). KDTREE is
%   a forest of kd-trees build by VL_KDTREEBUILD(). X is a
%   NUMDIMENSIONS x NUMDATA data matrix of class SINGLE or DOUBLE with
%   the data indexed by the kd-trees (it must be the same data matrix
%   passed to VK_KDTREEBUILD() to build the trees). Y is the
%   NUMDIMENSIONS x NUMQUERIES matrix of query points and must have
%   the same class of X. INDEX is a 1 x NUMQUERIES matrix of class
%   UINT32 with the index of the nearest column of X for each column
%   of Y. DIST is a 1 x NUMQUERIES vector of class SINGLE or DOUBLE
%   (depending on the class of X and Y) with the corresponding squared
%   Euclidean distances.
%
%   [INDEX, DIST] = VL_KDTREEQUERY(..., 'NUMNEIGHBORS', NN) can be
%   used to return the N nearest neighbors rather than just the
%   nearest one. In this case INDEX and DIST are NN x NUMQUERIES
%   matrices. Neighbors are returned by increasing distance.
%
%   VL_KDTREEQUERY(..., 'MAXNUMCOMPARISONS', NCOMP) performs at most
%   NCOMP comparisons for each query point. In this case the result is
%   only approximate (i.e. approximated nearest-neighbors, or ANNs)
%   but the speed can be greatly improved.
%
%   Options:
%
%   NumNeighbors::
%     Sets the number of neighbors to compute for each query point (by
%     default 1).
%
%   MaxNumComparisons::
%     Sets the maximum number of comparisons per query point. The
%     special value 0 means unbounded. The default is 0.
%
%   See also: VL_KDTREEBUILD(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

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


text:
% CANFIS-ART Feed-Forward operation.
function y = canfis_art_forward(x,u2,v2,gamma,ThetaL4)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        								 							     				                                       %
%   			                 	NETWORK FUNCTIONALITY SECTION					   %
%      																					                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NumInVars    = size(x,1);
NumInTerms = size(u2,2);
NumRules     = NumInTerms;

% LAYER 1 - MF NODES
% Each node in this layer acts as a one-dimensional membership function.
 In1mat = x*ones(1,NumInTerms);

 % Matrix containing a copy of Out1 for each Layer 2 Term Node.
 Out1 = 1 - g(In1mat-v2,gamma) - g(u2-In1mat,gamma);

% LAYER 2 - PRECONDITION MATCHING OF FUZZY LOGIC RULES
% NumRules == NumInTerms
 Out2 = prod(Out1);
 S_2 = sum(Out2);

% LAYER 3 - NORMALIZATION NODES - All Node Activity Adds-up to unity.
if S_2~=0
     Out3 = Out2/S_2;
end

% LAYERS 4 - 5: CONSEQUENT NODES - SUMMING NODE
 Aux1 = [x; 1]*Out3;

 % New Input Training Data shaped as a column vector.
 a = reshape(Aux1,(NumInVars+1)*NumRules,1);
 y = ThetaL4'*a;

author: Sable

repo: mcbench-benchmarks


text:
function y = FourierShift2D(x, delta)
%
% y = FourierShift(x, [delta_x delta_y])
%
% Shifts x by delta cyclically. Uses the fourier shift theorem.
%
% Real inputs should give real outputs.
%
% By Tim Hutt, 26/03/2009
% Small fix thanks to Brian Krause, 11/02/2010

% The size of the matrix.
[N, M] = size(x);

% FFT of our possibly padded input signal.
X = fft2(x);

% The mathsy bit. The floors take care of odd-length signals.
x_shift = exp(-1i * 2 * pi * delta(1) * [0:floor(N/2)-1 floor(-N/2):-1]' / N);
y_shift = exp(-1i * 2 * pi * delta(2) * [0:floor(M/2)-1 floor(-M/2):-1] / M);

% Force conjugate symmetry. Otherwise this frequency component has no
% corresponding negative frequency to cancel out its imaginary part.
if mod(N, 2) == 0
	x_shift(N/2+1) = real(x_shift(N/2+1));
end
if mod(M, 2) == 0
	y_shift(M/2+1) = real(y_shift(M/2+1));
end

Y = X .* (x_shift * y_shift);

% Invert the FFT.
y = ifft2(Y);

% There should be no imaginary component (for real input
% signals) but due to numerical effects some remnants remain.
if isreal(x)
    y = real(y);
end

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef ActionExecNodeEnum < matlab.mixin.SetGet
    %ActionExecNodeEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        BeforeProp('Execute Actions Before Propagation')
        AfterProp('Execute Actions After Propagation');
    end

    properties
        name char = '';
    end

    methods
        function obj = ActionExecNodeEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function [listBoxStr, enums] = getListBoxStr()
            m = enumeration('ActionExecNodeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
            enums = m;
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('ActionExecNodeEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('ActionExecNodeEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ axisAngle ] = Rot2AxisAngle( Rot )
%ROT2AXISANGLE Summary of this function goes here
%   Detailed explanation goes here

    theta = acos((trace(Rot) - 1) / 2);

    vec = 1.0/(2*sin(theta));
    vec = vec * [Rot(3,2) - Rot(2,3), Rot(1,3) - Rot(3,1), Rot(2,1) - Rot(1,2)];
    axisAngle = vec * theta;
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function YLim= visutil_selectYLim(h, varargin)

props = {'Policy',          'auto'          '!CHAR(auto tightest tight)';
         'TightenBorder'     0.03        	'DOUBLE';
         'Symmetrize'       0               '!BOOL';
         'SetLim'           1               '!BOOL';
         };

if nargin==0,
  YLim= props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);
misc_checkType(h,'!GRAPHICS');

switch(opt.Policy),
 case 'auto',
  YLim= get(h, 'YLim');
 case 'tightest',
  visutil_backaxes(h);
  axis('tight');
  YLim= get(h, 'YLim');
 case 'tight',
  visutil_backaxes(h);
  axis('tight');
  yl= get(h, 'YLim');
  % add border not to make it too tight:
  yl= yl + [-1 1]*opt.TightenBorder*diff(yl);
  % determine nicer limits
  dig= floor(log10(diff(yl)));
  if diff(yl)>1,
    dig= max(1, dig);
  end
  YLim= [util_trunc(yl(1),-dig+1,'floor') util_trunc(yl(2),-dig+1,'ceil')];
end

if opt.Symmetrize,
  ma= max(abs(YLim));
  YLim= [-ma ma];
end

if opt.SetLim,
  set(h, 'YLim',YLim);
end

if nargout==0,
  clear YLim;
end

author: PatternRecognition

repo: OpenBMI


text:
function z = GetVel()
%
%
persistent Velp Posp

if isempty(Posp)
  Posp = 0;
  Velp = 80;
end

dt = 0.1;

v  = 0 +  10*randn;

Posp = Posp + Velp*dt;     % true position
Velp = 80 + v;             % true speed

z = Velp;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function h = ref_lconv(f,g,ctype)
%REF_LCONV  Reference linear convolution
%   Usage:  h=ref_lconv(f,g)
%
%   PCONV(f,g) computes the linear convolution of f and g.

% AUTHOR: Jordy van Velthoven

Lf = length(f);
Lg = length(g);

Lh = Lf+Lg-1;

f = [f; zeros(Lh - Lf, 1)];
g = [g; zeros(Lh - Lg, 1)];

h = zeros(Lf+Lg-1, 1);

switch(lower(ctype))
	case {'default'}
    for ii = 0 : Lh-1
      for jj = 0 : Lh-1
        h(ii+1)=h(ii+1)+f(jj+1)*g(mod(ii-jj,Lh)+1);
      end
    end
  case {'r'}
    for ii=0:Lh-1
      for jj=0:Lh-1
	      h(ii+1)=h(ii+1)+f(jj+1)*conj(g(mod(jj-ii, Lh)+1));
      end;
   	end;
  case {'rr'}
    for ii=0:Lh-1
      for jj=0:Lh-1
	      h(ii+1)=h(ii+1)+conj(f(mod(-jj, Lh)+1))*conj(g(mod(jj-ii,Lh)+1));
      end;
    end;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function CRB = crb_general_sto_1d(design, wavelength, doas, P, noise_var, snapshot_count)
%CRB_GENERAL_STO_1D CRB for general 1D arrays based on the stochastic
%(unconditional) model, in radians.
%Inputs:
%   design - Array design.
%   wavelength - Wavelength.
%   doas - DOA vector in radians.
%   P - Source covariance matrix. If all sources are uncorrelated and
%       shares the same power, you can just pass in a scalar. If all
%       sources are uncorrelated but have different powers, you can just
%       pass in a vector.
%   noise_var - Noise power.
%   snapshot_count - (Optional) number of snapshots. Default is one.
%Reference:
% [1] P. Stoica and A. Nehorai, "Performance study of conditional and
%     unconditional direction-of-arrival estimation," IEEE Transactions on
%     Acoustics, Speech and Signal Processing, vol. 38, no. 10,
%     pp. 1783-1795, Oct. 1990.
if design.dim ~= 1
    error('1D array expected.');
end
if nargin <= 5
    snapshot_count = 1;
end
m = design.element_count;
k = length(doas);
P = unify_source_power_matrix(P, k);
[A, D] = steering_matrix(design, wavelength, doas);
R = A*P*A' + noise_var * eye(m);
H = D'*(eye(m) - A/(A'*A)*A')*D;
CRB = real(H .* ((P*A'/R*A*P).'));
CRB = eye(k) / CRB * (noise_var / snapshot_count / 2);
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------

function stats = get_single_hier_stats(params,hier_id,gt_set)

if nargin<3
    gt_set = params.gt_set_pareto;
    res_file = params.files.pareto_singles{hier_id};
else
    res_file = strrep(params.files.pareto_singles{hier_id},[params.gt_set_pareto '_pareto_single'],[gt_set '_pareto_single']);
end

% Sampled number of candidates from each hierarchy
n_cands = [10:10:100,200:100:1000,2000:1000:10000];

for ii=1:params.n_r_cand
    n_max_cands(ii) = max(n_cands); %#ok<*AGROW>
end

% Are the results already computed?
if exist(res_file, 'file')
    load(res_file)
    disp(['Loaded: ' res_file '.'])
    recompute = 0;
else
    disp(['RECOMPUTING: ' res_file '.'])
    recompute = 1;
end

if recompute
    % Load which images to consider from the database (train, val, etc.)
    im_ids = database_ids(params.database,gt_set);

    for ll=1:params.n_r_cand
        % Store number of regions and gt_set
        stats(ll).n_cands = n_cands;
        stats(ll).gt_set = gt_set;
        stats(ll).num_objects = 0;
        stats(ll).obj_classes = [];
    end

    % Sweep all images
    num_images = length(im_ids);
    batches = 20;
    whichOne = imresize(1:batches, [1 num_images], 'nearest');

    for iii = 1:batches,
        indIII = find(whichOne == iii);
        clear Jaccards Inters False_pos False_neg True_areas Obj_classes;

        parfor ii_i = 1:length(indIII),
            ii = indIII(ii_i);
            % Read the UCM as a hierarchy
            hier = ucm2hier(fullfile(params.hier_dirs{hier_id}, [im_ids{ii} '.mat']));
            n_regs = hier.ms_struct(end).parent;

            % Get all pairs of neighboring leave regions
            [~, idx_neighbors] = seg2gridbmap(hier.leaves_part);
            K = max(idx_neighbors.matrix_max(:)) + 1;
            neigh_pairs = unique(idx_neighbors.matrix_min+K*idx_neighbors.matrix_max);
            neigh_pairs(neigh_pairs==0) = [];
            neigh_pairs_min = mod(neigh_pairs,K);
            neigh_pairs_max = (neigh_pairs-neigh_pairs_min)/K;

            if isrow(neigh_pairs_min)
                neigh_pairs_min = neigh_pairs_min';
            end
            if isrow(neigh_pairs_max)
                neigh_pairs_max = neigh_pairs_max';
            end

            % Get the 'n_cands' top candidates from each hierarchy

            % Singletons
            all_cands = {};
            if n_regs<=n_max_cands(1)
                all_cands{1} = (n_regs:-1:1)';
            else
                all_cands{1} = (n_regs:-1:(n_regs-n_max_cands(1))+1)';
            end

            % Pairs, triplets, etc.
            all_cands(2:params.n_r_cand) = ...
            mex_get_tree_cands(double(hier.leaves_part)-1, double(hier.ms_matrix)-1,...
                                 neigh_pairs_min-1, neigh_pairs_max-1,...
                                 n_max_cands(2:end));

            % Load ground truth
            curr_gt = get_ground_truth( params.database, im_ids{ii} );
            for ll = 1:params.n_r_cand,
                [Jaccards{ii_i}{ll}, Inters{ii_i}{ll}, False_pos{ii_i}{ll}, False_neg{ii_i}{ll}, True_areas{ii_i}{ll}, Obj_classes{ii_i}{ll}] = eval_cands(hier, all_cands{ll}, curr_gt); %#ok<ASGLU>
            end
        end
        fprintf('.');

        for ii_i = 1:length(indIII),
            ii = indIII(ii_i);
            for ll = 1:params.n_r_cand,
                jaccards = Jaccards{ii_i}{ll};
                inters = Inters{ii_i}{ll};
                false_pos = False_pos{ii_i}{ll};
                false_neg = False_neg{ii_i}{ll};
                true_areas = True_areas{ii_i}{ll};
                obj_classes = Obj_classes{ii_i}{ll};

                % Get best candidate at different number of candidates
                for jj=1:length(stats(ll).n_cands)
                    curr_n_regs = min(stats(ll).n_cands(jj), size(inters,2));
                    to_consider = 1:curr_n_regs;
                    stats(ll).all_n_masks(ii,jj) = length(to_consider);
                    if (stats(ll).all_n_masks(ii,jj)>0)
                        for kk=1:size(inters,1)
                            [stats(ll).max_J(stats(ll).num_objects+kk,jj), which_one] = max(jaccards(kk,to_consider));
                            stats(ll).max_indicator(stats(ll).num_objects+kk,jj) = to_consider(which_one);
                            stats(ll).max_fp(stats(ll).num_objects+kk,jj)        = false_pos(kk,to_consider(which_one));
                            stats(ll).max_fn(stats(ll).num_objects+kk,jj)        = false_neg(kk,to_consider(which_one));
                            stats(ll).max_inters(stats(ll).num_objects+kk,jj)    = inters(kk,to_consider(which_one));
                        end
                    end
                end
                stats(ll).obj_classes = [stats(ll).obj_classes; obj_classes(1:size(inters,1))'];
                stats(ll).num_objects = stats(ll).num_objects + size(jaccards,1);
            end
        end
    end

    for ll=1:params.n_r_cand,
        stats(ll).jaccard_object = mean(stats(ll).max_J);
        stats(ll).mean_n_masks   = mean(stats(ll).all_n_masks);

        % ----- Compute jaccard at pixel level (J_p) ----
        class_ids = unique(stats(ll).obj_classes);

        % Compute per-class statistics and then mean
        for ii=1:length(class_ids)
            curr_class = class_ids(ii);
            stats(ll).per_class_results{ii}.num_objects = sum(stats(ll).obj_classes==curr_class);

            stats(ll).per_class_results{ii}.max_fp     = stats(ll).max_fp(logical(stats(ll).obj_classes==curr_class),:);
            stats(ll).per_class_results{ii}.max_fn     = stats(ll).max_fn(logical(stats(ll).obj_classes==curr_class),:);
            stats(ll).per_class_results{ii}.max_inters = stats(ll).max_inters(logical(stats(ll).obj_classes==curr_class),:);
            stats(ll).per_class_results{ii}.max_J      = stats(ll).max_J(logical(stats(ll).obj_classes==curr_class),:);
            stats(ll).per_class_results{ii}.meanmax    = mean(stats(ll).per_class_results{ii}.max_J);

            stats(ll).per_class_results{ii}.global_fp      = sum(stats(ll).per_class_results{ii}.max_fp,1);
            stats(ll).per_class_results{ii}.global_fn      = sum(stats(ll).per_class_results{ii}.max_fn,1);
            stats(ll).per_class_results{ii}.global_inters  = sum(stats(ll).per_class_results{ii}.max_inters,1);

            % Compute per-class total inters, fp, fn
            stats(ll).per_class_results{ii}.global_J = ...
                 stats(ll).per_class_results{ii}.global_inters ./...
                (stats(ll).per_class_results{ii}.global_inters+...
                 stats(ll).per_class_results{ii}.global_fp    +...
                 stats(ll).per_class_results{ii}.global_fn);
        end

        % Compute global mean
        tmp = [];
        for ii=1:length(class_ids)
            tmp = [tmp; stats(ll).per_class_results{ii}.global_J];
        end
        stats(ll).jaccard_class = mean(tmp);
    end

    % Store
    % fprintf('Done with get_single_hier_stats.\n');
    save(res_file,'stats');
end

author: s-gupta

repo: rcnn-depth

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


text:
function [anim,HEye] = setFirstPRtToId( anim )
% Sets the first P or (R,t) of an Animation to eye(3,4)
%
% The function changes all the S and SBasis accordingly to preserve the
% projections.
%
% if R and t is defined, the rigid transform H is computed such that
%   [R(:,:,1),t(:,1)]*H=eye(3,4)
%   H is of the form [rotation,translation; 0 0 0 1]
%
% if R is not defined but P is, the homography H is computed such that:
%   P(:,:,1)*H=eye(3,4)
%
% USAGE
%  anim = anim.setFirstRToId()
%
% INPUTS
%  anim     - Animation object (help Animation for details)
%
% OUTPUTS
%  anim     - modified Animation
%  HEye     - [3x4] rigid transform if (R,t) are defined, or [3x4]
%             projective transform is P is defined
%
% EXAMPLE
%
% See also GENERATETOYANIMATION
%
% Vincent's Structure From Motion Toolbox      Version 3.1.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

isTorresaniModel = anim.nBasis==size(anim.l,1)+1;
if isempty(anim.R)
  HEye=anim.P(:,:,1)\eye(3,3);
  [disc,disc,V]=svd(anim.P(:,:,1));
  HEye(:,4)=V(:,4);
  % apply the homography to P
  anim.P=multiTimes(anim.P,HEye,1);
else
  HEye=[anim.R(:,:,1)',zeros(3,1);0,0,0,1];
  if anim.nBasis~=0 || isTorresaniModel
    HEye(1:3,4)=-anim.R(:,:,1)'*anim.t(:,1);
    % Re-generate the rotations and translations
    anim=subsasgn(anim,struct('type','.','subs','t'),...
      bsxfun(@minus,anim.t,reshape(multiTimes(anim.R, ...
      anim.R(:,:,1)'*anim.t(:,1),1),3,anim.nFrame)));
  end
  anim=subsasgn(anim,struct('type','.','subs','R'),...
    multiTimes(anim.R,anim.R(:,:,1)',1));
end

if anim.nBasis~=0
  % Re-generate the basis
  invHEye = inv(HEye);
  if isTorresaniModel
    SBasis=multiTimes(invHEye(1:3,1:3),anim.SBasis,1.2);
    SBasis(:,:,1)=bsxfun(@plus,SBasis(:,:,1),invHEye(1:3,4));
  else
    SBasis=multiTimes(invHEye(1:3,1:3),anim.SBasis,1.2);
  end

  % use subsasgn so that S is modified at the same time too
  anim=subsasgn(anim,struct('type','.','subs','SBasis'),SBasis);
else
  anim.S=normalizePoint(multiTimes(inv(HEye),normalizePoint(anim.S,-4),1.2),4);
end

author: vrabaud


text:
% Copyright 2013 Oliver Johnson, Srikanth Patala
%
% This file is part of MisorientationMaps.
%
%     MisorientationMaps 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.
%
%     MisorientationMaps 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 MisorientationMaps.  If not, see <http://www.gnu.org/licenses/>.

function rgb = colormap432(v)

% this is to adjust to the "correct" fundamental sector
sR = sphericalRegion([vector3d(0,1,-1),vector3d(1,-1,0),vector3d.Z]);
cs = crystalSymmetry('m-3m');
v = project2FundamentalRegion(v,cs,sR.center);

pts = reshape(double(v),[],3);

kmax = sqrt(2)-1;
theta = atan2(pts(:,3),pts(:,2));
t2ind2 = find(pts(:,1) > 1/3 & theta > atan((1-2*pts(:,1))./(pts(:,1))));

tempvar2 = pts(t2ind2,2).*(1-pts(t2ind2,1));
tempvar3 = pts(t2ind2,1).*(pts(t2ind2,2)+pts(t2ind2,3))./(tempvar2 + 1*(tempvar2 == 0));
pts(t2ind2,:) = [pts(t2ind2,1) tempvar3.*pts(t2ind2,2) tempvar3.*pts(t2ind2,3)];
g1 = rotvec2mat([1 0 0 -3*pi/8]);
pts = pts*g1;
pts = [pts(:,1) - kmax pts(:,2) pts(:,3)];

tempvar = (1 + pts(:,2).*tan(pi/8)./(pts(:,3) + 1*(pts(:,3)==0)));
pts = [pts(:,1) pts(:,2).*tempvar pts(:,3).*tempvar];

%%%%% Def Step 3 %%%%%
pts = [pts(:,1) pts(:,2)*cos(pi/8)/tan(pi/8) (pts(:,3) - pts(:,1)./cos(pi/8))];

%%%%% Def Step 4 %%%%%%
theta = atan2(-pts(:,1),pts(:,2));
pts = [pts(:,1).*(sin(theta) + abs(cos(theta))) ...
    pts(:,2).*(sin(theta) + abs(cos(theta))) pts(:,3)];

%%%% Def Step 5 %%%%%
theta = atan2(-pts(:,1),pts(:,2));
rad = hypot(pts(:,2),pts(:,1));
pts = [-rad.*sin(2*theta) rad.*cos(2*theta) pts(:,3)];
kmax = (sqrt(2)-1);tempk = cos(pi/8)/tan(pi/8);
pts(:,1) = pts(:,1)/kmax;pts(:,2) = pts(:,2)/kmax;pts(:,3) = pts(:,3)*tempk;
g1 = rotvec2mat([0 0 1 -pi/6]);
pts = pts*g1;

% S= hsvconereg(pts);
initS = hsvconereg(pts);
g1 = rotvec2mat([0 1 0 pi]);
initS = initS*g1;
initS = [initS(:,1) + 1 initS(:,2) initS(:,3)+1];
rgb = [initS(:,2) initS(:,3) initS(:,1)];
rgb(rgb > 1 & rgb - 1 <= eps) = 1;
rgb(rgb < 0 & rgb >= -10*eps) = 0;

end

% testing code
%v = 0.2*equispacedS2Grid('points',1000);

%rgb = colormap432(v);
%plot(v,rgb)

author: mtex-toolbox

repo: mtex

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


text:
% KM_DEMO_AKCCA Demonstration of Alternating Kernel Canonical Correlation
% Analysis algorithm for blind equalization of single-input multiple-output
% Wiener systems.
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2011.
%
% The algorithm in this file is based on the following publication:
% S. Van Vaerenbergh, J. Via and I. Santamaria, "Blind Identification of
% SIMO Wiener Systems based on Kernel Canonical Correlation Analysis",
% accepted for publication in IEEE Transactions on Signal Processing, 2013.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

close all
clear
rs = 1; % seed for random generator
rng('default')
rng(rs)

fprintf('\nAlternating kernel CCA for blind equalization of ')
fprintf('SIMO Wiener systems.\n');

%% SETUP PARAMETERS
pars_setup.model_num = 30244; % source model, see generate_data.m
pars_setup.model_L = 5; % channel length
pars_setup.data_type = 'gaussian'; % bits, gaussian
pars_setup.data_N = 256; % number of data points
pars_setup.data_SNR = 20; % signal-to-noise ratio
pars_setup.zf_k = 15; % number of observations to estimate zero-forcing equalizers
pars_setup.p = 4; % number of branches in SIMO Wiener system
pars_setup.verbose = false;

%% ALT-KCCA PARAMETERS
pars = pars_setup; % copy some setup parameters
pars.it_max = 100; % maximum number of iterations
pars.it_stop = 1E-10; % stop iteration if change in cost is smaller than this
pars.kernel.type = 'gauss'; % kernel type
pars.kernel.par = @(x) km_silverman(x); % kernel parameter, either a scalar or a function to determine it
pars.m = 1E-8; % number of KPCA autovectors / precision of ICD, or fraction of discarded signal energy
pars.reg = 1E-5; % regularization
pars.decomp = 'ICD'; % ICD or KPCA
pars.identical_nonlin = 0;	% boolean indicating common nonlinearity for each channel

%% PROGRAM
tic

N = pars_setup.data_N; p = pars_setup.p; SNR = pars_setup.data_SNR;
% GENERATE DATA
switch pars_setup.data_type,
    case 'gaussian'
        s = randn(N,1); s = s-mean(s); % source signal
    case 'bits'
        s = 2*round(rand(N,1))-1; % source signal
end
y = cell(p,1); x = cell(p,1); sigpow = 0;
f = @(x) tanh(0.8*x)+0.1*x; % nonlinearity
B1 = [0.6172   -0.8601    2.1383    0.4269   -1.3153
    0.6247    0.1532    0.9686   -0.5820   -0.4584
    0.3373   -0.1888   -1.4263    0.8060   -0.1740
    -0.0349   -0.6264   -0.2486    1.1975    1.2195
    -3.2957    0.9985   -0.3768    0.6139   -1.2011]/1.5; % linear filter
B = cell(p,1);
for i=1:p, % simulate Wiener system
    B{i} = B1(:,i);
    y{i} = filter(B{i},1,s);
    x{i} = f(y{i});
    sigpow = sigpow + x{i}'*x{i}/p/N;
end
noisepow = 10^(-SNR/10)*sigpow;
for i=1:p, x{i} = x{i} + sqrt(noisepow)*randn(N,1); end % add noise
data.s = s; data.B = B; data.y = y; data.x = x;

%% AKCCA algorithm
[vars,eval] = km_akcca(pars,data);

toc
%% OUTPUT

% calculate MSE or BER
switch pars_setup.data_type
    case 'gaussian'
        % scale signals and compare
        s1 = s(vars.ind_s); s2 = eval.s_est;
        s1_norm = s1/norm(s1); sc = s1_norm'*s2/(s2'*s2);
        err_vec = s1_norm - sc*s2;

        MSE = err_vec'*err_vec;
        sc = sc*norm(s1);

        fprintf('MSE = %.2f dB\n',10*log10(MSE));
    case 'bits'
        sc = 1;
        fprintf('BER = %.4f\n',eval.finalresult);
end

N = pars_setup.data_N;

% draw estimates of nonlinearities
figure
for i=1:p
    subplot(1,p,i); hold all
    xi = vars.x{i};
    y_est = vars.y_est{i};
    sci = y{i}'*y_est/(y_est'*y_est);

    plot(xi,data.y{i},'.b')
    plot(xi,sci*y_est,'.r');
    legend(sprintf('x%d vs true, unknown y%d',i),sprintf('x%d vs estimated y%d',i))
    xlabel('x')
end

% draw estimates of linear filters
figure
for i=1:p
    subplot(1,p,i); hold all
    stem(B{i},'ob','LineWidth',2,'MarkerSize',8);
    stem(vars.h{i}*norm(B{i})/norm(vars.h{i}),'xr','LineWidth',2,...
        'LineStyle','none','MarkerSize',10);

    xlabel('l');ylabel(sprintf('h_%d',i));
    legend({'real channel','estimated channel'})
end

% draw equalization result
figure; hold all
plot(s(vars.ind_s))
plot(sc*eval.s_est);
legend('source signal','recovered signal')

author: epfl-lasa


text:
function p = protate ( p, phi )

%*****************************************************************************80
%
%% PROTATE rotates a set of points by a given angle.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input/output, real P(NP,2), a set of points to be rotated.
%
%    Input, real PHI, the rotation angle, in radians.
%

%
%  Compute the rotation matrix.
%
  A = [ cos(phi), -sin(phi); ...
        sin(phi),  cos(phi) ];
%
%  Multiply the points by the rotation matrix.
%
  p = p * A;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ abd, rcond, info ] = cpbco ( abd, lda, n, m )

%*****************************************************************************80
%
%% CPBCO factors a complex hermitian positive definite band matrix.
%
%  Discussion:
%
%    The routine also estimates the condition number of the matrix.
%
%    If RCOND is not needed, CPBFA is slightly faster.
%
%    To solve A*X = B, follow CPBCO by CPBSL.
%
%    To compute inverse(A)*C, follow CPBCO by CPBSL.
%
%    To compute determinant(A), follow CPBCO by CPBDI.
%
%  Band storage:
%
%    If A is a hermitian positive definite band matrix,
%    the following program segment will set up the input.
%
%      m = (band width above diagonal)
%      do j = 1, n
%        i1 = max ( 1, j-m )
%        do i = i1, j
%          k = i-j+m+1
%          abd(k,j) = a(i,j)
%        end
%      end
%
%    This uses M+1 rows of A, except for the M by M
%    upper left triangle, which is ignored.
%
%  Example:
%
%    If the original matrix is
%
%      11 12 13  0  0  0
%      12 22 23 24  0  0
%      13 23 33 34 35  0
%       0 24 34 44 45 46
%       0  0 35 45 55 56
%       0  0  0 46 56 66
%
%    then N = 6, M = 2 and ABD should contain
%
%       *  * 13 24 35 46
%       * 12 23 34 45 56
%      11 22 33 44 55 66
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input/output, complex ABD(LDA,N); on input, the matrix to be factored.
%    The columns of the upper triangle are stored in the columns of ABD,
%    and the diagonals of the upper triangle are stored in the rows of ABD.
%
%    Input, integer LDA, the leading dimension of ABD.
%    LDA must be at least M+1.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer M, the number of diagonals above the main diagonal.
%    0 <= M < N.
%
%    Output, complex ABD(LDA,N);an upper triangular matrix R, stored in band
%    form, so that A = hermitian(R) * R.  If INFO ~= 0, the factorization
%    is not complete.
%
%    Output, real RCOND, an estimate of the reciprocal condition of
%    the matrix.  For the system A*X = B, relative perturbations in A and B
%    of size EPSILON may cause relative perturbations in X of size
%    (EPSILON/RCOND).  If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%    Output, integer INFO.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is not
%    positive definite.
%
%  Local parameters:
%
%    Local, complex Z(N), a work vector whose contents are usually
%    unimportant.  If A is singular to working precision, then Z is
%    an approximate null vector in the sense that
%    norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ).
%    If INFO ~= 0, Z is unchanged.
%

%
%  Find the norm of A.
%
  for j = 1 : n

    l = min ( j, m + 1 );
    mu = max ( m + 2 - j, 1 );
    z(j) = scasum ( l, abd(mu:mu+l-1,j), 1 );
    k = j - l;

    for i = mu : m
      k = k + 1;
      z(k) = real ( z(k) ) + cabs1 ( abd(i,j) );
    end

  end

  anorm = 0.0;
  for j = 1 : n
    anorm = max ( anorm, real ( z(j) ) );
  end
%
%  Factor.
%
  [ abd, info ] = cpbfa ( abd, lda, n, m );

  if ( info ~= 0 )
    return
  end
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where hermitian(R)*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve hermitian(R)*W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  for k = 1 : n

    if ( cabs1 ( z(k) ) ~= 0.0 )
      ek = csign1 ( ek, -z(k) );
    end

    if ( real ( abd(m+1,k) ) < cabs1 ( ek - z(k) ) )
      s = real ( abd(m+1,k) ) / cabs1 ( ek - z(k) );
      z(1:n) = z(1:n) * s;
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = - ek - z(k);
    s = cabs1 ( wk );
    sm = cabs1 ( wkm );
    wk = wk / abd(m+1,k);
    wkm = wkm / abd(m+1,k);
    j2 = min ( k + m, n );
    i = m + 1;

    if ( k+1 <= j2 )

      for j = k+1 : j2
        i = i - 1;
        sm = sm + cabs1 ( z(j) + wkm * conj ( abd(i,j) ) );
        z(j) = z(j) + wk * conj ( abd(i,j) );
        s = s + cabs1 ( z(j) );
      end

      if ( s < sm )
        t = wkm - wk;
        wk = wkm;
        i = m + 1;
        for j = k+1 : j2
          i = i - 1;
          z(j) = z(j) + t * conj ( abd(i,j) );
        end
      end

    end

    z(k) = wk;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
%
%  Solve R * Y = W.
%
  for k = n : -1 : 1

    if ( real ( abd(m+1,k) ) < cabs1 ( z(k) ) )
      s = real ( abd(m+1,k) ) / cabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
    end

    z(k) = z(k) / abd(m+1,k);
    lm = min ( k - 1, m );
    la = m + 1 - lm;
    lb = k - lm;
    t = -z(k);
    z(lb:lb+lm-1) = z(lb:lb+lm-1) + t * transpose ( abd(la:la+lm-1,k) );

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = 1.0;
%
%  Solve hermitian(R)*V = Y.
%
  for k = 1 : n

    lm = min ( k - 1, m );
    la = m + 1 - lm;
    lb = k - lm;
    z(k) = z(k) - abd(la:la+lm-1,k)' * transpose ( z(lb:lb+lm-1) );

    if ( real ( abd(m+1,k) ) < cabs1 ( z(k) ) )
      s = real ( abd(m+1,k) ) / cabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
      ynorm = s * ynorm;
    end

    z(k) = z(k) / abd(m+1,k);

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;
%
%  Solve R * Z = W.
%
  for k = n : -1 : 1

    if ( real ( abd(m+1,k) ) < cabs1 ( z(k) ) )
      s = real ( abd(m+1,k) ) / cabs1 ( z(k) );
      z(1:n) = z(1:n) * s;
      ynorm = s * ynorm;
    end

    z(k) = z(k) / abd(m+1,k);
    lm = min ( k - 1, m );
    la = m + 1 - lm;
    lb = k - lm;
    t = -z(k);
    z(lb:lb+lm-1) = z(lb:lb+lm-1) + t * transpose ( abd(la:la+lm-1,k) );

  end
%
%  Make ZNORM = 1.
%
  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end


text:
function qualityMeasure = imageQuality_srm(sr, gt)

% Making sure that only functions in the sr-metric folder can be called
addpath(genpath('../algorithms/qualityMeasures/sr-metric/'));

path = pwd;
cd(['..', filesep, 'algorithms', filesep, 'qualityMeasures', filesep, 'sr-metric']);

% Making sure that only functions in the sr-metric folder can be called
addpath('external/matlabPyrTools','external/randomforest-matlab/RF_Reg_C');

qualityMeasure = quality_predict(sr);

% Making sure that only functions in the sr-metric folder can be called
rmpath('external/matlabPyrTools','external/randomforest-matlab/RF_Reg_C');

cd(path);

% Making sure that only functions in the sr-metric folder can be called
rmpath(genpath('../algorithms/qualityMeasures/sr-metric/'));

author: thomas-koehler

repo: SupER


text:
% Test file for chebfun constructor (splitting).

function pass = test_constructor_splitting(pref)

% Grab some preferences:
if ( nargin == 0 )
    pref = chebfunpref();
end

seedRNG(6178);
tol = 1e9;

% Test SQRT(X) on [0 1]:
F1 = @sqrt;
f1 = chebfun(F1, [0, 1], pref, 'splitting', 'on', 'blowup', 'off');
xx1 = linspace(f1.domain(1)+eps, f1.domain(end)-eps, 100);
pass(1) = norm(feval(f1, xx1) - feval(F1, xx1), inf) < ...
    tol*max(eps*vscale(f1));

% Test SQRT(1-X) on [0 1]:
F2 = @(x) sqrt(1-x);
f2 = chebfun(F2, [0, 1], pref, 'splitting', 'on', 'blowup', 'off');
xx2 = linspace(f2.domain(1)+eps, f2.domain(end)-eps, 100);
pass(2) = norm(feval(f2, xx2) - feval(F2, xx2), inf) < ...
    tol*max(eps*vscale(f2));

% Test SQRT(1-X^2) on [-1 1]:
F3 = @(x) sqrt(1-x.^2);
f3 = chebfun(F3, [-1, 1], pref, 'splitting', 'on', 'blowup', 'off');
xx3 = linspace(f3.domain(1)+eps, f3.domain(end)-eps, 100);
pass(3) = norm(feval(f3, xx3) - feval(F3, xx3), inf) < ...
    tol*max(eps*vscale(f3));

% Test array-valued construction.  (Check for GitHub Issue #4.)
F4 = @(x) [sin(x) sign(x)];
f4 = chebfun(F4, [-1 1], pref, 'splitting', 'on', 'blowup', 'off');
xx4 = linspace(f4.domain(1)+eps, f4.domain(end)-eps, 100);
pass(4) = (norm(f4.domain - [-1 0 1], inf) < 10*eps) && ...
    (norm(feval(f4, xx4) - feval(F4, xx4), inf) < ...
    tol*max(eps*vscale(f4)));

% Check for issue with call to merge on a function with multiple breakpoints.
F5 = @(x) sign(x - 0.1).*abs(x + 0.2).*sin(3*x);
f5 = chebfun(F5, [-1 1], pref, 'splitting', 'on', 'blowup', 'off');
xx5 = linspace(f5.domain(1)+eps, f5.domain(end)-eps, 100);
pass(5) = (norm(f5.domain - [-1 -0.2 0.1 1], inf) < 10*eps) && ...
    (norm(feval(f5, xx5) - feval(F5, xx5), inf) < ...
    tol*max(eps*vscale(f5)));

%% Test a logical function:
f = chebfun(@(x) x > 0, [-1 1], 'splitting', 'on', pref);
x = chebfun('x', [-1 1], pref);
h = heaviside(x);
pass(6) = norm(f - h) < 10*eps;

% Test use of breakpoint detection in conjunction with construction from a cell
% array of function handles. (See issue #1151 on GitHub.)
f = chebfun({@(x) abs(x - 0.25), 0}, [0 0.5 1], 'splitting', 'on');
xx1 = linspace(0 + eps, 0.5 - eps, 20).';
err1 = norm(feval(f, xx1) - abs(xx1 - 0.25), inf);
xx2 = linspace(0.5 + eps, 1 - eps, 20).';
err2 = norm(feval(f, xx2), inf);
tol = 10*vscale(f)*eps;
pass(7) = max(err1, err2) < tol;

%% Test for function defined on unbounded domain:

% Function defined on [0 Inf]:

% Specify the domain:
dom = [0 Inf];
domCheck = [0 100];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

op = @(x) 0.75+sin(10*x)./exp(x);
f = chebfun(op, dom, 'splitting', 'on');
x = sort(x);
fVals = feval(f, x);
fExact = op(x);
err = fVals - fExact;
pass(8) = norm(err, inf) < 1e6*eps*vscale(f);

%% Test SPLITTING ON with BLOWUP == 1:
op = @(x)tan(x);
f = chebfun(op, [-4 4], 'splitting', 'on', 'blowup', 1);

% Specify the domain:
dom = [-4 -pi/2 pi/2 4];

% Generate random points to use as test values:
x1 = diff(dom(1:2)) * rand(100, 1) + dom(1);
x2 = diff(dom(2:3)) * rand(100, 1) + dom(2);
x3 = diff(dom(3:4)) * rand(100, 1) + dom(3);

err1 = op(x1)-f(x1);
err2 = op(x2)-f(x2);
err3 = op(x3)-f(x3);

err = [err1; err2; err3];
pass(9) = ( norm(err, inf) < 1e5*eps*vscale(f) );

%% Test for splitting on + unbndfun:

op = @(x)(sin(100*x)./exp(x.^2)+1).*(x.^2);
dom = [-inf inf];
dom_test = [-200 200];
x = diff(dom_test) * rand(100, 1) + dom_test(1);
f = chebfun (op, dom, 'exps', [2 2], 'splitting', 'on');
vals = f(x);
exact = op(x);
pass(10) = ( norm(vals-exact, inf) < 1e5*eps*vscale(f) );

% % Test X*LOG(X) on [0 1]:
% F4 = @(x) x.*log(x);
% f4 = chebfun(F4, [0, 1], pref, 'splitting', 'on', 'blowup', 'off');
% xx4 = linspace(f4.domain(1)+eps, f4.domain(end)-eps, 100);
% pass(4) = norm(feval(f4, xx4) - feval(F4, xx4), inf) < ...
%     tol*max(eps*vscale(f4));
%
% % Test (-X)*LOG(-X) on [-1 0]:
% F5 = @(x) (-x).*log(-x);
% f5 = chebfun(F5, [-1, 0], pref, 'splitting', 'on', 'blowup', 'off');
% xx5 = linspace(f5.domain(1)+eps, f5.domain(end)-eps, 100);
% pass(5) = norm(feval(f5, xx5) - feval(F5, xx5), inf) < ...
%     tol*max(eps*vscale(f5));
%
% % Test (1-X)*LOG(1-X) on [0 1]:
% F6 = @(x) (1-x).*log(1-x);
% f6 = chebfun(F6, [0, 1], pref, 'splitting', 'on', 'blowup', 'off');
% xx6 = linspace(f6.domain(1)+eps, f6.domain(end)-eps, 100);
% pass(6) = norm(feval(f6, xx6) - feval(F6, xx6), inf) < ...
%     tol*max(eps*vscale(f6));
% pass(6) = pass(6) && numel(f6.funs) < 50;
%
% % Test (1-X^2)*LOG(1-X^2) on [-1 1]:
% F7 = @(x) (1-x.^2).*log(1-x.^2);
% f7 = chebfun(F7, [-1, 1], pref, 'splitting', 'on', 'blowup', 'off');
% xx7 = linspace(f7.domain(1)+eps, f7.domain(end)-eps, 100);
% pass(7) = norm(feval(f7, xx7) - feval(F7, xx7), inf) < ...
%     tol*max(eps*vscale(f7));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef PTKLobesFromFissurePlaneOblique < PTKPlugin
    % PTKLobesFromFissurePlaneOblique. Plugin which is part of the lobar segmentation.
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKLobesFromFissurePlaneOblique is an intermediate stage in segmenting the
    %     lobes.
    %
    %     For more information, see
    %     [Doel et al., Pulmonary lobe segmentation from CT images using
    %     fissureness, airways, vessels and multilevel B-splines, 2012]
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Lobes Oblique'
        ToolTip = ''
        Category = 'Lobes'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(application, reporting)

            left_and_right_lungs = application.GetResult('PTKLeftAndRightLungs');
            fissure_plane = application.GetResult('PTKFissurePlaneOblique');
            lung_mask = application.GetResult('PTKLeftAndRightLungs');
            left_lung_template = application.GetTemplateImage(PTKContext.LeftLung).BlankCopy;
            right_lung_template = application.GetTemplateImage(PTKContext.RightLung).BlankCopy;
            results_left = PTKLobesFromFissurePlaneOblique.GetLeftLungResults(left_lung_template, lung_mask.Copy, fissure_plane.Copy, reporting);
            results_right = PTKLobesFromFissurePlaneOblique.GetRightLungResults(right_lung_template, lung_mask, fissure_plane, reporting);

            results = PTKCombineLeftAndRightImages(application.GetTemplateImage(PTKContext.LungROI), results_left, results_right, left_and_right_lungs);
            results.ImageType = PTKImageType.Colormap;
        end

        function results = GenerateImageFromResults(results, ~, ~)
        end
    end

    methods (Static, Access = private)
        function left_results = GetLeftLungResults(lung_template, lung_mask, fissure_plane, reporting)

            lung_mask.ChangeRawImage(uint8(lung_mask.RawImage == 2));

            lung_mask.ResizeToMatch(lung_template);

            fissure_plane.ResizeToMatch(lung_template);
            fissure_plane = find(fissure_plane.RawImage(:) == 4);

            left_results = PTKDivideVolumeUsingScatteredPoints(lung_mask, fissure_plane, 5, reporting);
            left_results.ChangeColourIndex(1, 5);
            left_results.ChangeColourIndex(2, 6);
        end

        function results_right = GetRightLungResults(lung_template, lung_mask, fissure_plane, reporting)
            lung_mask.ChangeRawImage(uint8(lung_mask.RawImage == 1));

            lung_mask.ResizeToMatch(lung_template);
            fissure_plane.ResizeToMatch(lung_template);
            fissure_plane_o = find(fissure_plane.RawImage(:) == 3);

            results_right = PTKDivideVolumeUsingScatteredPoints(lung_mask, fissure_plane_o, 5, reporting);
            results_right.ChangeColourIndex(2, 4);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function models = loadUncModels(modPath, organisms, objre, printLevel)
% This function loads and unconstrains metabolic models from a specific folder
%
% USAGE:
%
%   models = loadUncModels(modPath, organisms, objre)
%
% INPUTS:
%   organisms:           nx1 cell array cell array with names of organisms in the study
%   modPath:             char with path of directory where models are stored
%   objre:               char with reaction name of objective function of organisms
%   printLevel:          Verbose level (default: printLevel = 1)
%
% OUTPUT:
%   models:              nx1 cell array cell array with models of organisms in the study
%
% .. Author: Federico Baldini 2017-2018

    if ~exist('modPath', 'var') || ~exist(modPath, 'dir')  % check if the modPath is set
        error('modPath is not defined. Please set the path of the model directory.');
    else
        if ~exist(modPath, 'dir')
            error(['modPath (' modPath ') does not exist.']);
        end
    end

    % check if the organisms are defined
    if ~exist('organisms', 'var') || isempty(organisms)
        error('organisms are not defined. Please check your input file.');
    end

    if ~exist('objre', 'var')
        objre = {'EX_biomass(e)'};
        warning(['The default objective (objre) has been set to ' objre{1}]);
    end

    % check if a printLevel is present
    if ~exist('printLevel', 'var')
        printLevel = 0;
    end

    % add an extra file separator if missing
    if strcmpi(modPath(end), filesep)
        modPath = modPath(1:end-1);
    end

    environment = getEnvironment();

    models = {[]};  % empty cell array to be filled with models
    for i = 1:length(organisms)

        restoreEnvironment(environment);

        % adding the file extension
        fileEnding = '';
        tmpOrganism = organisms(i, 1);
        tmpOrganism = tmpOrganism{1};
        if ~strcmpi(tmpOrganism(end-4:end), '.mat')
            fileEnding = '.mat';
        end

        ldm = readCbModel([modPath filesep tmpOrganism fileEnding]);

        % removing possible constraints of the bacs
        selExc = findExcRxns(ldm);
        Reactions2 = ldm.rxns(find(selExc));
        allex = Reactions2(strmatch('EX', Reactions2));
        biomass = allex(strmatch(objre, allex));
        finrex = setdiff(allex, biomass);
        ldm = changeRxnBounds(ldm, finrex, -1000, 'l');
        ldm = changeRxnBounds(ldm, finrex, 1000, 'u');

        % removing blocked reactions from the bacs
        %BlockedRxns = identifyFastBlockedRxns(ldm,ldm.rxns, printLevel);
        %ldm= removeRxns(ldm, BlockedRxns);
        %BlockedReaction = findBlockedReaction(ldm,'L2')

        % creating array with models as required as input from the following functions
        models(i, 1) = {ldm};
    end
end

author: opencobra

repo: cobratoolbox


text:
function [mri] = align_neuromag2acpc(mri, method, template)

% ALIGN_NEUROMAG2ACPC performs an approximate alignment of the anatomical
% volume from NEUROMAG towards ACPC coordinates. Only the homogenous transformation
% matrix is modified and the coordsys-field is updated.
%
% Use as
%   mri = align_neuromag2acpc(mri)
%   mri = align_neuromag2acpc(mri, method)
%   mri = align_neuromag2acpc(mri, method, template)
%
% The first input argument is a FieldTrip MRI-structure, and the second optional
% argument specifies how the registration is to be done:
%   method = 0: only an approximate coregistration
%   method = 1: an approximate coregistration, followed by spm_affreg
%   method = 2: an approximate coregistration, followed by spm_normalise (default)
%
% When method = 1 or 2, an optional template filename can be specified, which denotes
% the filename of the target volume. This is required when running in deployed
% mode.
%
% See also ALIGN_CTF2ACPC, ALIGN_FSAVERAGE2MNI

if nargin<2
  method = 2;
end

%--------------------------------------------------------------------------
% do a first round of approximate coregistration using the predefined voxel
% locations of the fiducials and landmarks in the template T1 image from
% SPM

acpcvox2acpchead = [
  2     0     0   -92
  0     2     0  -128
  0     0     2   -74
  0     0     0     1
  ];

% these are the voxel indices of some points in the SPM canonical T1
acpcvox_Ac           = [46 64  37  1]';     % the anterior commissure
acpcvox_Ori          = [46 48  10  1]';     % approximately between the ears in T1.mnc
acpcvox_Nas          = [46 106 13  1]';     % approximately the nasion in T1.mnc
acpcvox_Lpa_canal    = [ 5 48  10  1]';     % Left ear canal
acpcvox_Rpa_canal    = [87 48  10  1]';     % Right ear canal

acpchead_Ac           = acpcvox2acpchead * acpcvox_Ac  ;
acpchead_Ori          = acpcvox2acpchead * acpcvox_Ori ;
acpchead_Nas          = acpcvox2acpchead * acpcvox_Nas ;
acpchead_Lpa_canal    = acpcvox2acpchead * acpcvox_Lpa_canal ;
acpchead_Rpa_canal    = acpcvox2acpchead * acpcvox_Rpa_canal ;

neuromagvox2neuromaghead  = mri.transform;
acpchead2neuromaghead     = ft_headcoordinates(acpchead_Nas(1:3), acpchead_Lpa_canal(1:3), acpchead_Rpa_canal(1:3), 'neuromag');

%neuromagvox2acpchead =  inv(acpchead2neuromaghead) * neuromagvox2neuromaghead;
neuromagvox2acpchead  =  acpchead2neuromaghead \ neuromagvox2neuromaghead;

% change the transformation matrix, such that it returns approximate SPM head coordinates
mri.transform     = neuromagvox2acpchead;
mri.vox2headOrig  = neuromagvox2neuromaghead;
mri.vox2head      = neuromagvox2acpchead;
mri.head2headOrig = acpchead2neuromaghead;
mri.coordsys      = 'acpc';

if method>0
  % this requires SPM to be on the path. However, this is not the proper place to
  % choose between SPM versions. The user can either use cfg.spmversion in a high-level
  % function, or has to add the path to the desired SPM version by hand.
  ft_hastoolbox('spm', -1);
end

%--------------------------------------------------------------------------
% Do a second round of affine registration (rigid body) to get improved
% alignment with ACPC coordinate system. this is needed because there may be
% different conventions defining LPA and RPA. The affine registration may
% fail however, e.g. if the initial alignment is not close enough. In that
% case SPM will throw an error

if method==1
  % use spm_affreg

  switch lower(spm('ver'))
    case 'spm2'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==1'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.mnc');
      end

    case 'spm8'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==1'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.nii');
      end

    case 'spm12'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==1'); end
      else
        template = fullfile(spm('Dir'),'toolbox','OldNorm','T1.nii');
        if ~exist('spm_affreg', 'file')
          addpath(fullfile(spm('Dir'),'toolbox','OldNorm'));
        end
      end
      fprintf('using ''OldNorm'' affine registration\n');

    otherwise
      ft_error('unsupported SPM version');
  end
  mri2 = ft_read_mri(template);

  tname1 = [tempname, '.img'];
  tname2 = [tempname, '.img'];
  V1 = ft_write_mri(tname1, mri.anatomy,  'transform', mri.transform,  'spmversion', spm('ver'), 'dataformat', 'nifti_spm');
  V2 = ft_write_mri(tname2, mri2.anatomy, 'transform', mri2.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

  % the below, using just spm_affreg does not work robustly enough in some cases
  flags.regtype = 'rigid';
  [M, scale]    = spm_affreg(V1,V2,flags);

  % some juggling around with the transformation matrices
  neuromagvox2acpchead2  = M \ V1.mat;
  acpchead2neuromaghead2 = neuromagvox2neuromaghead / neuromagvox2acpchead2;

  % update the transformation matrix
  mri.transform     = neuromagvox2acpchead2;

  % this one is unchanged
  mri.vox2headOrig  = neuromagvox2neuromaghead;

  % these are new
  mri.vox2head      = neuromagvox2acpchead2;
  mri.head2headOrig = acpchead2neuromaghead2;

  % delete the temporary files
  delete(tname1); delete(strrep(tname1, 'img', 'hdr'));
  delete(tname2); delete(strrep(tname2, 'img', 'hdr'));

elseif method==2
  % use spm_normalise

  switch lower(spm('ver'))
    case 'spm2'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==2'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.mnc');
      end

    case 'spm8'
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==2'); end
      else
        template = fullfile(spm('Dir'),'templates','T1.nii');
      end

    case 'spm12'
      % this uses the 'OldNorm' functionality, so the path needs to be
      % added, can only be done if non-deployed.
      if isdeployed
        if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==2'); end
      else
        template = fullfile(spm('Dir'),'toolbox','OldNorm','T1.nii');
        if ~exist('spm_normalise', 'file')
          addpath(fullfile(spm('Dir'),'toolbox','OldNorm'));
        end
      end
      fprintf('using ''OldNorm'' normalisation\n');

    otherwise
      ft_error('unsupported SPM version');
  end
  mri2 = ft_read_mri(template);

  tname1 = [tempname, '.img'];
  tname2 = [tempname, '.img'];
  V1 = ft_write_mri(tname1, mri.anatomy,  'transform', mri.transform,  'spmversion', spm('ver'), 'dataformat', 'nifti_spm');
  V2 = ft_write_mri(tname2, mri2.anatomy, 'transform', mri2.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

  flags.nits       = 0; %set number of non-linear iterations to zero
  flags.regtype    = 'rigid';
  params           = spm_normalise(V2,V1,[],[],[],flags);
  acpchead2neuromaghead2 = acpchead2neuromaghead*V1.mat*params.Affine/V2.mat;
  neuromagvox2acpchead2  = acpchead2neuromaghead2\neuromagvox2neuromaghead;

  % update the transformation matrix
  mri.transform     = neuromagvox2acpchead2;

  % this one is unchanged
  mri.vox2headOrig  = neuromagvox2neuromaghead;

  % these are new
  mri.vox2head      = neuromagvox2acpchead2;
  mri.head2headOrig = acpchead2neuromaghead2;

  % delete the temporary files
  delete(tname1); delete(strrep(tname1, 'img', 'hdr'));
  delete(tname2); delete(strrep(tname2, 'img', 'hdr'));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% SCRIPT TEST TO LOAD ALL ROBOTS

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

close all

fprintf('\nTHE DEMO LOADS ALL ROBOTS AND PLOTS THEM')

%there are eight possible solutions for the

%load robot parameters
% you can try different robots
robot=load_robot('ABB', 'IRB52');
robot=load_robot('ABB', 'IRB140');
robot=load_robot('ABB', 'IRB1600_6_120');
robot=load_robot('ABB', 'IRB1600_X145_M2004');
robot=load_robot('ABB', 'IRB1600ID');
robot=load_robot('ABB', 'IRB2400');
%robot=load_robot('ABB', 'IRB4400');
robot=load_robot('ABB', 'IRB4600');
robot=load_robot('ABB', 'IRB6620');
robot=load_robot('ABB', 'IRB6620LX');
robot=load_robot('ABB', 'IRB6650S_125_350');
robot=load_robot('ABB', 'IRB7600_150');
robot=load_robot('ABB', 'IRB7600_400_255_m2000');
robot=load_robot('ABB', 'IRB7600_500_230');
%ADEPT
robot=load_robot('ADEPT', 'Viper_s1700D');
%EPSON
robot=load_robot('EPSON', 'Prosix_C3_A601C');
%FANUC
robot=load_robot('FANUC', 'LR_MATE_200iC');
%KUKA
robot=load_robot('KUKA', 'KR5_2ARC_HW');
robot=load_robot('KUKA', 'KR5_arc');
robot=load_robot('KUKA', 'KR5_scara_R350_Z200');
robot=load_robot('KUKA', 'KR5_sixx_R650');
robot=load_robot('KUKA', 'KR5_sixx_R850');
robot=load_robot('KUKA', 'KR6_2');
robot=load_robot('KUKA', 'KR30_jet');
robot=load_robot('KUKA', 'KR90_R2700_pro');
robot=load_robot('KUKA', 'KR90_R3100_EXTRA');
robot=load_robot('KUKA', 'KR_16_arc_HW');
robot=load_robot('KUKA', 'KR_30_L16_2');
robot=load_robot('KUKA', 'KR_1000_1300_TITAN');
%MITSUBISHI
robot=load_robot('MITSUBISHI', 'pa-10');
robot=load_robot('MITSUBISHI', 'rv-6s');
%STANFORD ARM
robot=load_robot('example', 'stanford');
%STAUBLI
robot=load_robot('STAUBLI', 'RX160L');
%UNIMATE
robot=load_robot('UNIMATE', 'PUMA560');

%EXAMPLE ROBOTS
robot=load_robot('example', '2dofplanar');
robot=load_robot('example', '3dofplanar');
robot=load_robot('example', 'scara');
robot=load_robot('example', 'stanford');

repo: ARTE


text:
function c = plus(obj1,obj2)
% Implements the plus (+) operator on image_vector objects across voxels.
% Requires that each object has an equal number of columns and voxels
%
% :Examples:
% ::
%
%    c = dat1 + dat2;
%
% ..
%    Programmer Notes:
%    Created 3/14/14 by Luke Chang
% ..

error('This method is deprecated.  Use image_math instead.');

%Check if image_vector object
if ~isa(obj1,'image_vector') || ~isa(obj2,'image_vector')
    error('Input Data is not an image_vector object')
end

%Check number of rows
if size(obj1.dat,1)~=size(obj2.dat,1)
    error('number of voxels is different between objects.')
end

%Check number of columns
if size(obj1.dat,2)~=size(obj2.dat,2)
    error('number of observations is different between objects.')
end

c = obj1;
c.dat = obj1.dat + obj2.dat;

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function output = GetCustomDefaults(current,property,value)

%GetCustomDefaults - Get custom default value for a given function property.
%
%  This is an 'internal' function used by FMAToolbox. You should not need
%  to use it, unless you are developping new functions for this toolbox.
%
%  See Example below.
%
%  WARNING
%
%    Custom defaults should be allowed with great care, because they *implicitly*
%    change the default behavior of functions: because the actual function call
%    does not make it explicit that default values are being overriden, the same
%    code can yield different results for different users, or even for the same
%    user before and after she/he changes the custom defaults.
%
%  USAGE
%
%    output = GetCustomDefaults(current,property,value)
%
%    current        current value
%    property       property name
%    value          default value if custom default value is undefined
%
%  EXAMPLE
%
%    The function GetPositions has optional property-value pairs, including:
%
%     - 'mode' (default value 'clean')
%     - 'coordinates' (default value 'normalized')
%     - 'pixel' (no default value)
%     - 'distances' (default value [0 Inf])
%
%    Users may want to change default values, e.g. have 'coordinates'
%    set to 'video' whenever they call GetPositions without an explicit
%    value. GetPositions therefore includes the following code:
%
%      % Default values (customizable defaults must be empty at this point)
%      c = '';
%
%      % Parse parameters
%      (possibly set a value for c, if supplied among optional input parameters)
%
%      % Customizable defaults
%      c = GetCustomDefaults(c,'coordinates','normalized');
%
%    GetCustomDefaults will first test if c was already set (non-empty), in which
%    case c will remain unaltered. Otherwise, it will test if there is a user-defined
%    default value for 'coordinates', and issue a warning and set c to this value if
%    it exists, or set it to 'normalized' otherwise.
%
%    To define a custom default value, users would have to e.g. add the following
%    lines to their startup.m file:
%
%     global SETTINGS;
%     SETTINGS.GetPositions.coordinates = 'video';
%
%  NOTE
%
%    Do not forget to document when the default value for a property can be customized.
%    See the code for GetPositions for an example.

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

if nargin == 2,
	% Backward compatibility for former usage: output = GetCustomDefaults(property,value)
	value = property;
	property = current;
	warning('Deprecated usage (type ''help GetCustomDefaults'' for details).');
else
	% Non-empty input
	if ~isempty(current),
		output = current;
		return
	end
end

% Default value
output = value;

% Find out calling function name
stack = dbstack;
if length(stack) < 2, error('Stack error (sorry, cannot be more explicit)'); end
functionName = stack(2).name;

% Look for custom default
global SETTINGS;
if ~exist('SETTINGS'), return; end
if ~isfield(SETTINGS,functionName), return; end
p = getfield(SETTINGS,functionName);
if ~isfield(p,property), return; end
output = getfield(p,property);

% Warn user
if isdvector(output),
	d =  sprintf('%g ',output);
	d = d(1:end-1);
	if length(output) ~= 1,
		d = ['[' d '] '];
	else
		d = [d ' '];
	end
elseif isstring_FMAT(output),
	d = ['''' output ''' '];
else
	d = '';
end
warning(['Using custom default value ' d 'for ''' property ''' in function ''' functionName '''.']);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function varargout = htmlcal(startdate,stopdate,htmlfile,days)
% HTMLCAL Create an HTML calendar for a range of months.
%    HTMLCAL(FIRSTMONTH,LASTMONTH) creates an HTML file with one or more
%    monthly calendars. Use 'mm/yyyy' or 'mm-yyyy' strings to specify the
%    months. You will be prompted for the name of the output file, whose
%    name will be returned as output if any output variable is given.
%
%    HTMLCAL(FIRSTMONTH,LASTMONTH,HTMLFILE) also specifies the name of the
%    output file.
%
%    HTMLCAL(FIRSTMONTH,LASTMONTH,HTMLFILE,DAYS) allows you to specify
%    which days of the week are included, as a length-7 logical vector
%    starting with Sunday. The default is [0 1 1 1 1 1 0], for M-F.
%
%    HTMLCAL by itself prompts for the first and last month as well as the
%    file name, and offers to show the calendar in a web browser.
%
%    See also DATENUM, DATEVEC.

% Copyright 2005 by Toby Driscoll.
% Created 20 January 2005.
%         30 August 2007 -- documentation added, style tweaked

if nargin < 2
  if nargin < 1
    dv = datevec(date);
    default = sprintf('%i/%i',dv(2),dv(1));
    startdate = input(['First month (',default,')? '],'s');
    if isempty(startdate), startdate = default; end
    default = startdate;
    stopdate = input(['Last month (',default,')? '],'s');
    if isempty(stopdate), stopdate = default; end
  else
    stopdate = startdate;
  end
end

if nargin < 3 || isempty(htmlfile)
  [fn,pn] = uiputfile( {'*.html','HTML Files'}, ...
                       'Select output file');
  if isequal(fn,0) || isequal(pn,0)
    error('No output file selected.')
  else
    htmlfile = fullfile(pn,fn);
  end
end

if nargin < 4 || isempty(days)
  days = logical([ 0 1 1 1 1 1 0 ]);
elseif ischar(days)
  % This allows the "command syntax" to work consistently.
  days = eval(days);
end
days = days(:).';

mth = {'January','February','March','April','May','June',...
      'July','August','September','October','November','December'};
dow = {'Sunday','Monday','Tuesday','Wednesday','Thursday','Friday','Saturday'};

% Parse dates.
[m0,startdate] = strtok(startdate,'/-');
[y0,startdate] = strtok(startdate,'/-');
[m1,stopdate] = strtok(stopdate,'/-');
[y1,stopdate] = strtok(stopdate,'/-');

m0 = str2double(m0);  m1 = str2double(m1);
y0 = str2double(y0);  y1 = str2double(y1);
if any(isnan([m0 y0 m1 y1]))
  error('Invalid month specification.')
end

% Open output file and write headers.
outfile = CreateHTMLFile(htmlfile);

% Walk through month by month.
m1 = m1 + 12*(y1-y0);			% account for year-end boundaries
numdays = sum(days);
for m = m0:m1
  Month = rem(m-1,12) + 1;
  Year = y0 + floor(m/13);
  MonthName = mth{Month};

  cal = calendar(Year,Month);

  select = false(size(cal));
  select(:,days) = 1;                   % include these days
  select(~cal) = 0;			% but not those outside the
                                        % actual month

  % Begin table in output.
  fmt = '
%s %i
\n';
  fprintf(outfile,fmt,MonthName,Year);
  fprintf(outfile,...
          '
\n');

  % Header row.
  fprintf(outfile,'
\n  ');
  for d=find(days)
    fprintf(outfile,'%s	',round(100/numdays),dow{d});
  end
  fprintf(outfile,'\n
\n');

  % Body.
  for week = 1:size(cal,1)
    if any(select(week,:))
      fprintf(outfile,'
\n');
      for day = find(days)
        fprintf(outfile,'  ');
        if select(week,day)
          fprintf(outfile,'
%2i
',cal(week,day));
        end
        fprintf(outfile,'	\n');
      end
      fprintf(outfile,'
\n');
    end
  end

  % End the month
  fprintf(outfile,'
\n\n');
end

fclose(outfile);

% Wrap up
if nargin==0
  msg = 'Do you want to view the calendar in a browser window?';
  answer = questdlg(msg,'Success!','Yes','No','Yes');
  if strcmp(answer,'Yes')
    web(htmlfile,'-browser');
  end
end

if nargout > 0
  varargout{1} = htmlfile;
end

% --------------------------------------------------------------------
function outfile = CreateHTMLFile(htmlfile)

% Create an HTML output file with a stylesheet header. Returns file id
% used by fprintf.

[outfile,msg] = fopen(htmlfile,'w');
if ~isempty(msg)
  error(msg)
end

headstr = {
    ''
    ''
    ''
          };

fprintf(outfile,'%s\n',headstr{:});
fprintf(outfile,'\n');

author: Sable

repo: mcbench-benchmarks


text:
function check = weibull_discrete_check ( a, b )

%*****************************************************************************80
%
%% WEIBULL_DISCRETE_CHECK checks the parameters of the discrete Weibull CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 <= A <= 1.0,
%    0.0 < B.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( a < 0.0 | 1.0 < a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WEIBULL_DISCRETE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  A < 0 or 1 < A.\n' );
    check = 0;
    return
  end

  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WEIBULL_DISCRETE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= 0.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end


text:
% std_substudy()  - create a sub-STUDY set by removing datasets, conditions, groups, or
%                   subjects.
% Usage:
%  >> [ STUDY ALLEEG ] = std_substudy(STUDY, ALLEEG, 'key', 'val', ...);
%
% Optional Inputs:
%   STUDY                - existing study structure.
%   ALLEEG               - vector of EEG dataset structures to be included in the STUDY.
%
% Optional Inputs:
%   'dataset'   - [integer array] indices of dataset to include in sub-STUDY
%                 Default is all datasets.
%   'subject'   - [cell array] name of subjects to include in sub-STUDY.
%                 Default is all subjects.%
%   'condition' - [cell array] name of conditions to include in sub-STUDY
%                 Default is all conditions.
%   'group'     - [cell array] name of gourps to include in sub-STUDY
%                 Default is all groups.
%   'rmdat'     - ['on'|'off'] actually remove datasets 'on', or simply
%                 remove all references to these datasets for channels and
%                 clusters ('off').
%
% Example:
%    % create a sub-STUDY using only the first 3 subjects
%    % WARNING: make sure your STUDY is saved before creating a sub-STUDY
%    [STUDY ALLEEG] = std_substudy(STUDY, ALLEEG, 'subject', STUDY.subjects(1:3));
%
% Authors: Arnaud Delorme, CERCO/CNSR & SCCN, INC, UCSD, 2009-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, arno@sccn.ucsd.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

function [ STUDY ALLEEG ] = std_substudy(STUDY, ALLEEG, varargin);

if nargin < 3
    help std_substudy;
    return;
end

opt = finputcheck(varargin, { 'condition' 'cell' {}      {};
                              'dataset'   'integer' {}   [];
                              'group'     'cell' {}      {};
                              'rmdat'     'string' { 'on' 'off' }      'on';
                              'subject'   'cell' {}      {} }, 'std_substudy');
if isstr(opt), return; end;

% find datasets to remove
% -----------------------
tagdel = [];
if ~isempty(opt.subject)
    for index = 1:length(STUDY.datasetinfo)
        if ~strmatch(STUDY.datasetinfo.subject, opt.subject, 'exact')
            tagdel = [ tagdel index ];
        end;
    end;
end;
if ~isempty(opt.condition)
    for index = 1:length(STUDY.datasetinfo)
        if ~strmatch(STUDY.datasetinfo.condition, opt.condition, 'exact')
            tagdel = [ tagdel index ];
        end;
    end;
end;
if ~isempty(opt.group)
    for index = 1:length(STUDY.datasetinfo)
        if ~strmatch(STUDY.datasetinfo.group, opt.group, 'exact')
            tagdel = [ tagdel index ];
        end;
    end;
end;
if ~isempty(opt.dataset)
    tagdel = [ tagdel setdiff([1:length(ALLEEG)], opt.dataset) ];
end;
tagdel = unique(tagdel);

% find new dataset indices
% ------------------------
alldats = [1:length(ALLEEG)];
if strcmpi(opt.rmdat, 'on')
    alldats(tagdel) = [];
    for index = 1:length(ALLEEG)
        tmp = find(alldats == index);
        if isempty(tmp), tmp = NaN; end;
        datcoresp(index) = tmp;
    end;
    ALLEEG(tagdel) = [];
    STUDY.datasetinfo(tagdel) = [];
    for index = 1:length(STUDY.datasetinfo)
        STUDY.datasetinfo(index).index = index;
    end;
else
    alldats(tagdel) = [];
    for index = 1:length(ALLEEG)
        tmp = find(alldats == index);
        if isempty(tmp), tmp = NaN; else tmp = index; end;
        datcoresp(index) = tmp;
    end;
end;

% check channel consistency
% -------------------------
for i = 1:length(STUDY.changrp)
    for c = 1:size(STUDY.changrp(i).setinds,1)
       for g = 1:size(STUDY.changrp(i).setinds,2)
           newinds = datcoresp(STUDY.changrp(i).setinds{c,g});
           nonnans = find(~isnan(newinds));
           STUDY.changrp(i).setinds{c,g} = newinds(nonnans);
           STUDY.changrp(i).allinds{c,g} = STUDY.changrp(i).allinds{c,g}(nonnans);
       end;
    end;
end;

% check cluster consistency
% -------------------------
for index = 1:length(STUDY.cluster)
    STUDY.cluster(index).sets(:) = datcoresp(STUDY.cluster(index).sets(:));
    for i = size(STUDY.cluster(index).sets,2):-1:1
        if all(isnan(STUDY.cluster(index).sets(:,i)))
            STUDY.cluster(index).sets(:,i) = [];
            STUDY.cluster(index).comps(:,i) = [];
        end;
    end;
    [tmp STUDY.cluster(index).setinds STUDY.cluster(index).allinds] = std_setcomps2cell(STUDY, STUDY.cluster(index).sets, STUDY.cluster(index).comps);
end;

STUDY = std_reset(STUDY);
STUDY = std_checkset(STUDY, ALLEEG);

author: PatternRecognition


text:
function [DF] = df_build(V,C)
% DF_BUILD Build a volumetric distance field and spatial index that can be used
% to speed up approximate NN queries in 3D This code is not optimized, it is
% very slow for big datasets
%
% [DF] = df_build(V,C)
%
% Inputs:
%   V  Coordinates of the points
%   C  Column vector that defines the number of cells for every
%      dimension
% Output:
%   DF  structure to be passed to df_query to execute queries

if ~exist('C','var')
    C=[10;10;10];
end

MIN = min(V,[],1);
MAX = max(V,[],1);

S = (MAX-MIN)./C';
S(S==0) = 0.0000001;

D = zeros(C'+1);
N = zeros(C'+1);

for x=1:size(D,1)
    progressbar(x,size(D,1));
    for y=1:size(D,2)
        for z=1:size(D,3)
            i = [x,y,z];
            p = MIN + (i-1).*S;
            [D(x,y,z) N(x,y,z)] = min(normrow(repmat(p,size(V,1),1)-V));
        end
    end
end

% Prepare DF struct
DF.MIN = MIN;
DF.MAX = MAX;
DF.D = D;
DF.N = N;
DF.C = C;
DF.S = S;
end

author: alecjacobson

repo: gptoolbox


text:
function [error_M,error_A,P] = compare_2_endmembers(M1,M2,A1,A2,m,n,...
    names,wl,show_figure,options)
%COMPARE_ENDMEMBERS Summary of this function goes here
%   Detailed explanation goes here
if nargin < 9
    show_figure = 0;
end

if nargin < 10
    options = [];
end

A1 = double(A1);

% use_uncertainty = 0;
var_dirs = [];
var_amts = [];

permute_criterion = 'auto';

if isstruct(options)
    arg_set = fieldnames(options);
    for i = 1:length(arg_set)
        eval([arg_set{i},'=options.',arg_set{i},';']);
    end
end

use_uncertainty = ~isempty(var_amts);

if isempty(M1) || isempty(A1) || size(A1,2) < size(A2,2)
    show_endmembers(M2,wl,names);
    show_abundances(A2,m,n);
    if use_uncertainty
        show_uncertainty_range([], M2, var_amts, var_dirs, wl, names, options);
    end

    error_M = [];
    error_A = [];
    P = eye(size(M2,1));
    return;
end

[M,B] = size(M1);

is_real_dataset = 0;
if length(unique(A1(:))) == 2
    is_real_dataset = 1;
end

if strcmp(permute_criterion, 'auto')
    if is_real_dataset
        permute_criterion = 'abundance';
    else
        permute_criterion = 'endmember';
    end
end

%% Permute the endmembers and abundances from the algorithm to accord with
% the ground truth ones.
if strcmp(permute_criterion, 'abundance')
    P = permute_abundances(A1,A2);
    M2_1 = P*M2;
elseif strcmp(permute_criterion, 'endmember')
    [P,M2_1] = permute_endmembers(M1,M2);
end
A2_1 = (P*A2')';

if use_uncertainty
    var_amts = P*var_amts;
    var_dirs = P*var_dirs;
end

%% Plot 2 groups of endmembers and show uncertainty range if possible
if show_figure
    show_uncertainty_range(M1, M2_1, var_amts, var_dirs, wl, names, options);
end

%% Calculate the error
error_A = calc_abundance_error(A1,A2_1);
error_M = nanmean(sqrt(nanmean(abs(M1-M2_1).^2, 2)));

if show_figure
    disp(['Error for A: ',num2str(error_A)]);
    disp(['Error for M: ',num2str(error_M)]);
end

if show_figure
    show_abundances(A2_1,m,n);
end

%% If size(A1,2) < size(A2,2), make the permutation matrix P square
if size(P,1) < size(P,2)
    missings = [];
    for i = 1:size(P,2)
        if all(P(:,i) == 0)
            missings = [missings,i];
        end
    end
    new_P1 = zeros(length(missings), size(P,2));
    for i = 1:length(missings)
        new_P1(i, missings(i)) = 1;
    end
    P = [P;new_P1];
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
% moving camera bundle adjustment problem

ncams = 10;  % number of cameras
npts = 110;  % number of landmark points
visprob = 0.65;  % probability of camera seeing a landmark
pixnoise = 0.5; % standard deviation of Gaussian noise added to camera projections
Tshift = SE3(-0.2, 0, 0);  % horizontal shift of camera at each view

% create a camera
cam = CentralCamera('default', 'noise', pixnoise)

% setup a bundle adjustment problem
ba = BundleAdjust(cam);

% create the camera nodes, all at the origin, and keep their handles ch(i)
for i=1:ncams
    if i == 1
        ch(i) = ba.add_camera( SE3, 'fixed' );
    else
        ch(i) = ba.add_camera( SE3 );
    end
end

% create a working volume containing npts random points
% x -3 -> 1
% y -2 -> 2
% z  4 -> 8
randinit
P = bsxfun(@plus, 2 * 2*(rand(3, npts) - 0.5), [-1, 0 , 6]');

% create the landmark nodes and keep their handles lh(j)
for j=1:numcols(P)
    lh(j) = ba.add_landmark( P(:,j) );
end

% slide the camera in the x-direciton
T = SE3;
for i=1:ncams
    % project all landmarks for this camera position
    [p, visible] = cam.project(P, 'pose', T);

    % find the subset of points that are visible
    for j=find(visible)'
        % add to the problem if visible
        if rand < visprob % with a probability
            ba.add_projection(ch(i), lh(j), p(:,j));
        end
    end
    T = T .* Tshift; % shift the camera
end

% display the problem summary
ba

figure(1)
ba.plot()

pause

% get the state vector
X = ba.getstate();

% get the initial error
ba.errors(X)

% display the Hessian
figure(2)
ba.spyH(X)
title('Hessian matrix sparsity')

% solve the problem
baf = ba.optimize(X);

figure(3)
baf.plot()

author: petercorke

repo: machinevision-toolbox-matlab


text:
function doseViewer3D(command, varargin);
%"doseFigure3D"
%   Creates a figure containing a 3D dose browser for plancheck.
%
%JRA 07/13/05
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.
%
%Usage:
%   3DDoseViewer

global planC;
global stateS;
indexS = planC{end};

if ~exist('command')
    command = 'init';
end

hFig = findobj('Tag', 'CERR_3DDoseViewer');

switch upper(command)
    case 'INIT'
        if ~isempty(hFig)
            delete(hFig);
        end

        % check if atleast one dose exists
        if length(planC{indexS.dose})<1
            msgbox('No dose exist. Please calculate and store dose first')
            return;
        else
            doseSet=1;
        end

        screenSize = get(0,'ScreenSize');
        units = 'pixels';
        y = 480;
        x = 640;

        dy = floor((y-30)/2);
        dx = floor((x-30)/2);

        %Create figure and UI controls.
        hFig = figure('Name','3D Viewer', 'units', 'pixels', 'position',[(screenSize(3)-x)/2 (screenSize(4)-y)/2 x y], 'MenuBar', 'none', 'NumberTitle', 'off', 'resize', 'off', 'Tag', 'CERR_3DDoseViewer', 'color', [0 0 0]);

        hAxis = axes('position', [0 0 .7 1], 'color', [0 0 0]);

        hFrame = uicontrol('style', 'frame', 'units', 'normalized', 'position', [.7 .65 .28 .33]);

        %Create controls for structure selection
        structures = {'Select...' planC{indexS.structures}.structureName};
        nStructs = length(structures);
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.71 .90 .1 .05], 'string', 'Structure:');
        uicontrol('style', 'popupmenu', 'units', 'normalized', 'position', [.82 .90 .15 .05], 'string', structures, 'value', 1, 'callback', 'doseViewer3D(''struct_dropdown'')');

        ud.surfaceHandles = cell(nStructs,1);
        ud.surfaceVisible = zeros(nStructs,1);

        doses = {planC{indexS.dose}.fractionGroupID};
        %doses = {doses{:} 'DoseDiff', 'DTA'};
        uicontrol('style', 'text', 'units', 'normalized', 'position', [.71 .84 .1 .05], 'string', 'Dose:');
        uicontrol('style', 'popupmenu', 'units', 'normalized', 'position', [.82 .84 .15 .05], 'string', doses, 'value', 1, 'callback', 'doseViewer3D(''dose_dropdown'')');

        %Alpha Controls

        %Colorbar

        %Snap view:
        views = {'Transverse', 'Sagittal', 'Coronal'};
        uicontrol('style', 'pushbutton', 'units', 'normalized', 'position', [.735 .78 .1 .05], 'string', 'Trans', 'callback', 'doseViewer3D(''view_trans'')');
        uicontrol('style', 'pushbutton', 'units', 'normalized', 'position', [.735 .72 .1 .05], 'string', 'Sag', 'callback', 'doseViewer3D(''view_sag'')');
        uicontrol('style', 'pushbutton', 'units', 'normalized', 'position', [.845 .78 .1 .05], 'string', 'Cor', 'callback', 'doseViewer3D(''view_cor'')');

        %Coordinate info.
        ud.handle.xPos = uicontrol('style', 'text', 'units', 'normalized', 'position', [.71 .66 .09 .05], 'string', 'x=', 'horizontalAlignment', 'left');
        ud.handle.yPos = uicontrol('style', 'text', 'units', 'normalized', 'position', [.80 .66 .09 .05], 'string', 'y=', 'horizontalAlignment', 'left');
        ud.handle.zPos = uicontrol('style', 'text', 'units', 'normalized', 'position', [.89 .66 .09 .05], 'string', 'z=', 'horizontalAlignment', 'left');

        %         uicontrol('style', 'text', 'units', 'normalized', 'position', [.71 .78 .1 .05], 'string', 'View');
        %         uicontrol('style', 'popupmenu', 'units', 'normalized', 'position', [.82 .78 .15 .05], 'string', views, 'value', 1, 'callback', 'doseViewer3D(''view_dropdown'')');
        %
        %Default to doseSet 1.
        %doseSet = 1;

        [xV, yV, zV] = getDoseXYZVals(planC{indexS.dose}(doseSet));
        dA = getDoseArray(doseSet);
        doseBrowser3D('init', hAxis, dA, xV, yV, zV);

        doseBrowser3D('SETPLANEMOTIONCALLBACK', hAxis, 'doseViewer3D(''PLANE_MOVED'')');
        doseBrowser3D('SETPLANEMOTIONDONECALLBACK', hAxis, 'doseViewer3D(''PLANE_MOTION_DONE'')');

        ud.hAxis = hAxis;
        colormap(hAxis, CERRColorMap(stateS.optS.doseColormap));
        set(hFig, 'userdata', ud);

        doseViewer3D('PLANE_MOVED');

    case 'PLANE_MOVED'
        ud = get(hFig, 'userdata');
        hAxis = ud.hAxis;
        [x,y,z] = doseBrowser3D('GETPLANECOORDINATES', hAxis);
        set(ud.handle.zPos, 'String', ['z=' num2str(z, 3)]);
        set(ud.handle.xPos, 'String', ['x=' num2str(x, 3)]);
        set(ud.handle.yPos, 'String', ['y=' num2str(y, 3)]);

    case 'PLANE_MOTION_DONE'
%         ud = get(hFig, 'userdata');
%         hAxis = ud.hAxis;
%         [x,y,z] = doseBrowser3D('GETPLANECOORDINATES', hAxis);

%         setAxisInfo(stateS.handle.CERRAxis(1), 'coord', z);
%         setAxisInfo(stateS.handle.CERRAxis(4), 'coord', x);
%         setAxisInfo(stateS.handle.CERRAxis(6), 'coord', y);
%         plancheckCallback('refresh');

%         for i=1:length(stateS.handle.CERRAxis)
%             viewAx = getAxisInfo(stateS.handle.CERRAxis(i),'view');
%             if strcmpi('transverse',viewAx)
%                 setAxisInfo(stateS.handle.CERRAxis(i), 'coord', z);
%             elseif strcmpi('sagittal',viewAx)
%                 setAxisInfo(stateS.handle.CERRAxis(i), 'coord', x);
%             elseif strcmpi('coronal',viewAx)
%                 setAxisInfo(stateS.handle.CERRAxis(i), 'coord', y);
%             end
%         end
%
%         sliceCallback('refresh')

    case 'VIEW_DROPDOWN'
        ud = get(hFig, 'userdata');
        hAxis = ud.hAxis;

        viewNum = get(gcbo, 'value');
        switch viewNum
            case 1
                view(hAxis, 0, 90);
            case 2
                view(hAxis, -90, 0);
            case 3
                view(hAxis, 0, 0);
            case 4
                %Do nothing: floating view.
        end

    case 'VIEW_TRANS'
        ud = get(hFig, 'userdata');
        hAxis = ud.hAxis;
        view(hAxis, 0, 90);

    case 'VIEW_SAG'
        ud = get(hFig, 'userdata');
        hAxis = ud.hAxis;
        view(hAxis, -90, 0);

    case 'VIEW_COR'
        ud = get(hFig, 'userdata');
        hAxis = ud.hAxis;
        view(hAxis, 0, 0);

    case 'STRUCT_DROPDOWN'
        ud = get(hFig, 'userdata');
        hAxis = ud.hAxis;

        %Get selected struct.
        structNum = get(gcbo, 'value') - 1;

        if structNum == 0
            return;
        end

        %Toggle surface visible.
        ud.surfaceVisible(structNum) = xor(ud.surfaceVisible(structNum),1);

        %Draw surface if necessary.
        if isempty(ud.surfaceHandles{structNum})
            p = isoStruct(structNum);
            ud.surfaceHandles{structNum} = p;
        else
            p = ud.surfaceHandles{structNum};
        end

        %Add/remove star to denote struct is on/off
        structList = get(gcbo, 'string');
        switch ud.surfaceVisible(structNum);
            case 0
                structList{structNum+1} = [planC{indexS.structures}(structNum).structureName];
                set(p, 'visible', 'off');
            case 1
                structList{structNum+1} = ['*' planC{indexS.structures}(structNum).structureName];
                set(p, 'visible', 'on');
        end
        set(gcbo, 'string', structList);
        set(hFig, 'userdata', ud);

    case 'DOSE_DROPDOWN'
        ud = get(hFig, 'userdata');
        hAxis = ud.hAxis;

        doseNum = get(gcbo, 'value');
        strings = get(gcbo, 'string');

        switch strings{doseNum}
            case 'DTA'
                dS = [getAxisInfo(stateS.handle.CERRAxis(1), 'doseSets') ...
                    getAxisInfo(stateS.handle.CERRAxis(2), 'doseSets')];
                if length(dS) == 2
                    [xV1, yV1, zV1] = getDoseXYZVals(planC{indexS.dose}(dS(1)));
                    [xV, yV, zV] = getDoseXYZVals(planC{indexS.dose}(dS(2)));
                    dose3M1 = getDoseArray(dS(1));
                    dose3M2 = getDoseArray(dS(2));
                    [dA, ptc] = calculateGammaFcn1(dose3M1, dose3M2, [abs(yV1(2)-yV1(1)) abs(xV1(2)-xV1(1)) abs(zV1(2)-zV1(1))]);
                    colormap(hAxis, CERRColorMap('DTA'));
                else
                    return;
                end
            case 'DoseDiff'
                if doseNum == 0
                    dS = [];
                else
                    dS = [getAxisInfo(stateS.handle.CERRAxis(1), 'doseSets') ...
                    getAxisInfo(stateS.handle.CERRAxis(2), 'doseSets')];
                end
                [xV1, yV1, zV1] = getDoseXYZVals(planC{indexS.dose}(dS(1)));
                [xV, yV, zV] = getDoseXYZVals(planC{indexS.dose}(dS(2)));
                dose3M1 = getDoseArray(dS(1));
                dose3M2 = getDoseArray(dS(2));
                dA = dose3M1 - dose3M2;

            otherwise
                dA = getDoseArray(doseNum, planC);
                [xV, yV, zV] = getDoseXYZVals(planC{indexS.dose}(doseNum));
                colormap(hAxis, CERRColorMap(stateS.optS.doseColormap));
        end
        doseBrowser3D('NEWDOSE', hAxis, dA, xV, yV, zV);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function noiseWriteParamsToFID(noise, FID)

% NOISEWRITEPARAMSTOFID Write the noise parameters to a stream.
% FORMAT
% DESC writes noise parameters to  a file stream.
% ARG noise : the noise structure that is being written.
% ARG FID : the file ID of the stream that is being written.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2008
%
% SEEALSO : noiseWriteToFID

% NOISE

writeIntToFID(FID, 'outputDim', noise.numProcess);
writeIntToFID(FID, 'numParams', noise.nParams);
if strcmp(noise.type, 'ncnm')
  writeIntToFID(FID, 'gammaSplit', noise.gammaSplit);
end
fhandle = str2func([noise.type 'NoiseExtractParam']);
params = fhandle(noise);
doubleMatrixWriteToFID(params, FID);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function hEllipse = ellipsedraw(hand,a,b,x0,y0,phi,lineStyle)
%ELLIPSEDRAW can draw an arbitrary ellipse with given parameters.
%   The properties of that ellipse plot can be customized
%   by setting the ellipse handle.
%
%       hEllipse = ellipsedraw(hand,a,b,x0,y0,phi,lineStyle)
%
%   Input parameters:
%       hand        Parent handle for ellipse
%       a           Value of the major axis
%       b           Value of the minor axis
%       x0          Abscissa of the center point of the ellipse
%       y0          Ordinate of the center point of the ellipse
%       phi         Angle between x-axis and the major axis
%       lineStyle   Definition of the plotted line style
%
%   Output:
%       hEllipse    Handle of the ellipse
%
%   Simple usage:
%       ellipsedraw(5,3);
%       ellipsedraw(5,3,'g--');
%       ellipsedraw(5,3,pi/4);
%
%   Complete usage:
%       h = ellipsedraw(5,3,1,-2,pi/4,'r-.');
%       set(h,'LineWidth',2);

% Designed by: Lei Wang, <WangLeiBox@hotmail.com>, 25-Mar-2003.
% Last Revision: 01-Apr-2003.
% Dept. Mechanical & Aerospace Engineering, NC State University.
% Copyright (c)2003, Lei Wang <WangLeiBox@hotmail.com>
%$Revision: 1.0 $  $ 4/1/2003 5:42:24 PM $

if (nargin < 3)||(nargin > 7),
    error('Too few or too many arguments.');

elseif nargin == 3
    x0 = 0;     y0 = 0;
    phi = 0;    lineStyle = 'b-';

elseif nargin == 4
    if ischar(x0) == 1
        lineStyle = x0;
        x0 = 0; y0 = 0;
        phi = 0;
    else
        phi = x0;
        x0 = 0; y0 = 0;
        lineStyle = 'b-';
    end

elseif nargin == 5
    phi = 0;    lineStyle = 'b-';

elseif nargin == 6
    lineStyle = 'b-';
end

theta = [-0.03:0.01:2*pi];

% Parametric equation of the ellipse
%----------------------------------------
 x = a*cos(theta);
 y = b*sin(theta);

% Coordinate transform
%----------------------------------------
 X = cos(phi)*x - sin(phi)*y;
 Y = sin(phi)*x + cos(phi)*y;
 X = X + x0;
 Y = Y + y0;

% Plot the ellipse
%----------------------------------------
 hEllipse = plot(hand,X,Y,lineStyle);

 %axis equal;

author: jramshur

repo: HRVAS

save_path: github-repos/MATLAB/jramshur-HRVAS


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

R = 0.3:1:20;
Tw = 0.1:0.1:4.5;
aa = a;

%maepi = [ -116.5 34.25 1992.48 6 6 7.2 10  ];
M = zeros(length(R),length(Tw))*nan; i1 = 0; i2 = 0;
maepi(1,3) = 1992.48;
an = maepi(1,:);
y = an(1,2),x = an(1,1); z = an(:,7); Te = an(:,3);
%figure_w_normalized_uicontrolunits(map)
axes(hs);
[x,y] = ginput(1)
Te = max(a.Date)
z = 5

l = a.Date < Te;
a = a.subset(l);

dT = max(a.Date) - min(a.Date);
di = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2 + (a.Depth-z).^2) ;

N =[];
for i1 = 1:1:length(R)
    ;
    % take first ni points
    %
    l = di <= R(i1);
    b = a.subset(l);      % ne
    for i2 = 1:1:length(Tw)

        l = b(:,3) >= Te-Tw(i2) & b(:,3) < Te;
        Q = length(b(l,3)) ;   % This is the annual rate
        B = length(b(:,1))*Tw(i2)/dT;
        P = poisscdf(Q,B);

        M(i1,i2) = P;

    end
    N =[N ;  length(b)];
    i2 = 0;
end
i1 = 0;

% plot the results
plotqmatri;

a = aa;

author: CelsoReyes


text:
function [rc,fval,it] = ARjones(ti,xi,rcinit,p)

nobs = length(xi);

opties = optimset('Display','off','TolX',.001/sqrt(nobs),'TolFun',.0001);

[rc_tan,fval,exitflag,output]= fminunc('Jonesfit',tan(.5*pi*rcinit),opties,ti,xi,p);
rc = 2/pi*atan(rc_tan);
it = output.iterations;


text:
function [ event ] = read_ah5_markers(hdr, filename )

if isempty(hdr)
    orig = read_ahdf5_hdr(filename);
    hdr.orig = orig;
    hdr.Fs = orig.channels(1).samplingRate;
    hdr.nChans = numel(orig.channels);
    hdr.nSamples = orig.numberOfSamples;
    hdr.nTrials = orig.numberOfBlocks;
    hdr.nSamplesPre = 0;
    hdr.label = orig.label;
end

marks = h5read(filename, '/Blocks/1/markers');
labels = marks.label';
positions = marks.position';
durations = marks.duration';
values = marks.value';

for i=1:size(labels, 1)
    event(i).type = strcat(labels(i, 1:256));
    event(i).value = values(i);
    event(i).sample = positions(i) * hdr.Fs;
    event(i).duration = durations(i) * hdr.Fs;
end
end


text:
function [result, ok] = import_ndk(nFunction, sFilename)
    % import NDK formatted catalog Global Centroid-Moment-Tensor (CMT) catalog,
    % formarly the Harvard CMT catalog
    %
    % see also NDK

    ok=false;
    if nFunction == FilterOp.getDescription
        % respond with brief summary of function, used in selection box
        result = 'Centroid-Moment-Tensor (CMT) Catalog, NDK file (formerly Harvard CMT catalog)';
    elseif nFunction == FilterOp.getWebpage
        % respond with an html file that should help
        % website last accessed 18-July-2018
        result = 'https://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/allorder.ndk_explained';
    else
        ndk = NDK.read(sFilename);
        result = ndk.toZmapCatalog();
        ok=true;
    end
end


text:
function rootPath=mraRootPath()
% Determine path to root of the mrAnatomy directory
%
%        rootPath = mraRootPath;
%
% This function MUST reside in the directory at the base of the mrAnatomy
% directory structure
%
% VISTASOFT, Stanford, Wandell

rootPath=which('mraRootPath');

[rootPath,fName,ext]=fileparts(rootPath);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [Data, para, vocab] = LoadWavLabel_CHiME4(para, step, dataset)
switch lower(dataset)
    case {'dt05','et05'}
        wavlist = ['../Kaldi/data/' dataset '_multi_noisy/wav.scp'];
        ali_file = [para.local.aliDir '_dt05/ali.txt'];
    case 'tr05'
        wavlist = ['../Kaldi/data/' dataset '_multi_noisy/wav.scp'];
        ali_file = [para.local.aliDir '/ali.txt'];
end

[ali, vocab] = LoadKaldiFrameLabel(ali_file);
wavlist = my_cat(wavlist);
wavlist = wavlist(step:step:end);

wavreader.name = 'wavfile';
wavreader.array = 1;
wavreader.multiArrayFiles = 1;

nCh = 6;
fs = 16000;
frame_size = fs*0.025;
frame_shift = fs*0.01;

% Load data file list
wav_noisy = {};  label = {};
for si = 1:length(wavlist)
    words = ExtractWordsFromString_v2(wavlist{si});
    curr_uttID = words{1};
    wavfile = words{2};
    PrintProgress(si, length(wavlist), 100, curr_uttID);

    % get alignment
    if ~strcmpi(para.local.data, 'mixed') && isempty(regexp(lower(curr_uttID), para.local.data))
        continue;
    end

    words = ExtractWordsFromString_v2(curr_uttID, '_');
    curr_uttID2 = [words{2} '_' words{3}(1:3) '_' words{4}];
    fieldname = ['U_' curr_uttID2];
    if ~isfield(ali, fieldname); continue; end    % if there is no label for current utterance, skip it.
    curr_label = ali.(fieldname);
    curr_label = curr_label(:)'+1;  % convert to row vector and the index starts with 1

    % read in the waveform
    % for the two channel track, we may choose a random pair, or use all
    % the channel pairs.
    if length(para.topology.useChannel)==1
        if para.topology.useChannel==6
            ch_idx = randperm(6);
        elseif para.topology.useChannel==5
            ch_idx = randperm(6);
            ch_idx(ch_idx==2) = [];     % we don't use channel 2
        elseif strcmpi(para.local.pair, 'randPair')
            ch_idx = randperm(nCh);
            ch_idx(ch_idx==2) = [];     % we don't use channel 2
            ch_idx = sort(ch_idx(1:para.topology.useChannel));
        elseif strcmpi(para.local.pair, 'allPair')
            ch_idx = [1 3; 1 4; 1 5; 1 6; 3 4; 3 5; 3 6; 4 5; 4 6; 5 6];
        end
    else
        ch_idx = para.topology.useChannel;
    end
    words = ExtractWordsFromString_v2(wavfile, '/');
    wavfileRoot = [para.local.wavroot_noisy '/' words{end-1} '/' words{end}(1:end-5)];

    for pi = 1:size(ch_idx,1)   % add all channel pairs to the training data
        clear wavfileArray
        for i=1:size(ch_idx,2)
            wavfileArray{i} = [wavfileRoot num2str(ch_idx(pi,i)) '.wav'];
        end
        [wav] = InputReader(wavfileArray, wavreader);
        wav = StoreWavInt16(wav);

        % synchronize the length of label and wav
        nFr_feat = enframe_decide_frame_number(size(wav,2), frame_size, frame_shift);
        nFr_label = length(curr_label);
        if nFr_feat>=nFr_label
            requiredLen = DecideWavLen4XFrames(nFr_label, frame_size, frame_shift);
            wav(:,requiredLen+1:end) = [];
        elseif nFr_feat<nFr_label
            curr_label = curr_label(1:nFr_feat);
        end

        if para.local.useFileName
            for i=1:size(ch_idx,2)
                wavfileArray{i} = sprintf('%s 0 %2.3f', wavfileArray{i}, size(wav,2)/fs);
            end
            wav_noisy{end+1} = wavfileArray;
        else
            wav_noisy{end+1} = wav;
        end
        label{end+1} = curr_label;
    end
end

Data(1).data = wav_noisy;
Data(2).data = label;

para.IO.inputFeature = [1 1];
para.IO.DataSyncSet{1} = [];
para.IO.frame_rate = [16000 100];
para.IO.isTensor = [1 1];
if para.local.useFileName
    para.IO.inputFeature(1) = 0;
    wavreader.precision = 'int16';
    para.IO.fileReader(1) = wavreader;
    para.IO.fileReader(2).name = '';
end
end

author: singaxiong

repo: SignalGraph


text:
function r = mtimes(p,q)
%TIMES this function multiplies the x, y, and z values of the structures that I use in
% the FVtool.
% reassign the scalar to new variable a and the struct to new variable b
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'FaceVariable')&&isa(q, 'FaceVariable'))
    error('FVMtool: Wrong use of mtimes for a face variable. Try using .* instead.');
elseif isa(p, 'FaceVariable')
    r=p;
    r.xvalue = p.xvalue*q;
    r.yvalue = p.yvalue*q;
    r.zvalue = p.zvalue*q;
else
    r=q;
    r.xvalue = p*q.xvalue;
    r.yvalue = p*q.yvalue;
    r.zvalue = p*q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function gx = p01_gx ( x )

%*****************************************************************************80
%
%% P01_GX evaluates the underlying function for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(2), the point at which the function is to
%    be evaluated.
%
%    Output, real GX(2), the value of the function at X.
%
  gx(1) = x(1) - ( ( x(2) - 5.0 ) * x(2) + 2.0 ) * x(2) - 13.0;
  gx(2) = x(1) + ( ( x(2) + 1.0 ) * x(2) - 14.0 ) * x(2) - 29.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% The  command invfreqs

%verification of the result of the command freqs
num=[ -2/j^2 7 0];
den=[3 0 2];
w=0:.1:10;
H=freqs(num,den,w);
freqs(num,den,w)
[num1,den1]=invfreqs(H,w,2,2)

%Equivalent frequency responses
[num2,den2]=invfreqs(H,w,3,3)
freqs(num2,den2,w)

figure
[num3,den3]=invfreqs(H,w,10,7);
freqs(num3,den3,w)


text:
% Demo of Hessenberg QR to RQ
clc;
close all;
clearvars;
n = 4;
% The binomial causes power method to oscillate
% A = gallery('binomial',n);
A = gallery('lotkin',n);

H = hess(A)
H1 = spx.la.hessenberg.qr_rq(H);

[Q, R] = qr(H);
H2 = R * Q;

H1
H2
fprintf('max diff ::: %f\n', max(max(abs(abs(H1)-abs(H2)))));

author: indigits

repo: sparse-plex