text:
% TESTFITPLANE - demonstrates RANSAC plane fitting
%
% Usage: testfitplane(outliers, sigma, t, feedback)
%
% Arguments:
%               outliers - Fraction specifying how many points are to be
%                          outliers.
%               sigma    - Standard deviation of inlying points from the
%                          true plane.
%               t        - Distance threshold to be used by the RANSAC
%                          algorithm for deciding whether a point is an
%                          inlier.
%               feedback - Optional flag 0 or 1 to turn on RANSAC feedback
%                          information.
%
%  Try using:  testfitplane(0.3, 0.05, 0.05)
%
% See also: RANSACFITPLANE, FITPLANE

% Copyright (c) 2003-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% June 2003

function testfitplane(outliers, sigma, t, feedback)
    if nargin == 3
	feedback = 0;
    end

    % Hard wire some constants - vary these as you wish

    npts = 100;  % Number of 3D data points

    % Define a plane  ax + by + cz + d = 0
    a = 10; b = -3; c = 5; d = 1;

    B = [a b c d]';
    B = B/norm(B);

    outsigma = 30*sigma;  % outlying points have a distribution that is
                          % 30 times as spread as the inlying points

    vpts = round((1-outliers)*npts);  % No of valid points
    opts = npts - vpts;               % No of outlying points

    % Generate npts points in the plane
    X = rand(1,npts);
    Y = rand(1,npts);
    Z = (-a*X -b*Y -d)/c;

    XYZ =  [X
	    Y
	    Z];

    % Add uniform noise of +/-sigma
    XYZ = XYZ + (2*rand(size(XYZ))-1)*sigma;

    % Generate opts random outliers

    n = length(XYZ);
    ind = randperm(n);  % get a random set of point indices
    ind = ind(1:opts);  % ... of length opts

    % Add uniform noise of outsigma to the points chosen to be outliers.
%    XYZ(:,ind) = XYZ(:,ind)  + (2*rand(3,opts)-1)*outsigma;

    XYZ(:,ind) = XYZ(:,ind)  +   sign(rand(3,opts)-.5).*(rand(3,opts)+1)*outsigma;

    % Display the cloud of points
    figure(1), clf, plot3(XYZ(1,:),XYZ(2,:),XYZ(3,:), 'r*');

    % Perform RANSAC fitting of the plane
    [Bfitted, P, inliers] = ransacfitplane(XYZ, t, feedback);

    fprintf('Original plane coefficients: ');
    fprintf('%8.3f ',B);
    fprintf('\nFitted plane coefficients:   ');
    fprintf('%8.3f ',Bfitted);
    fprintf('\n');

    % Display the triangular patch formed by the 3 points that gave the
    % plane of maximum consensus
    patch(P(1,:), P(2,:), P(3,:), 'g')

    box('on'), grid('on'), rotate3d('on')

    fprintf('\nRotate image so that planar patch is seen edge on\n');
    fprintf('If the fit has been successful the inlying points should\n');
    fprintf('form a line\n\n');

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function varargout = dartel3(varargin)
% DARTEL 3D image registration stuff
%
%_______________________________________________________________________
%
% FORMAT v = dartel3(v,g,f,param)
% v     - flow field n1*n2*n3*3 (single precision float)
% g     - first image n1*n2*n3*n4 (single precision float)
% f     - second image n1*n2*n3*n4 (single precision float)
% param - 9 parameters (settings)
%         - [1] Regularisation type, can take values of
%           - 0 Linear elasticity
%           - 1 Membrane energy
%           - 2 Bending energy
%         - [2][3][4] Regularisation parameters
%           - For "membrane energy", the parameters are
%             lambda, unused and id.
%           - For "linear elasticity", the parameters are
%             mu, lambda, and id
%           - For "bending energy", the parameters are
%             lambda, id1 and id2, such that regularisation is by
%             (-lambda*\grad^2 + id1)^2 + id2
%         - [5] Levenberg-Marquardt regularisation
%         - [6] Number of Full Multigrid cycles
%         - [7] Number of relaxation iterations per cycle
%         - [8] K, such that 2^K time points are used to
%               generate the deformations.  A value of zero
%               indicates a small deformation model.
%         - [9] code of 0, 1 or 2.
%               0 - asymmetric sums of squares objective function.
%               1 -  symmetric sums of squares objective function.
%               2 - assumes multinomial distribution, where template
%                   encodes the means and interpolation of temlate
%                   done using logs and softmax function.
%
% This is for performing a single iteration of the Dartel optimisation.
% All flow fields and images are represented by single precision floating
% point values. Images can be scalar fields, in which case the objective
% function is the sum of squares difference.  Alternatively, images can be
% vector fields, in which case the objective function is the sum of squares
% difference between each scalar field + the sum of squares difference
% between one minus the sum of the scalar fields.
%
%_______________________________________________________________________
%
% FORMAT v = dartel3('cgs',A, b, param)
% v     - the solution
% A     - parameterisation of 2nd derivatives
% b     - parameterisation of first derivatives
% param - 6 parameters (settings)
%         - [1] Regularisation type, can take values of
%           - 0 Linear elasticity
%           - 1 Membrane energy
%           - 2 Bending energy
%         - [2][3][4] Voxel sizes
%         - [5][6][7] Regularisation parameters
%           - For "membrane energy", the parameters are
%             lambda, unused and id.
%           - For "linear elasticity", the parameters are
%             mu, lambda, and id
%           - For "bending energy", the parameters are
%             lambda, id1 and id2.
%         - [8] Tolerance.  Indicates required degree of accuracy.
%         - [9] Maximum number of iterations.
%
% This is for solving a set of equations using a conjugate gradient
% solver. This method is less efficient than the Full Multigrid.
% v = inv(A+H)*b
% A, b and v are all single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT v = dartel3('fmg',A, b, param)
% v     - the solution n1*n2*n3*3
% A     - parameterisation of 2nd derivatives
% b     - parameterisation of first derivatives
% param - 6 parameters (settings)
%         - [1] Regularisation type, can take values of
%           - 0 Linear elasticity
%           - 1 Membrane energy
%           - 2 Bending energy
%         - [2][3][4] Voxel sizes
%         - [5][6][7] Regularisation parameters
%           - For "membrane energy", the parameters are
%             lambda, unused and id.
%           - For "linear elasticity", the parameters are
%             mu, lambda, and id
%           - For "bending energy", the parameters are
%             lambda, id1 and id2.
%         - [8] Number of Full Multigrid cycles
%         - [9] Number of relaxation iterations per cycle
%
% Solve equations using a Full Multigrid method.  See Press et al
% for more information.
% v = inv(A+H)*b
% A, b and v are all single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT [y,J] = dartel3('Exp', v, param)
% v - flow field
% J - Jacobian. Usually a tensor field of Jacobian matrices, but can
%     be a field of Jacobian determinants.
% param - 2 (or 3) parameters.
%         [1] K, the number of recursions (squaring steps), such
%             that exponentiation is done using an Euler-like
%             integration with 2^K time steps.
%         [2] a scaling parameter.
%         If there is a third parameter, and it is set to 1, then
%         the J will be the Jacobian determinants.
%
% A flow field is "exponentiated" to generate a deformation field
% using a scaling and squaring approach.  See the work of Arsigny
% et al, or Cleve Moler's "19 Dubious Ways" papers.
%
%_______________________________________________________________________
%
% FORMAT m = dartel3('vel2mom', v, param)
% v     - velocity (flow) field n1*n2*n3*3.
% param - 4 parameters (settings)
%         - [1] Regularisation type, can take values of
%           - 0 Linear elasticity
%           - 1 Membrane energy
%           - 2 Bending energy
%         - [2][3][4] Voxel sizes
%         - [5][6][7] Regularisation parameters
%           - For "membrane energy", the parameters are
%             lambda, unused and id.
%           - For "linear elasticity", the parameters are
%             mu, lambda, and id
%           - For "bending energy", the parameters are
%             lambda, id1 and id2.
% m       - `momentum' field n1*n2*n3*3.
%
% Convert a flow field to a momentum field by m = H*v, where
% H is the large sparse matrix encoding some form of regularisation.
% v and m are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT y3 = dartel3('comp',y1,y2)
% y1, y2 - deformation fields n1*n2*n3*3.
% y3     - deformation field field n1*n2*n3*3.
%
% Composition of two deformations y3 = y1(y2)
% y1, y2 and y3 are single precision floating point.
%
%
%
% FORMAT [y3,J3] = dartel3('comp', y1, y2, J1, J2)
% y1, y2 - deformation fields n1*n2*n3*3.
% y3     - deformation field n1*n2*n3*3.
% J1, J2 - Jacobian tensor fields n1*n2*n3*3*3.
% J3     - Jacobian tensor field n1*n2*n3*3*3.
%
% Composition of two deformations, with their Jacobian fields.
% All fields are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT f2 = dartel3('samp', f1, y)
% f1 - input image(s) n1*n2*n3*n4
% y  - points to sample n1*n2*n3*3
% f2 - output image n1*n2*n3*3
%
% Sample a function according to a deformation.
% f2 = f1(y)
% f1, f2 and y are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT v2 = dartel3('resize', v1, dim)
% v1  - input fields n1*n2*n3*n4
% v2  - output field dim1*dim2*dim3*n4
% dim - output dimensions
%
% Resize a field according to dimensions dim.  This is
% a component of the FMG approach.
%
%_______________________________________________________________________
%
% FORMAT v3 = dartel3('brc', v1, v2)
% v1, v2, v3 - flow fields n1*n2*n3*3
%
% Lie Bracket.  Useful for many things
% e.g. Baker-Campbell-Haussdorf series expansion.
% The Lie bracket is denoted by
% v3 = [v1,v2]
% and on scalar fields, is computed by
% v3 = J1*v2 - J2*v1, where J1 and J2 are the Jacobian
% tensor fields. For matrices, the Lie bracket is simply
% [A,B] = A*B-B*A
%
%_______________________________________________________________________
%
% Note that the boundary conditions are circulant throughout.
% Interpolation is trilinear, except for the resize function
% which uses a 2nd degree B-spline (without first deconvolving).
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: dartel3.m 5506 2013-05-14 17:13:43Z john $

%error('Not compiled for %s in MATLAB %s  (see make.m)\n', computer, version);

if nargin==0,
    error('Incorrect usage.');
elseif ~isa(varargin{1},'char'),
    param0 = varargin{4};
    switch param0(1),
    case 1,
        param = [param0(4) param0(2) 0  0 0];
    case 2,
        param = [param0(4) param0(3) param0(2) 0 0];
    otherwise
        param = [param0(4) 0 0 param0(2) param0(3)];
    end
    param = [param param0(5:end)];
    [varargout{1:nargout}] = spm_diffeo('dartel',varargin{1},varargin{2:3},param);
else
    switch varargin{1},
    case {'fmg','cgs'},
        param0 = varargin{4};
        vx2    = mean(param0(2:4).^2);
        switch param0(1),
        case 1,
            param = [param0(7) 0 0 param0(5) 0  0 0]*vx2;
        case 2,
            param = [param0(7) param0(6) param0(5) 0 0]*vx2;
        otherwise
            param = [param0(7) 0 0 param0(5) param0(6)]*vx2;
        end
        param = [param9(2:4) param param0(8:end)];
        [varargout{1:nargout}] = spm_diffeo('dartel',varargin{1:3},param);

    case {'vel2mom'},
        param0 = varargin{3};
        vx2    = mean(param0(2:4).^2);
        switch param0(1),
        case 1,
            param = [param0(7) param0(5) 0  0 0]*vx2;
        case 2,
            param = [param0(7) param0(6) param0(5) 0 0]*vx2;
        otherwise
            param = [param0(7) 0 0 param0(5) param0(6)]*vx2;
        end
        param = [param0(2:4) param param0(8:end)];
        [varargout{1:nargout}] = spm_diffeo(varargin{1:2},param);

    otherwise
        [varargout{1:nargout}] = spm_diffeo(varargin{:});
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function spm_check_registration(varargin)
% A visual check of image registration quality
% FORMAT spm_check_registration
% FORMAT spm_check_registration(images)
% Orthogonal views of one or more images are displayed. Clicking in
% any image moves the centre of the orthogonal views. Images are
% shown in orientations relative to that of the first selected image.
% The first specified image is shown at the top-left, and the last at
% the bottom right. The fastest increment is in the left-to-right
% direction (the same as you are reading this).
%__________________________________________________________________________
% Copyright (C) 1997-2019 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_check_registration.m 7759 2019-12-19 11:50:50Z guillaume $

SVNid = '$Rev: 7759 $';

%-Get input
%--------------------------------------------------------------------------
if nargin
    if nargin > 1 && iscellstr(varargin{2})
        images = varargin{1};
    elseif isstruct(varargin{1})
        images = [varargin{:}];
    else
        images = char(varargin);
    end
else
    [images, sts] = spm_select([1 24],'image','Select images');
    if ~sts, return; end
end

if ischar(images), images = spm_vol(images); end
if numel(images) > 24
    if ~isdeployed, addpath(fullfile(spm('Dir'),'spm_orthviews')); end
    img = cell(1,numel(images));
    for i=1:numel(images)
        img{i} = [images(i).fname ',' num2str(images(i).n(1))];
    end
    spm_ov_browser('ui',char(img));
    return
end
images = images(1:min(numel(images),24));

%-Print
%--------------------------------------------------------------------------
spm('FnBanner',mfilename,SVNid);                                        %-#
exactfname  = @(f) [f.fname ',' num2str(f.n(1))];
cmddispone  = 'spm_image(''display'',''%s'')';
cmddispall  = 'spm_check_registration(''%s'')';
if spm_platform('desktop')
    str     = '';
    for i=1:numel(images)
        str = [str sprintf('''%s'',',exactfname(images(i)))];
    end
    dispall = [' (' spm_file('all','link',sprintf(cmddispall,str(2:end-2))) ')  '];
else
    dispall = '        ';
end
for i=1:numel(images)
    if i==1,     fprintf('Display ');                                   %-#
    elseif i==2, fprintf('%s',dispall);
    else         fprintf('        '); end
    fprintf('%s\n',spm_file(exactfname(images(i)),'link',cmddispone));  %-#
end

%-Display
%--------------------------------------------------------------------------
fg = spm_figure('GetWin','Graphics');
spm_figure('Clear',fg);
spm_orthviews('Reset');
dm = get(fg,'Position');
mn = length(images);

minwasted = Inf;
m = 0; n = 0;
for m1=1:mn
    n1 = ceil(mn/m1);
    s  = max(dm(4)/m1,dm(3)/n1);
    wasted = (s*m1-dm(4))*dm(3) + (s*n1-dm(3))*dm(4) + (m1*n1-mn)*s^2*1.01;
    if wasted < minwasted
        minwasted = wasted;
        m = m1;
        n = n1;
    end
end

w  = 1/n;
h  = 1/m;
ds = (w+h)*0.02;
for ij=1:mn
    i = 1-h*(floor((ij-1)/n)+1);
    j = w*rem(ij-1,n);
    handle = spm_orthviews('Image', images(ij),...
        [j+ds/2 i+ds/2 w-ds h-ds]);
    if ij==1, spm_orthviews('Space'); end
    spm_orthviews('AddContext',handle);
end

%-Backward compatibility with spm_check_registration(images,captions)
%--------------------------------------------------------------------------
if nargin > 1 && iscellstr(varargin{2})
    spm_orthviews('Caption',varargin{2}, varargin{3:end});
    %for ij=1:numel(varargin{2})
    %    spm_orthviews('Caption', ij, varargin{2}{ij}, varargin{3:end});
    %end
end

author: spm


text:
function [ x, seed ] = sphere01_sample ( n, seed )

%*****************************************************************************80
%
%% SPHERE01_SAMPLE picks random points on the unit sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of samples.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real X(3,N), the sample points.
%
  x = zeros ( 3, n );

  for j = 1 : n
%
%  Pick a uniformly random VDOT, which must be between -1 and 1.
%  This represents the dot product of the random vector with the Z unit vector.
%
%  Note: this works because the surface area of the sphere between
%  Z and Z + dZ is independent of Z.  So choosing Z uniformly chooses
%  a patch of area uniformly.
%
    [ vdot, seed ] = r8_uniform_01 ( seed );
    vdot = 2.0 * vdot - 1.0;

    phi = r8_acos ( vdot );
%
%  Pick a uniformly random rotation between 0 and 2 Pi around the
%  axis of the Z vector.
%
    [ theta, seed ] = r8_uniform_01 ( seed );
    theta = 2.0 * pi * theta;

    x(1,j) = cos ( theta ) * sin ( phi );
    x(2,j) = sin ( theta ) * sin ( phi );
    x(3,j) = cos ( phi );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function feats = make_rank_feats(input)

feats = input;
feats(:, [17 18]) = 1./(1+exp(-input(:, [17 18]))); % These are logistic predictions
feats = [feats, feats.^2, ones(size(feats,1), 1)];

author: Cloud-CV

repo: object-proposals


text:
function [model, metsUnmapped, rxnsUpdated] = mediaConstraints(model, uptakeRates, uptakeInChi, uptakeNames)
    %this function finds the highest uptake rate for each metabolite in the
    %uptakeRatesTable and sets this value as a lower boudary in the
    %exchange reaction for that metabolite. If maximum uptake rate is higer
    %than 0 (metabolite is only secreted by the cells) then lower boundary
    %is set to 0
    %
    % Inputs:
    %   model - Cobra model
    %   uptakesRatesTable - Table consisting of 'Chemical_Name', 'InChICode',
    %                       and numeric columns with the experimental data
    %                       (uptake/secretion rates in mmol/gDW/h )
    % Output:
    %   model -  Cobra model with constrained uptake rates
    %
    % 20190617 Agnieszka Wegrzyn

    model_temp = model;
    mediaInformation = uptakeInChi;
    EXrxns_all = model_temp.rxns(findExcRxns(model_temp));
    EXrxns = EXrxns_all(contains(EXrxns_all, 'EX_'));
    EXmets = findMetsFromRxns(model_temp, model_temp.rxns(findExcRxns(model_temp)));
    Rates = uptakeRates;
    rxnsUpdated = [];
    metsUnmapped = [];
    for i=1:length(mediaInformation)
        mets_temp = model_temp.mets(ismember(model_temp.metInChIString, mediaInformation(i)));

        if ~isempty(mets_temp)
            rxn_temp = EXrxns(ismember(EXrxns, findRxnsFromMets(model_temp,mets_temp(ismember(mets_temp, EXmets)))));
            model_temp = changeRxnBounds(model_temp, rxn_temp, min(min(Rates(i,:)),0)*1e4,'l');
            rxnsUpdated = [rxnsUpdated; rxn_temp];
        else
            metsUnmapped = [metsUnmapped; uptakeNames(i)];
        end
    end

model = model_temp;

end

author: opencobra

repo: cobratoolbox


text:
function W = compute_point_weight(pts, type, rings, options)

% compute_point_weight - compute a weight matrix
%
%   W = compute_point_weight(pts, type, rings, options);
%
%   W is sparse weight matrix and W(i,j)=0 is vertex i and vertex j are not
%   connected.
%
%   type is either
%       'combinatorial': W(i,j)=1 is vertex i is conntected to vertex j.
%       'distance': W(i,j) = 1/d_ij^2 where d_ij is distance between vertex
%           i and j.
%       'spring': W(i,j) = 1/d_ij where d_ij is distance between vertex
%           i and j.
%       'conformal' or 'dcp': W(i,j) = cot(alpha_ij)+cot(beta_ij) where alpha_ij and
%           beta_ij are the adjacent angle to edge (i,j). Refer to Skeleton
%           Extraction by Mesh Extraction_08, and Intrinsic Parameterizations of Surface Meshes_02.
%       'Laplace-Beltrami': W(i,j) = (cot(alpha_ij)+cot(beta_ij))/2 where alpha_ij and
%           beta_ij are the adjacent angle to edge (i,j). Refer to Computing discrete minimal surfaces
%           andtheir conjugates_93, Lemma 2 of On the convergence of metric and geometric properties of
%           polyhedral surfaces_06, and Characterizing Shape Using
%           Conformal Factors_08, and ,
%       'mvc': W(i,j) = [tan(/_kij/2)+tan(/_jil/2)]/d_ij where /_kij and /_jil
%           are angles at i
%
%   If options.normalize=1, the the rows of W are normalize to sum to 1.
%   If M.ring is offered, we can avoid compute it.
%
%   NB: we adapt it for better numeric stablity for contracting a model by adding the following lines
%     tmp = sum(W(i,:));
%     if tmp>10000
%         W(i,:) = W(i,:)*10000/tmp;
%     end
%   Copyright (c) 2009 jjcao
options.null = 0;

if nargin<2
    type = 'conformal';
end

switch lower(type)
    case 'combinatorial'
        W = compute_point_weight_combinatorial(pts, rings);
    case 'distance'
       warning('not implemented!');
    case 'spring'
        W = compute_point_weight_spring(pts, rings);
    case {'conformal','dcp'} % conformal laplacian
        W = compute_point_weight_dcp(pts, rings);
    case 'laplace-beltrami' %
        W = compute_point_weight_dcp(pts, rings)*0.5;
    case 'mvc'% mvc laplacian
        W = compute_point_weight_mvc(pts, rings);
    otherwise
        error('Unknown type!!')
end

%#########################################################################
function W = compute_point_weight_combinatorial(points, rings)
n = length(points);
W = sparse(n,n);
for i = 1:n
    ring = rings{i};
    if ring(1) == ring(end)
        ring = ring(1,1:(end-1));
    end
    for j = ring
        W(i,j) = 1.0;
    end
end
function W = compute_point_weight_spring(points, rings)
n = length(points);
W = sparse(n,n);
for i = 1:n
    vi = points(i,:);
    ring = rings{i};
    if ring(1) == ring(end)
        ring = ring(1,1:(end-1));
    end
    for j = ring
        vj = points(j,:);
        W(i,j) = 1./sqrt(sum((vi-vj).^2));
    end
    tmp = sum(W(i,:));
    if tmp>10000
        W(i,:) = W(i,:)*10000/tmp;
    end
end
%#########################################################################
function W = compute_point_weight_dcp(points, rings)
n = length(points);
W = sparse(n,n);
for i = 1:n
    ring = rings{i};

    tmp = size(ring,2)-1;
    for ii = 1: tmp
        j = ring(ii); k = ring(ii+1);
        vi = points(i,:);
        vj = points(j,:);
        vk = points(k,:);

        % new % Oscar08 use this
        u = vk-vi; v = vk-vj;
        cot1 = dot(u,v)/norm(cross(u,v));
        W(i,j) = W(i,j) + cot1;
        u = vj-vi; v = vj-vk;
        cot2 = dot(u,v)/norm(cross(u,v));
        W(i,k) = W(i,k) + cot2;
%         % old
%         % angles
%         alpha = myangle(vk-vi,vk-vj);
%         beta = myangle(vj-vi,vj-vk);
%         % add weight
%         W(i,j) = W(i,j) + cot( alpha );
%         W(i,k) = W(i,k) + cot( beta );
    end

    tmp = abs(sum(W(i,:)));
    if tmp>10000
        W(i,:) = W(i,:)*10000/tmp;
    end
end
function W = compute_point_weight_mvc(points, rings)
n = length(points);
W = sparse(n,n);
for i = 1:n
    ring = rings{i};

    tmp = size(ring,2)-1;
    for ii = 1: tmp
        j = ring(ii); k = ring(ii+1);
        vi = points(i,:);
        vj = points(j,:);
        vk = points(k,:);

        % angles
        alpha = myangle(vi-vk,vi-vj);
        % add weight
        W(i,j) = W(i,j) + tan( 0.5*alpha )/sqrt(sum((vi-vj).^2));
        W(i,k) = W(i,k) + tan( 0.5*alpha )/sqrt(sum((vi-vk).^2));
    end

    tmp = sum(W(i,:));
    if tmp>10000
        W(i,:) = W(i,:)*10000/tmp;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function beta = myangle(u,v);

du = sqrt( sum(u.^2) );
dv = sqrt( sum(v.^2) );
du = max(du,eps); dv = max(dv,eps);
beta = acos( sum(u.*v) / (du*dv) );

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function [B,EV,ED] = biharmonic_embedding(V,F,dim,p);
  % [B,EV,ED] = biharmonic_embedding(V,F);
  %
  % Takes a mesh (V,F) and returns an embedding using the spectrum of the
  % biharmonic operator. Then the biharmonic distance between two points i and
  % j can be computed as the euclidean distance between B(i,:) and B(j,:),
  % namely:
  %   dist_ij = sqrt(sum((B(i,:)-B(j,:)).^2,2));
  %
  % Input:
  %   V  vertex list
  %   F  face list
  %   dim requested dimension of the embedding
  %   Optional:
  %     p  exponent above eigen values
  %       0.5  "semi-harmonic" embedding
  %       1  commute time embedding, "harmonic"
  %       2  biharmonic {default}
  %       3  "triharmonic" embedding
  % Output:
  %   B  biharmonic embedding
  %   EV  eigenvectors used in embedding
  %   ED  eigenvalues used in embedding
  %

  % if dimension is not specfied use 4
  if(~exist('dim','var'))
    dim = 4;
  end

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

  % get cotangent matrix
  L = cotmatrix(V,F);
  % This should be better, but yaron seemed to use barycentric
  %M = massmatrix(V,F);
  M = massmatrix(V,F,'barycentric');

  % get dim+1 smallest magnitude eigenvalues and corresponding vectors
  [EV,ED] = eigs(L,M,dim+1,'sm');
  EV = EV(:, 2:end);
  ED = ED(2:end, 2:end);

  % This is not exactly the same, essentially it removes the mass matrix and
  % multiplies everything by a factos of -2
  % % This also works, because of the sign change in the eigenvalues matlab
  % % reverses the output order so 0.0 is the last eigenvalue
  % [EV,ED] = eigs(-2*L,M./sum(M(:)),dim+1,'sm');
  % EV = EV(:, 1:end-1);
  % ED = ED(1:end-1, 1:end-1);

  %  divide each eigenvector by corresponding eigenvalue
  %  divide the power by 2 first because it will appear in the denominator of
  %  distance computation *outside* the squared difference see (4) and (11)
  B = EV * (inv(abs(ED))^(p/2));
end

author: alecjacobson

repo: gptoolbox


text:
function F=FTaylor(deltaT,xCur,curT,a,dadx,d2adx,method)
%%FTAYLOR Simulate a nonlinear continuous-time random process specified by
%         the Langevin equation forward in time by a step-size of deltaT
%         using a Taylor scheme.
%
%INPUTS: deltaT The size of the single step over which to generate the
%               state transition matrix.
%          xCur The initial target state at time curT.
%          curT The time of the initial state xCur.
%             a The drift function in the continuous-time stochastic
%               dynamic model. It takes the state and a time variable as
%               its arguments, a(x,t).
%          dadx A xDimXxDim matrix of the derivative of the drift function
%               with respect to the state at state xCur and time curT.
%               This can also be a function that takes the state and a time
%               variable as its arguments, dadx(x,t).
%         d2adx A xDimXxDimXxDim matrix of the second derivative of the
%               drift function with respect to the state at state xCur and
%               time curT. The value at point (m,k,l) represents
%               d2a(m)/dx(k)dx(l). This can also be a function that takes
%               the state and a time variable as its arguments,
%               d2adx(x,t). If not provided, this is assumed to be zero.
%        method Set to 0 for the shorter Euler-Maruyama expansion,
%               otherwise will use Taylor expansion (default).
%
%OUTPUTS: F The state transition matrix under a nonlinear continuous-time
%           random process specified by the Langevin equation forward in
%           time by a step-size of deltaT using a strong Taylor scheme.
%
%The state prediction matrix is derived from the stochastic order 1.5
%Taylor scheme described in 10.4 of [1].
%
%REFERENCES:
%[1] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic
%    Differential Equations. Berlin: Springer, 1999.
%
%April 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

aCur=a(xCur,curT);
xDim=size(aCur,1);

if(nargin<6 || isempty(d2adx))
    d2adx=zeros(xDim,xDim,xDim);
end
if(nargin<7||isempty(method))
    method=1;
end
if(isa(dadx,'function_handle'))
    dadx=dadx(xCur,curT);
end
if(isa(d2adx,'function_handle'))
    d2adx=d2adx(xCur,curT);
end

if(method==0)
    F=eye(xDim)+deltaT*dadx;
else
    %assumes 3rd derivative is zero and all cross-derivatives are zero
    term1=zeros(xDim);
    for n=1:xDim
        term1(:,n)=d2adx(:,:,n)*aCur;
    end
    F=eye(xDim)+deltaT*dadx+(deltaT^2/2)*(term1+dadx*dadx);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function out = conv3fft(z1,z2)

z1 = single(z1);
z2 = single(z2);

siz1 = size(z1);
siz2 = size(z2);
siz = siz1+siz2-1;

z1x=size(z1,1);
z1y=size(z1,2);
z2x=size(z2,1);
z2y=size(z2,2);

out=real(ifftn(fftn(z1,siz).*fftn(z2,siz)));

p = ((siz2-1)+mod((siz2-1),2))/2;

out=out(p(1)+1:p(1)+siz1(1),p(2)+1:p(2)+siz1(2),p(3)+1:p(3)+siz1(3));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function showmesh3(node,elem,expr,varargin)
%% SHOWMESH3 displays a tetrahedron mesh in 3-D.
%
%    showmesh3(node,elem) displays a 3-dimensional tetrahderon mesh given
%    by node and elem matrices; see meshdoc for the data structure:
%    node and elem.
%
%    showmesh3(node,elem,expr) displays parts of the mesh specificed by the
%    expression. For example, showmesh3(node,elem,'~(x>=0 & y>=0)') only
%    shows the tetrahedron not in the first quadrant.
%
%    showmesh3(node,elem,expr,'param','value','param','value'...) allows
%    additional patch param/value pairs to be used when displaying the
%    mesh. For example, the default transparency parameter is set to 0.5.
%    You can overwrite this value by using the param pair ('FaceAlpha',
%    value). The value has to be a number between 0 and 1. Other parameters
%    include: 'Facecolor', 'Edgecolor' etc.
%
%    For meshes with large data, the 3-D graphics is very slow. You may use
%    showboundary3 to display the boundary surface mesh only.
%
%   Example:
%     % A mesh for a cube
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     [node,elem] = uniformbisect3(node,elem);
%     subplot(1,2,1);
%     showmesh3(node,elem); pause(1)
%     subplot(1,2,2);
%     showmesh3(node,elem,'~(x>=0 & y>=0 & z>=0)','FaceAlpha',0.25);
%     axis on; view([59,20])
%
%   See also showboundary3, showsolution3, showmesh.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (nargin >= 3) && (any(expr))
    x = node(:,1);  y = node(:,2);  z = node(:,3); %#ok<NASGU>
    incl = find(eval(expr));
    elem = elem(any(ismember(elem,incl),2),:);
end
h = tetramesh(elem(:,1:4),node,ones(size(elem,1),1));
set(h,'facecolor',[0.5 0.9 0.45],'edgecolor','k');
if nargin > 3
    set(h,varargin{1:end})
else % default display properties
    set(h,'FaceAlpha',0.4);
end
view(3);
axis off; axis equal; axis tight

repo: ifem


text:
function [speedup dsize ksize] = testcudaconv()
%TESTCUDACONV Benchmarks GPU vs CPU convolution.
%   [SPEEDUP DSIZE KSIZE] = TESTCUDACONV() runs 2-D
%   convolutions on random data using both CONV2 and CUDACONV.
%
%   SPEEDUP contains the relative speed of CONV2 and CUDACONV.
%   Thus if SPEEDUP(i,j) = 10, CUDACONV completed the convolution
%   of random DSIZE(i) x DSIZE(i) and KSIZE(j) x KSIZE(j)
%   matrices in 1/10 the time of CONV2.
%
%   DSIZE and KSIZE contain the number of elements in the data
%   and kernel matrices, respectively, for each convolution.

data_sizes		= 100:100:800;
kernel_sizes 	= 10:20:90;
nrepeats		= 3;

gpu_results		= zeros(numel(data_sizes),numel(kernel_sizes));
cpu_results		= zeros(numel(data_sizes),numel(kernel_sizes));

for ds = 1:numel(data_sizes)
	disp(sprintf('Data size: %d elements\n',data_sizes(ds)^2));
	data = rand(data_sizes(ds));
	for ks = 1:numel(kernel_sizes)
		disp(sprintf(' Kernel size: %d elements\n',kernel_sizes(ks)^2));
		kernel = rand(kernel_sizes(ks));

		gpu_temp = zeros(1,nrepeats);
		cpu_temp = zeros(1,nrepeats);

		for t = [1 1:nrepeats]
			tic;
			temp = cudaconv(data,kernel);
			gpu_temp(t) = toc;

			tic;
			temp = conv2(data,kernel,'same');
			cpu_temp(t) = toc;
		end

		gpu_results(ds,ks) = mean(gpu_temp);
		cpu_results(ds,ks) = mean(cpu_temp);
	end
end

speedup = cpu_results./gpu_results;

dsize = data_sizes.^2;
ksize = kernel_sizes.^2;

[x y] = meshgrid(ksize,dsize);

figure;
[c h] = contourf(x,y,speedup,[1 5:5:100]); colorbar;
clabel(c,h,'color',[0 0 0],'FontSize',13,'LabelSpacing',400);
xlabel('Kernel size (# of elements)')
ylabel('Data size (# of elements)')
title('Relative speedup (GPU vs. CPU)')

%dlmwrite('performance2.dat',speedup,',');

author: Sable

repo: mcbench-benchmarks


text:
function el = export_structure_set_module_field(args)
%"export_structure_set_module_field"
%   Given a single scan, return a properly populated structure_set module
%   field for use with RT Structure set IODs.  See structure_set_module_tags.m.
%
%   For speed, tag must be a decimal representation of the 8 digit
%   hexidecimal DICOM tag desired, ie instead of '00100010', pass
%   hex2dec('00100010');
%
%   Arguments are passed in a structure, arg:
%       arg.tag         = decimal tag of field to fill
%       arg.data        = CERR structure(s) to fill from
%       arg.template    = an empty template of the module created by the
%                         function build_module_template.m
%
%   This function requires arg.data = {structureS};
%
%JRA 06/19/06
%
%Usage:
%   dcmobj = export_patient_module_field(args)
%
% 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/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag         = args.tag;
structuresS = args.data{1};
scansS      = args.data{2};
template    = args.template;

switch tag
    %Class 1 Tags -- Required, must have data.
    case 805699586  %3006,0002 Structure Set Label.
        data = 'CERR Export';
        el = data2dcmElement(data, tag);

    %Class 2 Tags -- Must be present, can be NULL.
    case 805699592  %3006,0008 Structure Set Date
        data = datestr(now, 29);
        el = data2dcmElement(data, tag);

    case 805699593  %3006,0009 Structure Set Time
        data = datestr(now, 13);
        el = data2dcmElement(data, tag);

    %Class 3 Tags -- presence is optional, currently undefined.
    case 805699588  %3006,0004 Structure Set Name
        %Currently unsupported.

    case 805699590  %3006,0006 Structure Set Description
        %Currently unsupported.

    case   2097171  %3006,0013 Instance Number
        %Currently unsupported.

    case 805699600  %3006,0010 Referenced Frame of Reference Sequence
        templateEl = template.getValue(tag);

        fHandle = @export_referenced_frame_of_reference_sequence;

        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        %Get unique frame of reference UIDs used by structures.
        [FORs,indV]    = unique({structuresS.Frame_Of_Reference_UID});
        nFORs   = length(FORs);

        for i=1:nFORs
            dcmobj = export_sequence(fHandle, templateEl, {FORs(i), scansS, structuresS(indV(i))});
            %dcmobj = export_sequence(fHandle, templateEl, {FORs(i), scansS});
            el.add(i-1, dcmobj);
        end

        el = el.getParent();

    case 805699616  %3006,0020 Structure Set ROI Sequence
        templateEl = template.getValue(tag);
        fHandle = @export_structure_set_ROI_sequence;

        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        nStructures = length(structuresS);

        for i=1:nStructures
            dcmobj = export_sequence(fHandle, templateEl, {structuresS(i), i});
            %dcmobj = export_sequence(fHandle, tag, {structuresS(i), i});
            el.add(i-1, dcmobj);
        end

        el = el.getParent();
    %Class 1C Tags -- presence is required under special circumstances

    %Class 2C Tags.

    otherwise
        warning(['No methods exist to populate DICOM structure_set module field ' dec2hex(tag,8) '.']);
end


text:
function lstOut = pyrLst1OP(lstIn, fun)
%
%
%        lstOut = pyrLst1OP(lstIn, fun)
%
%
%        Input:
%           -lstIn: a list of image pyramids
%           -fun: a function to apply to lstIn
%
%        Output:
%           -lstOut: the result of the function
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

n = length(lstIn);

lstOut=[];

for i=1:n
    p = fun(lstIn(i));
    lstOut = [lstOut, p];
end

end

author: banterle


text:
function str = cvx_replicate_structure( str, sz )
%CVX_REPLICATE_STRUCTURE Apply matrix structure to N-D arrays.

szs  = size( str );
nmat = prod( sz );
omat = ones( 1, nmat );
[ r, c, v ] = find( str );
nelm = length( r );
oelm = ones( 1, nelm );
r = r( : );
nvec = ( 0 : nmat - 1 ) * szs( 1 );
r = r( :, omat ) + nvec( oelm, : );
c = c( : );
c = c( :, omat );
nvec = ( 0 : nmat - 1 ) * szs( 2 );
c = c( :, omat ) + nvec( oelm, : );
v = v( : );
v = v( :, omat );
str = sparse( r, c, v, nmat * szs( 1 ), nmat * szs( 2 ) );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [success, C] = cornerDetection(I,saveFileName)
% function C = cornerDetection(moviefile,frameRange,showResult)
% Detects corners
% frameRange is the range of the frames to be used to detect corners.
% Show results shows the result of corner detection.

% [readfcn,nframes,fid,~] = get_readframe_fcn(...
%   moviefile);
% %% Read a frame.
%
% if isempty(frameRange),
%   frameRange = 1:nframes;
% end
%
% rfr = randsample(frameRange,1);
% I = readfcn(rfr);

%%
CM = cornermetric(I(:,:,1));
bw = CM>0.000005;
bw([1:10 end-10:end],:) = 0;
bw(:,[1:10 end-10:end]) = 0;
bwl = bwlabel(bw);
bwcenter = regionprops(bwl,'Centroid','Area');
small2rem = find([bwcenter.Area]<3);

xc = []; yc = [];
for ndx = 1:numel(bwcenter)
xc(ndx) = bwcenter(ndx).Centroid(1);
yc(ndx) = bwcenter(ndx).Centroid(2);
end

dd = zeros(numel(xc),numel(xc));
for ndx = 1:numel(xc)
  dd(ndx,:) = sqrt( (xc - xc(ndx)).^2 + (yc-yc(ndx)).^2);
end

Isz = size(I,1);

success = false;
selidx = [];
for tol = 1:0.5:7
  for sz = round(0.7*Isz):(tol/2):round(Isz)
      ff = abs(dd-sz)<tol;
      if nnz(sum(ff)>1) >= 4
        curidx = find(sum(ff)>1);
        allcombs = nchoosek(curidx,4);
        for combn = 1:size(allcombs,1);
          selidx = allcombs(combn,:);
          curm = dd(selidx,selidx);
          if (nnz( abs(curm-sz)<tol) == 8) && (nnz( curm>1.4*sz)==4) ,
            success = true;
            break
          end
        end
        if success, break; end
      end
  end
  if success; break; end

end

if ~success,
  C = [];
  return;
end

topleft = selidx( (xc(selidx)<Isz/2) & (yc(selidx)<Isz/2));
topright = selidx( (xc(selidx)>Isz/2) & (yc(selidx)<Isz/2));
bottomleft = selidx( (xc(selidx)<Isz/2) & (yc(selidx)>Isz/2));
bottomright = selidx( (xc(selidx)>Isz/2) & (yc(selidx)>Isz/2));

if isempty(topleft) || isempty(topright) || isempty(bottomleft) || isempty(bottomright)
 success = false;
 C = [];
 return;
end

%{
xc(small2rem) = size(I,2)/2;
yc(small2rem) = size(I,1)/2;
minx = min(xc);
miny = min(yc);
maxx = max(xc);
maxy = max(yc);

tol = 10;
topleft = find( (abs(xc-minx)<tol) & (abs(yc-miny)<tol));
topright = find( (abs(xc-maxx)<tol) & (abs(yc-miny)<tol));
bottomleft = find( (abs(xc-minx)<tol) & (abs(yc-maxy)<tol));
bottomright = find( (abs(xc-maxx)<tol) & (abs(yc-maxy)<tol));

% C(1,:) = bwcenter(topleft).Centroid;
% C(2,:) = bwcenter(topright).Centroid;
% C(3,:) = bwcenter(bottomleft).Centroid;
% C(4,:) = bwcenter(bottomright).Centroid;

if isempty(topleft) || isempty(topright) || isempty(bottomleft) || isempty(bottomright),
continue;
end

if numel(topleft)>1 || numel(topright) > 1 || numel(bottomleft)>1 || numel(bottomright)>1,
  topleft = topleft(argmax([bwcenter(topleft).Area]));
  topright = topright(argmax([bwcenter(topright).Area]));
  bottomleft = bottomleft(argmax([bwcenter(bottomleft).Area]));
  bottomright = bottomright(argmax([bwcenter(bottomright).Area]));
end
success = true; break;
end
%}

if ~success,
  C = [];
  return;
end

%{
% selImg = (bwl == topleft) | (bwl == bottomleft) | (bwl == topright) | (bwl == bottomright);
%
%
% CM(~selImg) = 0;
% corner_peaks = imregionalmax(CM);
% corner_idx = find(corner_peaks == true);
% [xx yy] = ind2sub(size(I(:,:,1)),corner_idx);
% C = [];
%}
C(:,1) = yc([topleft topright bottomleft bottomright]);
C(:,2) = xc([topleft topright bottomleft bottomright]);

for ndx = 1:4
  I(round(C(ndx,1) + [-1:1]),round(C(ndx,2)+[-1:1]),1) = 256;
  I(round(C(ndx,1) + [-1:1]),round(C(ndx,2)+[-1:1]),2) = 0;
  I(round(C(ndx,1) + [-1:1]),round(C(ndx,2)+[-1:1]),3) = 0;
end
imwrite(I,saveFileName);

%% Detect corners based on square constraint.

%{

%}

%% Detect corners

%{
C = corner(I(:,:,1),'Harris',4,'SensitivityFactor',0.001);

%% Show Corners
imshow(I);
hold on
plot(C(:,1), C(:,2), 'r*');

%% Corner metric

hx = [];
figure;
hx(1) = subplot(1,3,1);
imshow(I);
title('Original Image');
CM = cornermetric(I(:,:,1));

CM_adjusted = imadjust(CM);
hx(2) = subplot(1,3,2);
imshow(CM_adjusted);
title('Corner Metric');

corner_peaks = imregionalmax(CM);
corner_idx = find(corner_peaks == true);
[r g b] = deal(I(:,:,1));
r(corner_idx) = 255;
g(corner_idx) = 255;
b(corner_idx) = 0;
RGB = cat(3,r,g,b);
hx(3) = subplot(1,3,3);
imshow(RGB);
title('Corner Points');
linkaxes(hx);
%}

author: kristinbranson

repo: JAABA


text:
diffusivityUnitStr = '(\mum^2/msec)';

if(~exist('dataDir','var'))
    dataDir = 'dti';
end

if(ispc)
    baseDir = '//171.64.204.10/biac2-wandell2/data/reading_longitude';
else
    baseDir = '/biac2/wandell2/data/reading_longitude';
end

subDir = fullfile(baseDir,dataDir,'*0*');
outFileName = ['mniRoiAnalysis_' dataDir];
outDir = fullfile(baseDir,'roiAnalyses');
if(~exist(outDir,'dir')) mkdir(outDir); end
dataSumFile = fullfile(outDir,[outFileName '_sum.mat']);
logFileName = fullfile(outDir,[outFileName '_log_' datestr(now,'yymmdd_HHMM') '.txt']);

if(~exist(dataSumFile,'file'))
    % Klingberg:
    % "The VOI had a volume of 960 mm^3, and was located within x = -36 to -26,
    % y = -50 to -10, and z = 0 to 32 mm relative the to anterior commissure."
    % (These are presumably MNI coordinates; peak at [-28 -20 28])
    %kCoordsMni = [-31 -20 26];
    % Beaulieu: -28 -14 24
    % Klingberg: -28 -20 28 (reported in Tal: -28 -22 29)
    % Deutsch: -28 -28 24 (reported in Tal: -28 -26 23)
    % Nagy: -27 -29 5 (reported in Tal: -27 -28 6)
    % Niogy: hard to tell. The slice shown has a coord above it, probably Tal, -28 -11 26, which is in MNI: -28 -12 28.
    clear c mniRois;
    c(1).name = 'Klingberg'; c(1).coords = [-28 -20 28];
    c(2).name = 'Deutsch'; c(2).coords = [-28 -28 24];
    c(3).name = 'Nagy'; c(3).coords = [-27 -29 5];
    c(4).name = 'Niogi'; c(4).coords = [-28 -12 28];
    for(ii=1:length(c))
        mniRois(ii,1) = dtiNewRoi([c(ii).name '03']); mniRois(ii,1).coords = dtiBuildSphereCoords(c(ii).coords, 3);
        mniRois(ii,2) = dtiNewRoi([c(ii).name '05']); mniRois(ii,2).coords = dtiBuildSphereCoords(c(ii).coords, 5);
    end
    mniRois = mniRois(:);

    excludeSubs = {'dh040607','an041018','hy040602','lg041019','mh040630','tk040817'};
    [f,sc] = findSubjects(subDir,'*_dt6_noMask',excludeSubs);
    N = length(f);

    goodSubs = ones(1,N);
    clear data;
    s=license('inuse'); userName = s(1).user;
    analysisDate = datestr(now,31);
    analyzedBy = userName;
    logFile = fopen(logFileName,'wt');
    fprintf(logFile, '* * * Analysis run by %s on %s * * *\n', userName, analysisDate);
    nSubs = 0;
    for(ii=1:N)
        tic;
        fname = f{ii};
        disp(['Processing ' fname '...']);
        fprintf(logFile, ['Processing ' fname ': ']);
        roiPath = fullfile(fileparts(fname), 'ROIs');
        fiberPath = fullfile(fileparts(fname), 'fibers');

        dt = load(f{ii});
        % We need to invert the spatial normalization
        def = dt.t1NormParams(1);
        if(isfield(def,'name')&&~strcmp(def.name,'MNI'))
            error(['Problem with MNI sn for ' fname '.']);
        end
        [def.deformX, def.deformY, def.deformZ] = mrAnatInvertSn(def.sn);
        def.inMat = inv(def.sn.VF.mat); % xform from acpc space to deformation field space

        data(ii).sc = sc{ii};
        % WARP MNI ROIS TO THIS BRAIN & GET THE TENSORS
        for(jj=1:length(mniRois))
            roi = mniRois(jj);
            roi.coords = mrAnatXformCoords(dt.t1NormParams.sn, roi.coords);
            roi.coords = unique(round(roi.coords),'rows');
            roi = dtiRoiClean(roi, 3, {'fillHoles', 'removeSat'});
            %dtiWriteRoi(roi, fullfile(roiPath, roi.name), 1, 'acpc');
            data(ii).rois(jj) = roi;
            % get the tensors for this ROI
            data(ii).dt6{jj} = dtiGetValFromTensors(dt.dt6, roi.coords, inv(dt.xformToAcPc), 'dt6', 'nearest');
        end
        fprintf('(%0.1f seconds) ',toc);
    end
    fclose(logFile);
    save(dataSumFile, 'data', 'f', 'mniRois', 'analysisDate', 'analyzedBy');
else
    disp(['Loading ' dataSumFile '...']);
    load(fullfile(outDir,[outFileName '_sum.mat']));
end

behavDataFile = fullfile(baseDir,'read_behav_measures.csv');
[bd,bdColNames] = dtiGetBehavioralData({data(:).sc},behavDataFile);
ageInd = strmatch('Age',bdColNames);
paInd = strmatch('Phonological Awareness',bdColNames);
widInd = strmatch('Word ID ss',bdColNames);
readerTypeInd = strmatch('Type of Reader',bdColNames);
goodReaders = bd(:,readerTypeInd)'>=0;

%goodSubs=[1:43,45:length(f)];
goodSubs=[1:length(f)];

clear fa md rd ad pdd;
for(ii=1:length(data))
  for(jj=1:length(data(ii).dt6))
    [vec,val] = dtiEig(data(ii).dt6{jj});
    val = val./1000;
    fa(ii,jj) = mean(dtiComputeFA(val));
    md(ii,jj) = mean(val(:));
    rd(ii,jj) = mean((val(:,2)+val(:,3))./2);
    ad(ii,jj) = mean(val(:,1));
    pdd(ii,jj,:) = mean(squeeze(vec(:,:,1)));
  end
end

logFile = fopen(logFileName,'wt');
fprintf(logFile, '\nSTATS:\n');
roiNames = {mniRois.name};
bdInds = [1:3,5,7:20];
roiInds = [1:8];
valNames = {'fa','md','rd'};
for(valNum=1:length(valNames))
    eval(['val=' valNames{valNum} ';']);
    fprintf('\n\n%s:\n',upper(valNames{valNum}));
    fprintf(logFile,'\n\n%s:\n',upper(valNames{valNum}));
    for(ii=bdInds)
        [r,p] = corrcoef([bd(goodSubs,ii),val(goodSubs,:)]);
        r = r(2:end,1); p = p(2:end,1);
        fprintf('%031s:',bdColNames{ii});
        fprintf(logFile,'%031s:',bdColNames{ii});
        for(jj=roiInds)
            if(p(jj)<0.01) sig='***'; elseif(p(jj)<0.05) sig='** '; elseif(p(jj)<0.1) sig='*  '; else sig='   '; end
            fprintf('\t%s r=%+0.2f %s',roiNames{jj},r(jj),sig);
            fprintf(logFile,'\t%s r=%+0.2f %s',roiNames{jj},r(jj),sig);
        end
        fprintf('\n'); fprintf(logFile,'\n');
    end
end
fclose(logFile);

bdInds = [3,5,7,8,9,10];
for(jj=bdInds)
    figure; set(gcf,'name',bdColNames{jj});
    for(ii=1:4)
        subplot(2,2,ii);
        x = bd(goodSubs,jj);
        y = fa(goodSubs,ii);
        plot(x,y,'k.');
        xlabel(bdColNames{jj});
        ylabel(['FA (' roiNames{ii} ')']);
        [r,p] = corr(x,y);
        title(sprintf('r=%0.2f (p=%0.4f)',r,p));
    end
end

for(jj=bdInds)
    for(ii=1:4)
        x = bd(goodSubs,jj);
        y = fa(goodSubs,ii);
        xlabel(bdColNames{jj});
        ylabel(['FA (' roiNames{ii} ')']);
        [r,p] = corr(x,y);
        title(sprintf('r=%0.2f (p=%0.4f)',r,p));
    end
end

meanPdd = dtiDirMean(permute(pdd,[2 3 1]));
error('stop here');

%% Whole-brain analysis

diffusivityUnitStr = '(\mum^2/msec)';
dataDir = 'dti';
baseDir = '/biac2/wandell2/data/reading_longitude';
subDir = fullfile(baseDir,dataDir,'*0*');
outFileName = ['mniRoiAnalysis_' dataDir];
outDir = fullfile(baseDir,'roiAnalyses');
if(~exist(outDir,'dir')) mkdir(outDir); end
dataSumFile = fullfile(outDir,[outFileName '_sum.mat']);
logFileName = fullfile(outDir,[outFileName '_log_' datestr(now,'yymmdd_HHMM') '.txt']);
excludeSubs = {'dh040607','an041018','hy040602','lg041019','mh040630','tk040817'};
[f,sc] = findSubjects(subDir,'*_dt6_noMask',excludeSubs);
N = length(f);

addpath /home/bob/matlab/stats/

behavDataFile = fullfile(baseDir,'read_behav_measures.csv');
[bd,bdColNames] = dtiGetBehavioralData(sc,behavDataFile);
ageInd = strmatch('Age',bdColNames);
paInd = strmatch('Phonological Awareness',bdColNames);
widInd = strmatch('Word ID ss',bdColNames);
readerTypeInd = strmatch('Type of Reader',bdColNames);
goodReaders = bd(:,readerTypeInd)'>=0;
goodSubs=[1:length(f)];

ni = niftiRead(fullfile('/silver/scr1/mniWarpAnalysis/data/',[sc{1} '_FA_sn.nii.gz']));
xform = ni.qto_xyz;
fa = zeros([size(ni.data),N]);
fa(:,:,:,1) = ni.data;
for(ii=2:N)
   ni = niftiRead(fullfile('/silver/scr1/mniWarpAnalysis/data/',[sc{ii} '_FA_sn.nii.gz']));
   fa(:,:,:,ii) = ni.data;
end

mnFa = mean(fa,4);
%mask = mnFa>=0.25;
mask = dtiCleanImageMask(all(fa>=0.15,4));
maskInds = find(mask);
showMontage(mask);
faInd = zeros(N,length(maskInds));
for(ii=1:N)
   tmp = fa(:,:,:,ii);
   faInd(ii,:) = tmp(maskInds);
end

bdInd = paInd;
valName = 'FA';
%[p,s] = myStatTest(bd(:,bdInd),faInd,'r');
[p,s] = myStatTest(bd(:,bdInd),faInd,'k');

imP=zeros(size(mask)); imP(maskInds)=p;

slAc=[10:2:40];
slIm = inv(xform)*[ones(2,length(slAc));slAc;ones(1,length(slAc))];
slIm = round(slIm(3,:));

% Remove 26-connected regions samller than clustThresh voxels
clustThresh = 10;
pThresh = 0.05;

imR=zeros(size(mask)); imR(maskInds)=s;
threshIm = imP<pThresh & imR<0; threshIm(~mask) = 0;
threshIm = bwareaopen(threshIm,clustThresh,26);
imR(~threshIm) = 0;
mrAnatOverlayMontage(-imR, xform, mnFa, xform, autumn(256), [0.1 1], slAc, fullfile(outDir,'FA_PA_negCorr'), 3, 1, true, 1);
set(gcf,'name','negative correlations');

imR=zeros(size(mask)); imR(maskInds)=s;
threshIm = imP<pThresh & imR>0; threshIm(~mask) = 0;
threshIm = bwareaopen(threshIm,clustThresh,26);
imR(~threshIm) = 0;
%showMontage(imR);
mrAnatOverlayMontage(imR, xform, mnFa, xform, autumn(256), [0.1 1], slAc, fullfile(outDir,'FA_PA'), 3, 1, true, 1);
set(gcf,'name','positive correlations');

% Show stats for each cluster
[clustLabel,clustN] = bwlabeln(threshIm, 26);
fprintf('%s vs. %s: %d clusters passed the size threshold of %d.\n',bdColNames{bdInd},valName,clustN,clustThresh);
plotBest = false;
for(jj=1:clustN)
    % Get the original data for each cluster and summarize it.
    curClust = clustLabel==jj;
    nvox = sum(curClust(:));
    clustInds = find(curClust(:));
    [inMask,maskLoc] = ismember(clustInds,maskInds);
    [x,y,z] = ind2sub(size(mask),clustInds);
    bestVox = find(p(maskLoc)==min(p(maskLoc)));
    mniCoords = mrAnatXformCoords(xform, [x,y,z]);
    mnMni = mean(mniCoords);
    bestMni = mniCoords(bestVox,:);
    if(mean(s(maskLoc))<0) sSum = min(s(maskLoc));
    else sSum = max(s(maskLoc)); end
    fprintf('   Cluster %d: size = %d voxels; tau = %0.2f; p = %0.5f; mean MNI: [%d %d %d]\n',...
        jj, nvox, sSum, min(p(maskLoc)), round(mnMni));
    % Plot the mean or best data within the cluster vs. behavior
    if(plotBest)
        mnDataVal = faInd(:,maskLoc(bestVox(1)));
    else
        mnDataVal = mean(faInd(:,maskLoc),2);
    end
    bdSum = bd(:,bdInd);
    notNan = ~isnan(bdSum);
    [curP,curR,curDf] = myStatTest(mnDataVal(notNan),bdSum(notNan),'r');
    [curKP,curKR] = myStatTest(mnDataVal(notNan),bdSum(notNan),'k');
    figure;
    set(gcf,'Name',sprintf('%d: %s vs. %s Scatter',jj,valName,bdColNames{bdInd}));
    plot(mnDataVal(notNan),bdSum(notNan),'ko');
    axis square;
    %fprintf('[%d %d %d]: r=%0.2f (p=%0.1e), tau=%0.2f (p=%0.1e), n=%d\n',round(mnMni),curR,curP,curKR,curKP,nvox);
    xlabel(valName);
    ylabel(bdColNames{bdInd});
end

[sigN,sigInd] = fdr(p,0.2,'original','mean');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Example 6.4: Regressor selection problem
% Section 6.3.1, Figure 6.7
% Original by Lieven Vandenberghe
% Adapted for CVX Argyris Zymnis - 10/2005
%
% Solves
%        minimize   ||A*x-b||_2
%        subject to card(x) <= k
%
% where card(x) denotes the number of nonzero elements in x,
% by first solving (for some value of alpha close to ||x_ln||_1)
%        minimize   ||A*x-b||_2
%        subject to ||x||_1 <= alpha
%
% and iteratively decreasing alpha so as to get card(x) = k
% The sparsity pattern is then fixed in A and b and
%        minimize   ||A*x-b||_2
%
% is solved

rand('state',0);
randn('state',0);

m = 10;
n = 20;

A = randn(m,n);
b = A*[randn(round(m/2),1); zeros(n-round(m/2),1)];
b = b + 0.1*norm(b)*randn(m,1);

if (1) %%%%%%%%%%%%

% tradeoff curve for heuristic
%
% min.  ||Ax-b||_2
% s.t.  ||x||_1 <= alpha

xln = A'*((A*A')\b);
lnorm = norm(xln,1);
nopts = 100;
alphas = linspace(0,lnorm,nopts);
residuals_heur = norm(b);
card_heur = 0;

for k=2:(nopts-1)
  alpha = alphas(k);

  cvx_begin quiet
    variable x(n)
    minimize(norm(A*x-b))
    subject to
        norm(x,1) <= alpha; %#ok
  cvx_end

  x(abs(x) < 1e-3*max(abs(x))) = 0; %#ok
  ind = find(abs(x));
  sparsity = length(ind);
  fprintf(1,'Current sparsity pattern k = %d \n',sparsity);
  x = zeros(n,1);  x(ind) = A(:,ind)\b;
  card_heur = [card_heur, sparsity]; %#ok
  residuals_heur = [residuals_heur, norm(A*x-b)]; %#ok
end;

obj1 = norm(b);
obj2 = 0;
obj1 %#ok

i=1;
for k=1:m-1
  if any(card_heur == k)
     obj2(i+1) = k; %#ok
     obj1(i+1) = min(residuals_heur(card_heur ==k));
     i=i+1;
  end;
end;
obj2(i) = m;  obj1(i) = 0;

end; %%%%%%%%%%%%%%%%%%%

% globally optimal tradeoff

if (1) %%%%%%%%%%%%%

bestx = zeros(n,m);
bestres = zeros(1,m);

for k=1:m-1
  k %#ok
  % enumerate sparsity patterns with exactly k nonzeros
  bestres(k) = Inf;
  ind = 1:k %#ok
  nocases = 1;
  done = 0;
  while ~done
     done = 1;
     for i=0:k-1
       if (ind(k-i) < n-i),
          ind(k-i:k) = ind(k-i)+(1:i+1);
          done = 0;
          break;
       end;
     end;
     if done, break; end;
     x = zeros(n,1);
     x(ind) = A(:,ind)\b;
     if (norm(A*x-b) < bestres(k)),
        bestres(k) = norm(A*x-b);
        bestx(:,k) = x;
     end;
     nocases = nocases + 1;
  end;
  nocases %#ok
  factorial(n)/(factorial(n-k)*factorial(k)) %#ok
end;

x = A\b;
bestres(m) = norm(A*x-b);
bestres = [norm(b) bestres];

end; %%%%%%%%%

figure
hold off
obj1dbl =[];
obj2dbl =[];
for i=1:length(obj1)-1
  obj1dbl = [obj1dbl, obj1(i), obj1(i)]; %#ok
  obj2dbl = [obj2dbl, obj2(i), obj2(i+1)]; %#ok
end;
obj1dbl = [obj1dbl, obj1(length(obj1))];
obj2dbl = [obj2dbl, obj2(length(obj1))];

bestobj1 = bestres;
bestobj2 = 0:m;
bestobj1dbl =[];
bestobj2dbl =[];
for i=1:length(bestobj1)-1
  bestobj1dbl = [bestobj1dbl, bestobj1(i), bestobj1(i)]; %#ok
  bestobj2dbl = [bestobj2dbl, bestobj2(i), bestobj2(i+1)]; %#ok
end;
bestobj1dbl = [bestobj1dbl, bestobj1(length(bestobj1))];
bestobj2dbl = [bestobj2dbl, bestobj2(length(bestobj1))];

plot(obj1dbl,obj2dbl,'-', bestobj1dbl, bestobj2dbl,'--');
hold on
plot(obj1,obj2,'o', bestobj1, bestobj2,'o');
axis([0 ceil(2*norm(b))/2 0 m+1])
xlabel('x');
ylabel('y');
hold off

%print -deps sparse_regressor_global_helv.eps
%save regressor_results

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [dat] = get(dt,ind,fInd)

%   data=get(data)                      returns the data of the data object
%   data=get(data,examIndex,featIndex)  returns the data of the data object
%                                       for given indexes

empty=0;
if nargin>1
    empty=isempty(ind);
end;

if (nargin==1 | empty)
    ind=[1:size(dt.Y,1)];  %<--- return all examples
end;
if nargin<3
    fInd=[1:size(dt.X,2)]; %<--- return all features
end;

y=[];
if ~isempty(dt.Y)
    y=dt.Y(ind,:);
end;

x=[];
if ~isempty(dt.X)
        x=dt.X(ind,fInd);
end;

dat  = data(dt.algorithm.name,x,y);
dat.index = dt.index(ind);
dat.findex = dt.findex(fInd);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function pass = test_merge(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

d = domain.merge([-inf 1 10 inf], [-inf 1 5 inf]);
pass = all(d == [-Inf     1     5    10   Inf]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [B0 phi0 G I_o omega f0 upsilon0]=stvol3prior(ar,arvar,lambda1,lambda3,lambda4,n,m,p,T,k,q,gamma,priorexo)

% obtain B0
% start first by obtaining beta0
beta0=zeros(q,1);
for ii=1:n
beta0((ii-1)*k+ii,1)=ar(ii,1);
end

% if a prior for the exogenous variables is selected put it in here:
for ii=1:n
    beta0(k*ii)=priorexo(ii,1);
end

% reshape the vector to obtain the matrix B0
B0=reshape(beta0,k,n);

% next compute phi0
% set first phi0 as a k*k matrix of zeros
phi0=zeros(k,k);
% set the variance for coefficients on lagged values
for ii=1:n
   for jj=1:p
   phi0((jj-1)*n+ii,(jj-1)*n+ii)=(1/arvar(ii,1))*(lambda1/jj^lambda3)^2;
   end
end
% set the variance for exogenous variables
for ii=1:m
phi0(k-m+ii,k-m+ii)=(lambda1*lambda4(ii))^2;
end

% compute the G matrix
G=speye(T)-sparse(diag(gamma*ones(T-1,1),-1));
% set the value for omega
omega=10000;
% compute I_omega
I_o=sparse(diag([1/omega;ones(T-1,1)]));

% compute the series of f0 vector and upsilon0 matrices
f0=cell(n,1);
upsilon0=cell(n,1);
for ii=2:n
f0{ii,1}=zeros(ii-1,1);
upsilon0{ii,1}=10000*eye(ii-1);
end

author: european-central-bank

repo: BEAR-toolbox


text:
function s = ymd_to_s_julian ( y, m, d )

%*****************************************************************************80
%
%% YMD_TO_S_JULIAN writes a Julian YMD date into a string.
%
%  Format:
%
%    AD YYYY/MM/DD
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, the YMD date.
%
%    Output, string S, a representation of the date.
%
  if ( y < 0 )
    s = sprintf ( 'BC %d/%02d/%02d', -y, m, d );
  else
    s = sprintf ( 'AD %d/%02d/%02d', y, m, d );
  end

  return
end


text:
function [x,state] = struct_band(z,task,size_mat,band)
%STRUCT_BAND Band matrix.
%   [x,state] = struct_band(z,[],size_mat,band) generates a band matrix x
%   of size size_mat, where the diagonals band(1) to band(2) are filled
%   column by column with entries from the vector z. For example, if x is a
%   square matrix of order n, the vector z should have length
%   sum(n-abs(band(1):band(2))). The structure state stores information
%   which is reused in computing the right and left Jacobian-vector
%   products.
%
%   struct_band(z,task,size_mat,band) computes the right or left
%   Jacobian-vector product of this transformation, depending on the
%   structure task. Use the structure state and add the field 'r' of the
%   same shape as z or the field 'l' of the same shape as x to obtain the
%   structure task for computing the right and left Jacobian-vector
%   products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_diag, struct_tridiag, struct_tril, struct_triu.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end
if nargin < 3 || ~isvector(size_mat)
    error('struct_band:size_mat','Missing integer matrix order.');
end
if nargin < 4 || ~isvector(band)
    error('struct_band:band','Missing definition of band.');
end

if isempty(task) || (isempty(task.l) && isempty(task.r))
    state.idx = bsxfun(@plus,(size_mat(1)-1:-1:0).', ...
        (0:size_mat(2)-1)-size_mat(1)+1);
    state.idx = find(state.idx >= band(1) & state.idx <= band(2));
    x = zeros(size_mat);
    x(state.idx) = z;
elseif ~isempty(task.r)
    x = zeros(size_mat);
    x(task.idx) = task.r;
    state = [];
elseif ~isempty(task.l)
    x = task.l(task.idx);
    state = [];
end

end

author: andrewssobral

repo: lrslibrary


text:
function [ num1, num2 ] = myconvert( str )
% Function to convert the vector string to numbers
% By Kyriakos Tsourapas
% You may contact me through the Mathworks site
% University of Essex 2002

% CREATE THE NUMBERS FROM THE STRING
sign1    = str(1);
intpart1 = sprintf('%d%d', str(2:3));
decpart1 = sprintf('%d%d%d%d%d%d%d%d%d%d', str(4:13));
sign2    = str(14);
intpart2 = sprintf('%d%d', str(15:16));
decpart2 = sprintf('%d%d%d%d%d%d%d%d%d%d', str(17:26));

% CONVERT THEM TO DECIMAL FROM BINARY
num1 = bin2dec(intpart1) + bin2dec(decpart1)/1000;
num2 = bin2dec(intpart2) + bin2dec(decpart2)/1000;

if sign1 == 1
    num1 = - num1;
end

if sign2 == 1
    num2 = - num2;
end


text:
function output = ConvertPowerToDecibel(x)
output = 10 * log10(x);
end

author: google


text:
function v = grad(f)
%GRAD Gradient of a BALLFUN in cartesian coordinates.
%   GRAD(F) is the gradient of the BALLFUN F expressed in
%   cartesian coordinates.
%
%   This is shorthand for the command GRADIENT.
%
% See also DIV, CURL

% Copyright 2019 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

v = gradient(f);
end


text:
function test_suite=equation_test
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
    test_functions=localfunctions();
catch % no problem; early Matlab versions can use initTestSuite fine
end
initTestSuite;

function TestSetup
setenv('ISDISPLAY','0') % go faster! Fit only 2 voxels in FitData.m

function test_equation
disp('testing equations...')
MethodList = list_models;
for im = 1:length(MethodList)
    Model = str2func(MethodList{im}); Model = Model();
    if ~Model.voxelwise, continue; end
    disp(class(Model))
    ModelOpt = button2opts(Model.buttons,1);
    clear st Smodel
    for iopt=1:length(ModelOpt) % try all model options
        Model.options = ModelOpt(iopt);
        disp(['Testing ' class(Model) ' option:'])
        disp(Model.options)
        try Model = Model.UpdateFields; end
        try st{iopt} = Model.st; catch, try st{iopt} = mean([Model.lb(:),Model.ub(:)],2); catch, st{iopt} = ones(length(Model.xnames),1); end; end
        Smodel{iopt} = Model.equation(st{iopt});
        % CHECK CONSITENSY WITH PREVIOUS VERSIONS:
        if exist(['value_' class(Model) '.mat'],'file')
            % compare with Ground Truth
            GT = load(['value_' class(Model) '.mat']);
            [~,ModelOpttest,GTModelOpttest]=comp_struct(ModelOpt,GT.ModelOpt);
            if ~isempty(ModelOpttest) || ~isempty(GTModelOpttest)
                msg = [MethodList{im} ' buttons/options has changed' evalc('ModelOpttest, GTModelOpttest')];
            else
                msg = '';
            end
            assertVectorsAlmostEqual(st{iopt},GT.st{iopt},'relative',1e-4,[MethodList{im} ' starting point (st) changed... ' msg ' value_' MethodList{im} '.mat has to be regenerated.'])
            assertVectorsAlmostEqual(Smodel{iopt},GT.Smodel{iopt},'relative',1e-2,['Synthetic signal obtained from ' MethodList{im} ' equation is not consistent with previous versions... ' msg])
        elseif iopt == length(ModelOpt)
            save(['value_' class(Model) '.mat'],'Smodel','st','ModelOpt')
        end
    end
end

function TestTeardown
setenv('ISDISPLAY','') % go faster! Fit only 2 voxels in FitData.m

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [imageCha, ssim_map, metrics] = algo_ADMM_proxA_2D( obj,input )
% 2D variant of ADMM
% based on Afonso et al. paper on SALSA
%
% input:
% obj           CS reconstruction object (holding all parameters)
% input         struct containing recon parameters and image
%
% output:
% imageCha      reconstructed channel individual image
% ssim_map      structural similarity map
% metrics       evaluation metrics
%
% (c) Marc Fischer, Thomas Kuestner, May 2015
% -------------------------------------------------------------------------

%%
timer_proxA = tic;

%% variables:
% internal flags
flag_wavetree = true;
flag_fast = true;
flag_extendImage = true;

% internal variables:
L = 1;
t_old = 1;
NLTV_struct.kernelratio = 3;
NLTV_struct.windowratio = 6;
NLTV_struct.nThreads = 1;  % mind this option if used with -singleCompThread on BWCluster
itrNLTV = obj.iNINNER - 20;
% chambolle tv:
parsin.MAXITER=100; parsin.tv='iso'; % 'iso' or 'l1'

% from obj:
% maxitr = obj.maxitr;
maxitr = obj.iNINNER;
% n1 = obj.measPara.dim(1);
% n2 = obj.measPara.dim(2);
n1 = input.n1;
n2 = input.n2;
% nSlices = obj.measPara.dim(3);
nCha = obj.measPara.dim(5);
lambdaWave = obj.lambda;
lambdaTV = obj.lambdaTV;
lambdaGroup = obj.lambdaGroup;
NLTV_struct.filterstrength = obj.lambdaNLTV_h; % 0.03 % converted from NLTV h = 0.01 %old: used: 3e-10
lambdaNLTV = obj.lambdaNLTV;
lambdaNLTV_h = obj.lambdaNLTV_h;
mue = obj.mue;
regularizerWeights = obj.regularizerWeights;
flagTV = obj.flagTV;
flagTV_iso = obj.flagTV_iso;
flagWave = obj.flagWave;
flagGroup = obj.flagGroup;
flagNLTV = obj.flagNLTV;
flagSBNLTV = obj.flagSBNLTV;
flagRealAndImag = obj.flagRealAndImag;
waveletStages = obj.trafo.waveletStages;
waveletFilterName_l1 = obj.trafo.waveletFilterName_l1;
waveletFilterName_l12 = obj.trafo.waveletFilterName_l12;

% from input:
b=input.b;
mask = input.mask;
G_prox = input.G_prox;
Gt_prox = input.Gt_prox;
groupnorm_index = input.groupnorm_index;
waveS_l1 = input.waveS_l1;
waveS_l12 = input.waveS_l12;
waveS_l12proxA = input.waveS_l12proxA;
proxA_extend_y = waveS_l12proxA(waveletStages+2,1) - waveS_l12(waveletStages+2,1);
proxA_extend_x = waveS_l12proxA(waveletStages+2,2) - waveS_l12(waveletStages+2,2);
x_proxA = cell(1,nCha);

% im_ref = input.im_ref;
% im_ref_full = zeros(n1,n2);
% for j = 1:nCha
%     im_ref_full = im_ref_full + abs(im_ref{1,j}).^2;
% end;
% im_ref_full = sqrt(im_ref_full);
clear input

% initialize cells/vectors
for j=1:nCha
    FTb{1,j} = iFFT2D(b{1,j});
    x_cmplx{1,j} = real(FTb{1,j});
    x_cmplx{1,j+nCha} = imag(FTb{1,j});
end;
% starting point:
x = FTb; % (x0 = FTb)
rhs = FTb; % different start (x0 = F^-1 K^-1 F FTb); %used atm
% rhs: right hand side

x_g_helper = cell(2,nCha);
x_wave =  cell(1,2*nCha);
x_g_proxA = x_wave;
x_wave_helper = x_wave;
v_wave = x_wave;
d_wave = x_wave;
x_tv = x_wave;
v_tv = x_wave;
d_tv = x_wave;
x_group = x_wave;
v_group = x_wave;
d_group = x_wave;
x_nltv = x_wave;
v_nltv = x_wave;
d_nltv = x_wave;
K_inv = x_wave;
for j = 1:2*nCha
    x_wave{1,j} = 0;
    v_wave{1,j} = 0;
    d_wave{1,j} = 0;
    x_tv{1,j} = zeros(n1,n2);
    v_tv{1,j} = 0;
    d_tv{1,j} = 0;
    x_group{1,j} = 0;
    v_group{1,j} = 0;
    d_group{1,j} = 0;
    x_nltv{1,j} = 0;
    v_nltv{1,j} = 0;
    d_nltv{1,j} = 0;
end;

%% MAD dependent lambdas:
flag_MAD = true;
if flag_MAD
    x_wavedec = cell(1,nCha);
    threshold = zeros(1:nCha);
    if flagRealAndImag
        for j=1:nCha
            x_wavedec{1,j} = wavedec2(abs(x{1,j}),waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
            x_wave_fine_scale = size(x_wavedec{1,j},2) - (3*waveS_l1(waveletStages+1,1)*waveS_l1(waveletStages+1,2));
            threshold(j) = mad(x_wavedec{1,j}(x_wave_fine_scale:end),1);
        end;
    else
        for j=1:nCha
            x_wavedec{1,j} = wavedec2(real(x{1,j}),waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
            x_wavedec{1,j+nCha} = wavedec2(imag(x{1,j}),waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
            x_wave_fine_scale_real = size(x_wavedec{1,j},2) - (3*waveS_l1(waveletStages+1,1)*waveS_l1(waveletStages+1,2));
            threshold(j) = mad(x_wavedec{1,j}(x_wave_fine_scale_real:end),1);
            x_wave_fine_scale_imag = size(x_wavedec{1,j},2) - (3*waveS_l1(waveletStages+1,1)*waveS_l1(waveletStages+1,2));
            threshold(j+nCha) = mad(x_wavedec{1,j}(x_wave_fine_scale_imag:end),1);
        end;
    end;
    clear x_wavedec
else
    threshold(1:nCha) = 1;
end;

if flagRealAndImag
    threshold_group = 0;
    for j = 1:nCha
        threshold_wave(j) = lambdaWave * threshold(j)/2 * 2/L;
        threshold_TV(j) = lambdaTV * threshold(j) * 2/L;
        threshold_TV(j+nCha) = lambdaTV * threshold(j) * 2/L;
        threshold_group = threshold_group + lambdaGroup * threshold(j) * 2/L;
        threshold_NLTV(j) = lambdaNLTV; % * threshold(j) * 2/L;
        threshold_NLTV(j+nCha) = lambdaNLTV;
        threshold_NLTV_h(j) = lambdaNLTV_h; %*threshold(j); % adjust carefully or NLTV won't find a solution. lambdaNLTV_h should be < 0.01
        threshold_NLTV_h(j+nCha) = lambdaNLTV_h;
    end;
    threshold_group = threshold_group /3;
else
    threshold_group_real = 0;
    threshold_group_imag = 0;
    for j = 1:nCha
        threshold_wave(j) = lambdaWave * threshold(j) * 2/L;
        threshold_wave(j+nCha) = lambdaWave * threshold(j+nCha) * 2/L;
        threshold_TV(j) = lambdaTV * threshold(j) * 2/L;
        threshold_TV(j+nCha) = lambdaTV * threshold(j+nCha) * 2/L;
        threshold_group_real = threshold_group_real + lambdaGroup * threshold(j) * 2/L;
        threshold_group_imag = threshold_group_imag + lambdaGroup * threshold(j) * 2/L;
        threshold_NLTV(j) = lambdaNLTV; % * threshold(j) * 2/L;
        threshold_NLTV(j+nCha) = lambdaNLTV;
        threshold_NLTV_h(j) = lambdaNLTV_h; %*threshold(j); % adjust carefully or NLTV won't find a solution. lambdaNLTV_h should be < 0.01
        threshold_NLTV_h(j+nCha) = lambdaNLTV_h;
    end;
    threshold_group_real = threshold_group_real /3;
    threshold_group_imag = threshold_group_imag /3;
end;

%% initialize metrics:
    itr = 0;
   	metrics.xtime(itr+1)= 0;
%     [metrics, ssim_map{1,1}] = get_metrics_itr( im_ref, im_ref_full, x, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );
    ssim_map = [];

%% recon
% inverse for x - step:
for j = 1:nCha
    K_inv{1,j} = 1./(mue*mask + (flagWave*regularizerWeights(1) + flagTV*regularizerWeights(2) +  flagGroup*regularizerWeights(3))*ones(n1,n2)); % flagNLTV changes K_inv -> in separate
end;

dispProgress('ADMM', 0, maxitr);
for itr = 1:maxitr  % total iter counter

%% x - step
    if itr > 1 % for i = 1 already set (x{1,j} = FTb{1,j})
        for j = 1:2*nCha
            % rhs part:
            if flagWave
                x_wave{1,j} = waverec2((v_wave{1,j} - d_wave{1,j}),waveS_l1,waveletFilterName_l1); end;
            if flagTV
                x_tv{1,j} = v_tv{1,j} - d_tv{1,j}; end;
            if flagGroup
                x_g_proxA{1,j} = waverec2((v_group{1,j} - d_group{1,j}),waveS_l12proxA,waveletFilterName_l12);
                x_group{1,j} =  x_g_proxA{1,j}(1:end-proxA_extend_y,1:end-proxA_extend_x);
            end;
            if flagNLTV
                if itr == itrNLTV
                    K_inv{1,j} = 1./(mue*mask + (flagWave*regularizerWeights(1) + flagTV*regularizerWeights(2) +  flagGroup*regularizerWeights(3) + flagNLTV*regularizerWeights(4))*ones(n1,n2));
                    x_nltv{1,j} = x{1,j};
                elseif itr >= itrNLTV
                    x_nltv{1,j} = v_nltv{1,j} - d_nltv{1,j};
                end;
            end;
        end;
        for j = 1:nCha
            rhs{1,j} = mue*FTb{1,j} + ( flagWave*regularizerWeights(1).*(x_wave{1,j}+1i*x_wave{1,j+nCha}) + flagTV*regularizerWeights(2).*(x_tv{1,j}+1i*x_tv{1,j+nCha}) + flagGroup*regularizerWeights(3).*(x_group{1,j}+1i*x_group{1,j+nCha}) + flagNLTV*regularizerWeights(4).*(x_nltv{1,j}+1i*x_nltv{1,j+nCha}));
            x{1,j} = iFFT2D(K_inv{1,j}.*FFT2D(rhs{1,j}));
            x_cmplx{1,j} = real(x{1,j});
            x_cmplx{1,j+nCha} = imag(x{1,j});
        end;
    end;

%% auxiliary variables:
%% l1-Wavelet
    if flagWave
        for j = 1:2*nCha
            x_wave_helper{1,j} = wavedec2(x_cmplx{1,j},waveletStages,waveletFilterName_l1);
        end;
        if flagRealAndImag
            for j = 1:nCha
                [v_wave{1,j}, v_wave{1,j+nCha}] = groupthresh_2vec(x_wave_helper{1,j} + d_wave{1,j},x_wave_helper{1,j+nCha} + d_wave{1,j+nCha},threshold_wave(j));
            end;
        else
            for j = 1:2*nCha
                v_wave{1,j} = softthresh_real(x_wave_helper{1,j} + d_wave{1,j},threshold_wave(j));
            end;
        end;
        for j = 1:2*nCha
            d_wave{1,j} = d_wave{1,j} + x_wave_helper{1,j} - v_wave{1,j};
        end;
    end;

%% TV:
    if flagTV
        if ~flagTV_iso
            for j = 1:2*nCha
                    v_tv{1,j} = MTV_2D((x_cmplx{1,j} + d_tv{1,j}), threshold_TV(j), n1, n2);
            end;
        else
            for j = 1:2*nCha
                if (itr==1)
                    [v_tv{1,j}, P]=denoise_TV_One((x_cmplx{1,j} + d_tv{1,j}), threshold_TV(j),-inf,inf,[],parsin);
                else
                    [v_tv{1,j}, P]=denoise_TV_One((x_cmplx{1,j} + d_tv{1,j}), threshold_TV(j),-inf,inf,P,parsin);
                end;
            end;
        end;
        d_tv{1,j} = d_tv{1,j} + x_cmplx{1,j} - v_tv{1,j};
    end;

%% NLTV
    if flagNLTV
        if itr >= itrNLTV
            if flagSBNLTV
                for j = 1:2*nCha
                    v_nltv{1,j} = SB_NLTVfunc_slim_rescale( (x_cmplx{1,j} + d_nltv{1,j}),n1,n2, threshold_NLTV(j), threshold_NLTV_h(j) );
                end;
            else
                for j = 1:2*nCha
                    % if mod(itr,5) == 0 || itr == 1
                    v_nltv{1,j} = NLMF((x_cmplx{1,j} + d_nltv{1,j}),NLTV_struct);
                    v_nltv{1,j} = (L.*(x_cmplx{1,j} + d_nltv{1,j}) + 2*threshold_NLTV(j)*v_nltv{1,j})./(L+2*threshold_NLTV(j));
                    % end;
                end;
            end;
            d_nltv{1,j} = d_nltv{1,j} + x_cmplx{1,j} - v_nltv{1,j};
        else
            v_nltv{1,j} = 0;
            d_nltv{1,j} = 0;
        end;

    end;

%% l12-Wavelet
    if flagGroup
        for j = 1:2*nCha
            if flag_extendImage
                x_proxA{1,j} = extend_image(x_cmplx{1,j}, waveS_l12proxA, waveletStages, proxA_extend_y, proxA_extend_x);
                x_g_helper{1,j} = wavedec2(x_proxA{1,j},waveletStages,waveletFilterName_l12) + d_group{1,j};
            else
                x_g_helper{1,j} = wavedec2(x_cmplx{1,j},waveletStages,waveletFilterName_l12) + d_group{1,j};
            end;

            if flag_wavetree
                x_g_helper{2,j} = (G_prox*x_g_helper{1,j}')';
            else
                x_g_helper{2,j} = zeros(1,size(x_g_helper{1,j},2));
            end;
        end;
        if flagRealAndImag
            v_group(1,:) = softthresh_proxA_cha(x_g_helper,threshold_group,2*nCha,waveS_l12proxA,groupnorm_index);
        else
            v_group(1,1:nCha) = softthresh_proxA_cha(x_g_helper(:,1:nCha),threshold_group_real,nCha,waveS_l12proxA,groupnorm_index);
            v_group(1,nCha+1:2*nCha) = softthresh_proxA_cha(x_g_helper(:,1:nCha),threshold_group_imag,nCha,waveS_l12proxA,groupnorm_index);
        end;
        for j = 1:2*nCha
            d_group{1,j} = d_group{1,j} + x_g_helper{1,j} - v_group{1,j};
        end;
    end;

%% metrics of current itr:
%     disp(itr);
    dispProgress('ADMM', itr/maxitr);

   	metrics.xtime(itr+1)= toc(timer_proxA);
%     [metrics, ssim_map{1,2}] = get_metrics_itr( im_ref, im_ref_full, x, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );

end;
dispProgress('ADMM', 'Close');

    imageCha = x;
    for j = 1:nCha
        imageCha{1,j} = turn_image( imageCha{1,j} );
    end;
%     for j = 1:nCha+1
%         ssim_map{1,1}{1,j} = turn_image( ssim_map{1,1}{1,j} );
%         ssim_map{1,2}{1,j} = turn_image( ssim_map{1,2}{1,j} );
%     end;

end

author: thomaskuestner


text:
function varargout = process_select_tag( varargin )
% PROCESS_SELECT_TAG: Keep only the files within the current selection that have the specified tag in their filename.
%
% USAGE:  sProcess = process_select_tag('GetDescription')
%                    process_select_tag('Run', sProcess, sInputs)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Select files: By tag';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'File';
    sProcess.Index       = 1013;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/SelectFiles';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data', 'results', 'timefreq', 'matrix', 'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.OutputTypes = {'raw', 'data', 'results', 'timefreq', 'matrix', 'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 0;
    % Definition of the options
    % === TARGET
    sProcess.options.tag.Comment = 'Search: ';
    sProcess.options.tag.Type    = 'text';
    sProcess.options.tag.Value   = '';
    % === FILENAME / COMMENT
    sProcess.options.label1.Comment = 'Where to look for:';
    sProcess.options.label1.Type    = 'label';
    sProcess.options.search.Comment = {'Search the file paths', 'Search the file names', 'Search the names of the parent file'};
    sProcess.options.search.Type    = 'radio';
    sProcess.options.search.Value   = 2;
    % === SELECT / IGNORE
    sProcess.options.label2.Comment = 'What to do with these files:';
    sProcess.options.label2.Type    = 'label';
    sProcess.options.select.Comment = {'Select only the files with the tag', 'Ignore the files with the tag'};
    sProcess.options.select.Type    = 'radio';
    sProcess.options.select.Value   = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    % Option: Tag
    tag = strtrim(sProcess.options.tag.Value);
    if isempty(tag)
        tag = 'Not defined';
    end
    % Option: Filename/comment
    if isfield(sProcess.options, 'search') && isfield(sProcess.options.search, 'Value') && ~isempty(sProcess.options.search.Value)
        switch (sProcess.options.search.Value)
            case 1,  Method = 'filename';
            case 2,  Method = 'comment';
            case 3,  Method = 'parent';
        end
    else
        Method = 'comment';
    end
    % Option: Ignore/select
    if isfield(sProcess.options, 'select') && isfield(sProcess.options.select, 'Value') && ~isempty(sProcess.options.select.Value)
        isSelect = isequal(sProcess.options.select.Value, 1);
    else
        isSelect = 1;
    end
    % Assemble comment
    if isSelect
        Comment = 'Select';
    else
        Comment = 'Ignore';
    end
    switch (Method)
        case 'filename',  Comment = [Comment ' file paths with tag: ' tag];
        case 'comment',   Comment = [Comment ' file names with tag: ' tag];
        case 'parent',    Comment = [Comment ' parent names with tag: ' tag];
    end
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Get new tag
    tag = sProcess.options.tag.Value;
    if isempty(tag)
        bst_report('Error', sProcess, [], 'Tag is not defined.');
        OutputFiles = {};
        return
    end
    % Error: no input files
    if isempty(sInputs)
        bst_report('Error', sProcess, [], 'No input files.');
        OutputFiles = {};
        return;
    end
    % Option: Filename/comment
    if isfield(sProcess.options, 'search') && isfield(sProcess.options.search, 'Value') && ~isempty(sProcess.options.search.Value)
        switch (sProcess.options.search.Value)
            case 1,  Method = 'filename';
            case 2,  Method = 'comment';
            case 3,  Method = 'parent';
        end
    else
        Method = 'comment';
    end
    % Option: Ignore/select
    if isfield(sProcess.options, 'select') && isfield(sProcess.options.select, 'Value') && ~isempty(sProcess.options.select.Value)
        isSelect = isequal(sProcess.options.select.Value, 1);
    else
        isSelect = 0;
    end

    % Search filenames/comments
    switch (Method)
        case 'filename'
            isTag = ~cellfun(@(c)isempty(strfind(upper(c),upper(tag))), {sInputs.FileName});
        case 'comment'
            isTag = ~cellfun(@(c)isempty(strfind(upper(c),upper(tag))), {sInputs.Comment});
        case 'parent'
            isTag = zeros(1, length(sInputs));
            upTag = upper(tag);
            % Search all the parent files one by one
            for i = 1:length(sInputs)
                % Do not go further if there is no parent
                if isempty(sInputs(i).DataFile)
                    continue;
                end
                % What to check depends on the file type
                switch (file_gettype(sInputs(i).DataFile))
                    case 'data'
                        % Find the file in database
                        [sStudy, iStudy, iData] = bst_get('DataFile', sInputs(i).DataFile, sInputs(i).iStudy);
                        % Check the comment
                        isTag(i) = ~isempty(strfind(upper(sStudy.Data(iData).Comment), upTag));
                    case {'results', 'link'}
                        % Find the file in database
                        [sStudy, iStudy, iResult] = bst_get('ResultsFile', sInputs(i).DataFile, sInputs(i).iStudy);
                        % Check the comment
                        isTag(i) = ~isempty(strfind(upper(sStudy.Result(iResult).Comment), upTag));
                        % If the file is not found but there is another parent level
                        if ~isTag(i) && ~isempty(sStudy.Result(iResult).DataFile) && strcmpi(file_gettype(sStudy.Result(iResult).DataFile), 'data')
                            [sStudy, iStudy, iData] = bst_get('DataFile', sStudy.Result(iResult).DataFile, sInputs(i).iStudy);
                            isTag(i) = ~isempty(strfind(upper(sStudy.Data(iData).Comment), upTag));
                        end
                    case 'matrix'
                        % Find the file in database
                        [sStudy, iStudy, iMatrix] = bst_get('MatrixFile', sInputs(i).DataFile, sInputs(i).iStudy);
                        % Check the comment
                        isTag(i) = ~isempty(strfind(upper(sStudy.Matrix(iMatrix).Comment), upTag));
                end
            end
    end
    % Ignore or select
    if isSelect
        iFiles = find(isTag);
    else
        iFiles = find(~isTag);
    end
    % Warning: nothing found
    if isempty(iFiles)
        bst_report('Error', sProcess, [], ['No files found for tag: ' tag]);
        OutputFiles = {};
        return;
    end
    % Report information
    strReport = sprintf('Files selected for tag "%s": %d/%d', tag, length(iFiles), length(sInputs));
    bst_report('Info', sProcess, [], strReport);
    % Return only the filenames that have the specific tag
    OutputFiles = {sInputs(iFiles).FileName};
end

author: brainstorm-tools


text:
function [Hdraw,HvarsDraw,phi_Hdraw,Adraw,h0,Kh0]=sampleH_update(yData,Psi,A,B,phi_H,HvarsOld,priorValues,dataValues,h0,Kh0)
% this function samples the stochastic covariance matrix H. The steps are
% (1) the states s, (2) the log variances, (3) the phi parameters, and (4)
% the matrix A.
%yData=YData
%Psi = PsiDraw_prop; %draw for the local mean
%A=Adraw;            %draw for the A matrix (constant part of H=VCV_VAR),
%H=A^-1*Lambda*A^-1'
%B = Bdraw; %VAR coefficients
%phi_H =phi_Hdraw; %volatitilty of innovations in random walk for the elemnts of Lambda (stochastic volatility parameters)
%HvarsOld = HvarsDraw; %H from previous iteration
%% Initialize
offset_c=priorValues.offset_c; %constant for log transformation

[T,M]=size(yData);
n=M;
p = size(B,1)/M; %lags of the B vector that has dimensions MxM*p (no constant), each column is one regression

% obtain prior data
%startMeanVector=priorValues.mean_ln_h0;  %log mean of the initial state (variance scaling parameters) as residuals from an AR(4) in the training sample
%startVarVector=priorValues.var_ln_h0;    %variance of the initial state for the elements of lambda

% obtain prior data
b0=priorValues.mean_ln_h0;  %log mean of the initial state (variance scaling parameters) as residuals from an AR(4) in the training sample
a0=eye(M)/priorValues.var_ln_h0(1,1);    %variance of the initial state for the elements of lambda

priorPhi_H=priorValues.phi_h;            %centering parameter for the inverse gamma of the variance of the innovations governing the random walk for lambda
priorD_H  =priorValues.d_h;              %scaling parameter for the inverse gamma distribution for the variance of the innovations governing the random walk for lambda

S_h = priorPhi_H*ones(n,1);
nu_h = priorD_H*ones(n,1);

%% Prepare Data

Y_Psi=yData-Psi;  %subtract the local mean
Y_Psi(1:p,:)=yData(1:p,:)-ones(p,1)*mean(yData(1:p,:)); %also generate initial conditions for the construction of the lagmatrix
% X_Psi = lagmatrix(Y_Psi,1:p); %create RHS of the VAR part
X_Psi = bear.lagx(Y_Psi,p-1);
X_Psi = X_Psi(1:end-1,:);
% X_Psi = X_Psi(p+1:end,:);     %remove the first p rows of RHS
Y_Psi=Y_Psi(p+1:end,:);       %and do so for LHS

E=Y_Psi-X_Psi*B;              %VAR residuals
Escaled=E*A';                 %transform such that VAR residuals have variance Lambda

Hvars=zeros(T,M);             %sampled states (diagonal elements of lambda)
%phi_Hdraw=zeros(1,M);         %variance of the random walk process governing the elements of lambda

%% Sample (States,Vars,phi)

for i=1:M
    % prepare for each i
    residsTemp=Escaled(:,i);
    yStar=log(residsTemp.^2+offset_c); %transform scaled residuals
    phi=phi_H(i); %previous draw for the variance of this particular variable
    lnSigma2=log(HvarsOld(p+1:T,i)); % note log => h = ln sigma2

    startMean=h0(i); %variable specific mean of the initial condition for the  volatility
    startVar=Kh0(i,i);   %variable specific variance of the disturbances of the random walk for the volatility process

    % sample states
    [yStarAdj, Ht] = bear.statesMix(yStar,lnSigma2);

    % sample log variances
    [logVarsDraw_Hi]=bear.KF_CKsimSV(yStarAdj,Ht,phi,startMean,startVar);
    Hvars(:,i)=exp([zeros(1,p) logVarsDraw_Hi']');

%     % sample phi
%     [phiDraw_Hi]=samplePhi(logVarsDraw_Hi,priorD_H,priorPhi_H);
%     phi_Hdraw(i)=phiDraw_Hi;
end

  %sample phi
%   for kk=1:n
%    [phiDraw_Hi]=samplePhi([h0(kk,1); log(Hvars(p+1:end,kk))],priorD_H,priorPhi_H);
%    phi_Hdraw(1,kk)=phiDraw_Hi;
%   end

for kk=1:n
logSVseries=[h0(kk,1); log(Hvars(p+1:end,kk))];
  %% Initialize
T1=length(logSVseries);
vDiffSV=(logSVseries(2:T1)-logSVseries(1:T1-1));

%% posteror parameter values
a1=priorD_H*priorPhi_H+sum(vDiffSV.^2); %scale parameter
b1=priorD_H+(T1-1); %shape parameter

%% Take a draw form the posterior distribution
gammaDraw=bear.grandn(b1/2,2/a1);
phiDrawkk=1/gammaDraw;
phi_Hdraw(1,kk)=phiDrawkk;
end
    % sample phi

%diff = (log(Hvars(p+1:end,:)) - [h0'; log(Hvars(p+1:end-1,:))]).^2;
%phi_Hdraw = 1./gamrnd(nu_h + size(diff,1)/2, 1./(S_h + sum(diff)'/2))';
%phi_Hdraw = 1./gamrnd(S_h + size(diff,1)/2, 1./((nu_h.*S_h) + sum(diff)'/2))';

%
% % sample the initial value
Kh0 = a0 + sparse(1:n,1:n,1./phi_Hdraw');
h0_hat = Kh0\(a0*b0 + log(Hvars(p+1,:)')./phi_Hdraw');
h0 = h0_hat + chol(Kh0,'lower')'\randn(n,1);

    %h0 = priorValues.mean_ln_h0; %initialize initial variance
    %Kh0 = diag(priorValues.var_ln_h0); %initialize initial variance

%% Sample A
[Adraw]=bear.sampleA_H(yData,Psi,B,Hvars,T,priorValues,dataValues);
AdrawInv=Adraw\eye(M);

%% Construct H
Hdraw=zeros(M,M,T);

for t=1:T
    Hdraw(:,:,t)=AdrawInv*diag(Hvars(t,:))*AdrawInv';
end

HvarsDraw=Hvars;

end


text:
function res = rcnn_all(task, model_typ, trainset, testset)
  global RUNNAME
  RUNNAME = 'release';
  p = get_paths(RUNNAME);
  c = benchmarkPaths();
  NYU_ROOT_DIR = c.dataDir;
  REGIONDIR = fullfile(p.box_dir, 'release');
  SALT = 'release';
  MAX_BOXES = 2000;

  imdb = imdb_from_nyud2(NYU_ROOT_DIR, trainset, task, REGIONDIR, SALT, MAX_BOXES);
  imdb.roidb_func = @roidb_from_nyud2;

  feat_opts = struct('featDir', '');
  feat_opts(1) = [];
  global RCNN_CONFIG_OVERRIDE;
  switch model_typ
    case 'rgb_hha',
      conf_override.sub_dir = sprintf('rgb_hha_%d_%s', 30000, trainset);
      feat_opts(end+1) = struct('featDir', fullfile(p.cnnF_cache_dir, 'rgb_30000', imdb.dataset_name));
      feat_opts(end+1) = struct('featDir', fullfile(p.cnnF_cache_dir, 'hha_30000', imdb.dataset_name));
    case 'rgb',
      conf_override.sub_dir = sprintf('rgb_%d_%s', 30000, trainset);
      feat_opts(end+1) = struct('featDir', fullfile(p.cnnF_cache_dir, 'rgb_30000', imdb.dataset_name));
    case 'hha',
      conf_override.sub_dir = sprintf('hha_%d_%s', 30000, trainset);
      feat_opts(end+1) = struct('featDir', fullfile(p.cnnF_cache_dir, 'hha_30000', imdb.dataset_name));
    case 'disparity',
      conf_override.sub_dir = sprintf('disparity_%d_%s', 30000, trainset);
      feat_opts(end+1) = struct('featDir', fullfile(p.cnnF_cache_dir, 'disparity_30000', imdb.dataset_name));
  end
  RCNN_CONFIG_OVERRIDE = @() conf_override;
  rcnn_train_opts = {'feat_opts', feat_opts};

  rcnn_models = rcnn_train(imdb, rcnn_train_opts{:});

  if(isstr(testset)), testset = {testset}; end
  for i = 1:length(testset),
    imdb = imdb_from_nyud2(NYU_ROOT_DIR, testset{i}, task, REGIONDIR, SALT, MAX_BOXES);
    imdb.roidb_func = @roidb_from_nyud2;
    res{i} = rcnn_test(rcnn_models, imdb, '');
  end
end

author: s-gupta

repo: rcnn-depth

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


text:
function sVF = approximation(v, y, varargin)
%
% Syntax
%   sVF = S2VectorField.quadrature(v, value)
%   sVF = S2VectorField.quadrature(v, value, 'bandwidth', bw)
%
% Input
%   value - @vector3d
%   v - @vector3d
%
% Output
%   sVF - @S2VectorFieldHarmonic
%
% Options
%   bw - degree of the spherical harmonic (default: 128)
%

y = y.xyz;
sF = S2FunHarmonic.quadrature(v, y, varargin{:});

sVF = S2VectorFieldHarmonic(sF);

end

author: mtex-toolbox

repo: mtex

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


text:
%MDL_PHANTOMX Create model of PhantomX pincher manipulator
%
% MDL_PHANTOMX is a script that creates the workspace variable px which
% describes the kinematic characteristics of a PhantomX Pincher Robot, a 4
% joint hobby class  manipulator by Trossen Robotics.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%
% Notes::
% - Uses standard DH conventions.
% - Tool centrepoint is middle of the fingertips.
% - All translational units in mm.
%
% Reference::
%
% - http://www.trossenrobotics.com/productdocs/assemblyguides/phantomx-basic-robot-arm.html

% MODEL: Trossen Robotics, PhantomX Pincher, 4DOF, standard_DH

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

clear L
L(1) = Revolute('d', 40, 'alpha', -pi/2);
L(2) = Revolute('a', -105, 'alpha', pi, 'offset', pi/2);
L(3) = Revolute('a', -105);
L(4) = Revolute('a', -105);

% Note alpha_2 = pi, needed to account for rotation axes of joints 3 and 4 having
% opposite sign to joint 2.
%
% s='Rz(q1) Tz(L1) Ry(q2) Tz(L2) Ry(q3) Tz(L3) Ry(q4) Tz(L4)'
% DHFactor(s)

px = SerialLink(L, 'name', 'PhantomX', 'manufacturer', 'Trossen Robotics');
qz = [0 0 0 0];
px.tool = trotz(-pi/2) * trotx(pi/2);

author: Allopart

repo: rbpf-gmapping


text:
function setfotnsizefigure(fontsize)
set(get(gca,'XLabel'),'FontSize',fontsize)
set(get(gca,'YLabel'),'FontSize',fontsize)
set(gca,'FontSize',fontsize)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

for i=1:size(X, 1)
    distances = zeros(size(centroids, 1), 1);
    for k=1:K
        distances(k) = sumsq(X(i,:) - centroids(k, :));
    end;
    [minDistance, minIndex] = min(distances);
    idx(i) = minIndex;
end;

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

end

author: zsiciarz

repo: ml-coursera


text:
function [cost grad] = soft_cost(theta, instances, labels, lambda)

[s1 s2] = size(instances);
[~, s4] = size(labels);

mat = reshape(theta(1:end-s4), s4, s2)';
b = theta(end-s4+1:end);
vec = mat'*instances';
if s4 > 1
    sm = softmax(vec + b(:,ones(1,s1)));
else
    sm = sigmoid(vec + b(:,ones(1,s1)));
end
lbl = labels';
lbl_sm = sm - lbl;

cost = 1/2*sum(sum(((lbl_sm').*(lbl_sm'))))/s1 + 1/2*lambda*(theta(1:end-s4)'*theta(1:end-s4));

if s4 > 1
    del = (sm.*lbl_sm) - bsxfun(@times,sum(sm.*lbl_sm), sm);
    gradW = del*instances;
    gradb = sum(del,2);
else
    del = (lbl_sm).*sigmoid_prime(sm);
    gradW = del*instances;
    gradb = sum(del);

end
grad = [gradW(:); gradb]/s1 + [lambda*theta(1:end-s4); zeros(s4,1)];

author: jacoxu


text:
function [ U, S, V ] = SubspaceIter( U1, V1, U0, V0, spa, bi, R, maxIter )

if(size(U1, 1) > size(V1, 1))
    V1 = V1*(1 + bi);
else
    U1 = U1*(1 + bi);
end

if(size(U0, 1) > size(V0, 1))
    V0 = V0*(-bi);
else
    U0 = U0*(-bi);
end

V = R;
for i = 1:maxIter
    U = AfuncAcc( U1, V1, U0, V0, spa, V);
    [U, ~] = qr(U, 0);

    V = AtfuncAcc( U1, V1, U0, V0, spa, U);
    [V, ~] = qr(V, 0);
end

S = AfuncAcc( U1, V1, U0, V0, spa, V);
S = U'*S;

[Us, S, Vs] = svd(S, 'econ');

U = U*Us;
V = V*Vs;

end

author: HKUST-KnowComp

repo: FMG

save_path: github-repos/MATLAB/HKUST-KnowComp-FMG