text:
% % See also plotbezierpatch3D.m

% % function to plot 3D cubic Bezier surface in many ways.

% % A Bezier surface is composed of one or more patches.
% % A patch is defined by 16 control points arranged in 4 x 4 x 3 matrix.
% % A control point has three coordinates (x,y,z)

% % Details
% % -> Input Matrix S stores all the control points of all the patches of
% %    a Bezier surface such that
% % -> S(:,:,:,k) holds control points of kth patch,
% %    where k=1..np and np is number of patches in the surface
% % -> Size of S(:,:,:,k) is 4 x 4 x 3, i.e., 16 control points and each
% %    control point has three coordinates (x,y,z)
% %    S(:,:,1,k): x-coordates of control points of kth patch as 4 x 4 matrix
% %    S(:,:,2,k): y-coordates of control points of kth patch as 4 x 4 matrix
% %    S(:,:,3,k): z-coordates of control points of kth patch as 4 x 4 matrix
% % -> Input Matrix Q stores interpolated values between control points
% %    Q is similar to S in format but has more values, i.e., it stores
% %    end control points and interpolated values.
% %    see the function bezierpatchinterp.m for more details

function plotbeziersurface3D(S,Q)

[r c dim np]=size(S);
% % np: number of patches

if dim > 3 or dim < 3
    error ('plotbeziersurface3D.m function handles only 3-D points')
end

az=21;  %azimuth
el=19;  %elevation.

lw=1; %plotting linewidth

str1='\bf Control Point';
str2='\bf Control Polygon';
str3='\bf Surface (bi-directional Bezier curve)';

% %-----------------------------------------------
% % For plotting it is convenient if we reshape the data
% % into vecotor format and separate (X,Y,Z) coordinates

xS=[]; yS=[]; zS=[];
for k=1:np
     xS =horzcat(xS, reshape(S(:,:,1,k),1,[]));
     yS =horzcat(yS, reshape(S(:,:,2,k),1,[]));
     zS =horzcat(zS, reshape(S(:,:,3,k),1,[]));
end

xQ=[]; yQ=[]; zQ=[];
for k=1:np
     xQ =horzcat(xQ, reshape(Q(:,:,1,k),1,[]));
     yQ =horzcat(yQ, reshape(Q(:,:,2,k),1,[]));
     zQ =horzcat(zQ, reshape(Q(:,:,3,k),1,[]));
end
% %------------------------------------------------
% % Wireframe Plot of Bezier Surface using interpolation of control points
% % Control Points and Control Polygon are marked
figure, hold on
plot3(xS,yS,zS,'ro','LineWidth',lw)
plot3(xS,yS,zS,'g','LineWidth',lw)
plot3(xQ,yQ,zQ,'b','LineWidth',lw)
legend(str1,str2,str3);
title('\bf Wireframe Plot of a Bezier Surface with Control Points and Control Polygon')
view(3); box;  view(az,el)
% %-----------------------------------------------
% % Wireframe Plot of Bezier Surface using interpolation of control points
figure, hold on
plot3(xQ,yQ,zQ,'b','LineWidth',lw)
legend(str3);
title('\bf Wireframe Plot of a Bezier Surface')
view(3); box;  view(az,el)
% %-----------------------------------------------
% % Plot of Surface using control points
figure, hold on
for k=1:np
    surface(S(:,:,1,k),S(:,:,2,k),S(:,:,3,k),'FaceColor','green')
end
title('\bf Bezier Surface using Control Points');
view(3); box;  view(az,el)
% %-----------------------------------------------
% % Plot of Surface using interpolation of control points
figure, hold on
for k=1:np
    surface(Q(:,:,1,k),Q(:,:,2,k),Q(:,:,3,k),'FaceColor','green')
end
title('\bf Bezier Surface using Interpolated Points');
view(3); box;  view(az,el)

% %-----------------------------------------------
% % Surface with interpolated shading
figure, hold on
for k=1:np
    surface(Q(:,:,1,k),Q(:,:,2,k),Q(:,:,3,k))
end
shading interp
title('\bf Bezier Surface using Interpolated Shading');
view(3); box;  view(az,el)
% %-----------------------------------------------%
% % Surface with faceted shading
figure, hold on
for k=1:np
    surface(Q(:,:,1,k),Q(:,:,2,k),Q(:,:,3,k))
end
shading faceted
title('\bf Bezier Surface using Faceted Shading');
view(3); box;  view(az,el)
% %-----------------------------------------------%

% % --------------------------------
% % This program or any other program(s) supplied with it does not provide any
% % warranty direct or implied.
% % This program is free to use/share for non-commerical purpose only.
% % Kindly reference the author.
% % Author: Dr. Murtaza Khan
% % URL : http://www.linkedin.com/pub/dr-murtaza-khan/19/680/3b3
% % --------------------------------


text:
function [stat, cfg] = ft_statistics_montecarlo(cfg, dat, design, varargin)

% FT_STATISTICS_MONTECARLO performs a nonparametric statistical test by calculating
% Monte-Carlo estimates of the significance probabilities and/or critical values from
% the permutation distribution. This function should not be called directly, instead
% you should call the function that is associated with the type of data on which you
% want to perform the test.
%
% Use as
%   stat = ft_timelockstatistics(cfg, data1, data2, data3, ...)
%   stat = ft_freqstatistics    (cfg, data1, data2, data3, ...)
%   stat = ft_sourcestatistics  (cfg, data1, data2, data3, ...)
%
% where the data is obtained from FT_TIMELOCKANALYSIS, FT_FREQANALYSIS or
% FT_SOURCEANALYSIS respectively, or from FT_TIMELOCKGRANDAVERAGE,
% FT_FREQGRANDAVERAGE or FT_SOURCEGRANDAVERAGE respectively
% and with cfg.method = 'montecarlo'
%
% The configuration options that can be specified are:
%   cfg.numrandomization = number of randomizations, can be 'all'
%   cfg.correctm         = string, apply multiple-comparison correction, 'no', 'max', cluster', 'tfce', 'bonferroni', 'holm', 'hochberg', 'fdr' (default = 'no')
%   cfg.alpha            = number, critical value for rejecting the null-hypothesis per tail (default = 0.05)
%   cfg.tail             = number, -1, 1 or 0 (default = 0)
%   cfg.correcttail      = string, correct p-values or alpha-values when doing a two-sided test, 'alpha','prob' or 'no' (default = 'no')
%   cfg.ivar             = number or list with indices, independent variable(s)
%   cfg.uvar             = number or list with indices, unit variable(s)
%   cfg.wvar             = number or list with indices, within-cell variable(s)
%   cfg.cvar             = number or list with indices, control variable(s)
%   cfg.feedback         = string, 'gui', 'text', 'textbar' or 'no' (default = 'text')
%   cfg.randomseed       = string, 'yes', 'no' or a number (default = 'yes')
%
% If you use a cluster-based statistic, you can specify the following options that
% determine how the single-sample or single-voxel statistics will be thresholded and
% combined into one statistical value per cluster.
%   cfg.clusterstatistic = how to combine the single samples that belong to a cluster, 'maxsum', 'maxsize', 'wcm' (default = 'maxsum')
%                          the option 'wcm' refers to 'weighted cluster mass', a statistic that combines cluster size and intensity;
%                          see Hayasaka & Nichols (2004) NeuroImage for details
%   cfg.clusterthreshold = method for single-sample threshold, 'parametric', 'nonparametric_individual', 'nonparametric_common' (default = 'parametric')
%   cfg.clusteralpha     = for either parametric or nonparametric thresholding per tail (default = 0.05)
%   cfg.clustercritval   = for parametric thresholding (default is determined by the statfun)
%   cfg.clustertail      = -1, 1 or 0 (default = 0)
%
% To include the channel dimension for clustering of channel level data, you should specify
%   cfg.neighbours       = neighbourhood structure, see FT_PREPARE_NEIGHBOURS
% If you specify an empty neighbourhood structure, clustering will only be done
% over frequency and/or time and not over neighbouring channels.
%
% The statistic that is computed for each sample in each random reshuffling
% of the data is specified as
%   cfg.statistic       = 'indepsamplesT'           independent samples T-statistic,
%                         'indepsamplesF'           independent samples F-statistic,
%                         'indepsamplesregrT'       independent samples regression coefficient T-statistic,
%                         'indepsamplesZcoh'        independent samples Z-statistic for coherence,
%                         'depsamplesT'             dependent samples T-statistic,
%                         'depsamplesFmultivariate' dependent samples F-statistic MANOVA,
%                         'depsamplesregrT'         dependent samples regression coefficient T-statistic,
%                         'actvsblT'                activation versus baseline T-statistic.
% or you can specify your own low-level statistical function.
%
% You can also use a custom statistic of your choice that is sensitive to the
% expected effect in the data. You can implement the statistic in a "statfun" that
% will be called for each randomization. The requirements on a custom statistical
% function is that the function is called ft_statfun_xxx, and that the function returns
% a structure with a "stat" field containing the single sample statistical values.
% Have a look at the functions in the fieldtrip/statfun directory (e.g.
% FT_STATFUN_INDEPSAMPLEST) for the correct format of the input and output.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS, FT_SOURCESTATISTICS,
% FT_STATISTICS_ANALYTIC, FT_STATISTICS_STATS, FT_STATISTICS_MVPA,
% FT_STATISTICS_CROSSVALIDATE

% Undocumented local options:
%   cfg.resampling       permutation, bootstrap
%   cfg.computecritval   yes|no, for the statfun
%   cfg.computestat      yes|no, for the statfun
%   cfg.computeprob      yes|no, for the statfun
%   cfg.voxelstatistic   deprecated
%   cfg.voxelthreshold   deprecated
%   cfg.precondition     before|after|[], for the statfun

% Copyright (C) 2005-2015, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% do a sanity check on the input data
assert(isnumeric(dat),    'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');
assert(isnumeric(design), 'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',     {'factor',           'ivar'});
cfg = ft_checkconfig(cfg, 'renamed',     {'unitfactor',       'uvar'});
cfg = ft_checkconfig(cfg, 'renamed',     {'repeatedmeasures', 'uvar'});
cfg = ft_checkconfig(cfg, 'renamedval',  {'clusterthreshold', 'nonparametric', 'nonparametric_individual'});
cfg = ft_checkconfig(cfg, 'renamedval',  {'correctm', 'yes', 'max'});
cfg = ft_checkconfig(cfg, 'renamedval',  {'correctm', 'none', 'no'});
cfg = ft_checkconfig(cfg, 'renamedval',  {'correctm', 'bonferoni', 'bonferroni'});
cfg = ft_checkconfig(cfg, 'renamedval',  {'correctm', 'holms', 'holm'});
cfg = ft_checkconfig(cfg, 'required',    {'statistic'});
cfg = ft_checkconfig(cfg, 'forbidden',   {'ztransform', 'removemarginalmeans', 'randomfactor', 'voxelthreshold', 'voxelstatistic'});
cfg = ft_checkconfig(cfg, 'renamedval',  {'statfun', 'depsamplesF', 'ft_statfun_depsamplesFmultivariate'});
cfg = ft_checkconfig(cfg, 'renamedval',  {'statfun', 'ft_statfun_depsamplesF', 'ft_statfun_depsamplesFmultivariate'});

% set the defaults for the main function
cfg.alpha        = ft_getopt(cfg, 'alpha',      0.05);
cfg.tail         = ft_getopt(cfg, 'tail',       0);
cfg.correctm     = ft_getopt(cfg, 'correctm',   'no');
cfg.resampling   = ft_getopt(cfg, 'resampling', 'permutation');
cfg.feedback     = ft_getopt(cfg, 'feedback',   'text');
cfg.ivar         = ft_getopt(cfg, 'ivar',       'all');
cfg.uvar         = ft_getopt(cfg, 'uvar',       []);
cfg.cvar         = ft_getopt(cfg, 'cvar',       []);
cfg.wvar         = ft_getopt(cfg, 'wvar',       []);
cfg.correcttail  = ft_getopt(cfg, 'correcttail',  'no');
cfg.precondition = ft_getopt(cfg, 'precondition', []);

% explicit check for option 'yes' in cfg.correctail.
if strcmp(cfg.correcttail, 'yes')
  ft_error('cfg.correcttail = ''yes'' is not allowed, use either ''prob'', ''alpha'' or ''no''')
end

if strcmp(cfg.correctm, 'tfce')
  % TODO this could require some better defaults
  cfg.connectivity = ft_getopt(cfg, 'connectivity', []);
  cfg.tfce_h0      = ft_getopt(cfg, 'tfce_h0', 0);
  cfg.tfce_H       = ft_getopt(cfg, 'tfce_H',  2);
  cfg.tfce_E       = ft_getopt(cfg, 'tfce_E',  0.5);
  cfg.tfce_nsteps  = ft_getopt(cfg, 'tfce_nsteps', 100);
else
  % these options only apply to tfce, to ensure appropriate configs they are forbidden when _not_ clustering
  cfg = ft_checkconfig(cfg, 'unused', {'tfce_h0', 'tfce_H', 'tfce_E', 'tfce_nsteps'});
end

if strcmp(cfg.correctm, 'cluster')
  % set the defaults for clustering
  cfg.clusterstatistic = ft_getopt(cfg, 'clusterstatistic', 'maxsum');
  cfg.clusterthreshold = ft_getopt(cfg, 'clusterthreshold', 'parametric');
  cfg.clusteralpha     = ft_getopt(cfg, 'clusteralpha',     0.05);
  cfg.clustercritval   = ft_getopt(cfg, 'clustercritval',   []);
  cfg.clustertail      = ft_getopt(cfg, 'clustertail',      cfg.tail);
  cfg.connectivity     = ft_getopt(cfg, 'connectivity',     []); % the default is dealt with below
else
  % these options only apply to clustering, to ensure appropriate configs they are forbidden when _not_ clustering
  cfg = ft_checkconfig(cfg, 'unused', {'clusterstatistic', 'clusteralpha', 'clustercritval', 'clusterthreshold', 'clustertail'});
end

if any(strcmp(cfg.correctm, {'cluster' 'tfce'}))
  % these options might require a spatial neighbourhood matrix

  % deal with the neighbourhood of the channels/triangulation/voxels
  if isempty(cfg.connectivity)
    if isfield(cfg, 'dim') && ~isfield(cfg, 'channel') && ~isfield(cfg, 'tri')
      % input data can be reshaped into a 3D volume, use bwlabeln/spm_bwlabel rather than clusterstat
      ft_info('using connectivity of voxels in 3-D volume\n');
      cfg.connectivity = nan;
    elseif isfield(cfg, 'tri')
      % input data describes a surface along which neighbours can be defined
      ft_info('using connectivity of vertices along triangulated surface\n');
      cfg.connectivity = triangle2connectivity(cfg.tri);
      if isfield(cfg, 'insideorig')
        cfg.connectivity = cfg.connectivity(cfg.insideorig, cfg.insideorig);
      end
    elseif isfield(cfg, 'avgoverchan') && istrue(cfg.avgoverchan)
      % channel dimension has been averaged across, no sense in clustering across space
      cfg.connectivity = true(1);
    elseif isfield(cfg, 'channel')
      cfg.neighbours   = ft_getopt(cfg, 'neighbours', []);
      cfg.connectivity = channelconnectivity(cfg);
    else
      % there is no connectivity in the spatial dimension
      cfg.connectivity = false(size(dat,1));
    end
  else
    % use the specified connectivity: this is not fully robust because
    % there is no guarantee that the order of the spatial elements in the
    % data is the same as the order of the spatial elements in the
    % adjacency matrix
  end
end

% for backward compatibility and other warnings relating correcttail
if isfield(cfg,'correctp') && strcmp(cfg.correctp,'yes')
  ft_warning('cfg.correctp has been renamed to cfg.correcttail and the options have been changed')
  disp('setting cfg.correcttail to ''prob''')
  cfg.correcttail = 'prob';
  cfg = rmfield(cfg,'correctp');
elseif isfield(cfg,'correctp') && strcmp(cfg.correctp,'no')
  cfg = ft_checkconfig(cfg, 'renamed', {'correctp', 'correcttail'});
end
if strcmp(cfg.correcttail,'no') && cfg.tail==0 && cfg.alpha==0.05
  ft_warning('Doing a two-sided test without correcting p-values or alpha-level, p-values and alpha-level will reflect one-sided tests per tail. See http://bit.ly/2YQ1Hm8')
end

% for backward compatibility
if size(design,2)~=size(dat,2)
  design = transpose(design);
end

if ischar(cfg.ivar) && strcmp(cfg.ivar, 'all')
  cfg.ivar = 1:size(design,1);
end

% fetch function handle to the low-level statistics function
statfun = ft_getuserfun(cfg.statistic, 'statfun');
if isempty(statfun)
  ft_error('could not locate the appropriate statistics function');
else
  ft_info('using "%s" for the single-sample statistics\n', func2str(statfun));
end

% construct the resampled design matrix or data-shuffling matrix
ft_info('constructing randomized design\n');
resample = resampledesign(cfg, design);
Nrand = size(resample,1);

% most of the statfuns result in this warning, which is not interesting
ws = ft_warning('off', 'MATLAB:warn_r14_stucture_assignment');

if strcmp(cfg.correctm, 'cluster')
  % determine the critical value for cluster thresholding
  if strcmp(cfg.clusterthreshold, 'nonparametric_individual') || strcmp(cfg.clusterthreshold, 'nonparametric_common')
    ft_info('using a nonparametric threshold for clustering\n');
    cfg.clustercritval = [];  % this will be determined later
  elseif strcmp(cfg.clusterthreshold, 'parametric') && isempty(cfg.clustercritval)
    ft_info('computing a parametric threshold for clustering\n');
    tmpcfg = cfg; % the next line does not pass on non-standard options that a statfun might use
    % tmpcfg = keepfields(cfg, {'dim' 'dimord' 'clusteralpha' 'clustertail' 'ivar' 'uvar' 'cvar' 'wvar' 'contrastcoefs'});
    tmpcfg.computecritval = 'yes';  % explicitly request the computation of the crtitical value
    tmpcfg.computestat    = 'no';   % skip the computation of the statistic
    tmpcfg.alpha          = cfg.clusteralpha; % the statfun uses cfg.alpha most likely
    try
      cfg.clustercritval    = getfield(statfun(tmpcfg, dat, design), 'critval');
    catch
      disp(lasterr);
      ft_error('could not determine the parametric critical value for clustering');
    end
  elseif strcmp(cfg.clusterthreshold, 'parametric') && ~isempty(cfg.clustercritval)
    ft_info('using the specified parametric threshold for clustering\n');
    cfg.clusteralpha = [];
  end
end

% compute the statistic for the observed data
ft_progress('init', cfg.feedback, 'computing statistic');
% get an estimate of the time required per evaluation of the statfun
time_pre = cputime;

try
  % the nargout function in MATLAB 6.5 and older does not work on function handles
  num = nargout(statfun);
catch
  num = 1;
end

if num==1
  % only the statistic is returned
  [statobs] = statfun(cfg, dat, design);
elseif num==2
  % both the statistic and the (updated) configuration are returned
  [statobs, cfg] = statfun(cfg, dat, design);
elseif num==3
  % both the statistic and the (updated) configuration and the (updated) data are returned
  tmpcfg = cfg;
  if strcmp(cfg. precondition, 'before'), tmpcfg.preconditionflag = 1; end
  [statobs, tmpcfg, dat]  = statfun(tmpcfg, dat, design);
  tmpcfg.preconditionflag = 0;
  cfg = tmpcfg;
end

if isstruct(statobs)
  % remember all details for later reference, continue to work with the statistic
  statfull = statobs;
  statobs  = statobs.stat;
end

% remember the statistic for later reference, continue to work with the statistic
statfull.stat = statobs;

time_eval = cputime - time_pre;
ft_info('estimated time per randomization is %.2f seconds\n', time_eval);

% pre-allocate some memory
if strcmp(cfg.correctm, 'cluster')
  statrand = zeros(size(statobs,1), size(resample,1), class(dat)); % this reduces the memory footprint, requires the user to use ft_struct2single on the input data
else
  prb_pos   = zeros(size(statobs));
  prb_neg   = zeros(size(statobs));
end

if strcmp(cfg.precondition, 'after')
  tmpcfg = cfg;
  tmpcfg.preconditionflag = 1;
  [tmpstat, tmpcfg, dat] = statfun(tmpcfg, dat, design);
end

if any(strcmp(cfg.correctm, {'tfce' 'max'}))
  % pre-allocate the memory to hold the distribution of most extreme positive (right) and negative (left) statistical values
  posdistribution = nan(1,Nrand);
  negdistribution = nan(1,Nrand);
end

% compute the statistic for the randomized data and count the outliers
for i = 1:Nrand
  ft_progress(i/Nrand, 'computing statistic %d from %d\n', i, Nrand);
  if strcmp(cfg.resampling, 'permutation')
    tmpdesign = design(:,resample(i,:));     % the columns in the design matrix are reshufled by means of permutation
    tmpdat    = dat;                        % the data itself is not shuffled
    if size(tmpdesign,1)==size(tmpdat,2)
      tmpdesign = transpose(tmpdesign);
    end
  elseif strcmp(cfg.resampling, 'bootstrap')
    tmpdesign = design;                     % the design matrix is not shuffled
    tmpdat    = dat(:,resample(i,:));        % the columns of the data are resampled by means of bootstrapping
  end
  if any(strcmp(cfg.correctm, {'cluster' 'tfce'}))
    % keep each randomization in memory for cluster postprocessing
    dum = statfun(cfg, tmpdat, tmpdesign);
    if isstruct(dum)
      statrand(:,i) = dum.stat;
    else
      statrand(:,i) = dum;
    end
  else
    % do not keep each randomization in memory, but process them on the fly
    statrand = statfun(cfg, tmpdat, tmpdesign);
    if isstruct(statrand)
      statrand = statrand.stat;
    end

    % the following line is for debugging
    % stat.statkeep(:,i) = statrand;
    if strcmp(cfg.correctm, 'max')
      % compare each data element with the maximum statistic
      prb_pos = prb_pos + (statobs<max(statrand(:)));
      prb_neg = prb_neg + (statobs>min(statrand(:)));
      posdistribution(i) = max(statrand(:));
      negdistribution(i) = min(statrand(:));
    else
      % compare each data element with its own statistic
      prb_pos = prb_pos + (statobs<statrand);
      prb_neg = prb_neg + (statobs>statrand);
    end
  end
end
ft_progress('close');

if strcmp(cfg.correctm, 'cluster')
  % do the cluster postprocessing
  [stat, cfg] = clusterstat(cfg, statrand, statobs);
elseif strcmp(cfg.correctm, 'tfce')
  [stat, cfg] = tfcestat(cfg, statrand, statobs);
else
  if ~isequal(cfg.numrandomization, 'all')
    % in case of random permutations (i.e., montecarlo sample, and NOT full
    % permutation), the minimum p-value should not be 0, but 1/N
    prb_pos = prb_pos + 1;
    prb_neg = prb_neg + 1;
    Nrand = Nrand + 1;
  end
  switch cfg.tail
    case 1
      clear prb_neg  % not needed any more, free some memory
      stat.prob = prb_pos./Nrand;
    case -1
      clear prb_pos  % not needed any more, free some memory
      stat.prob = prb_neg./Nrand;
    case 0
      % for each observation select the tail that corresponds with the lowest probability
      prb_neg = prb_neg./Nrand;
      prb_pos = prb_pos./Nrand;
      stat.prob = min(prb_neg, prb_pos); % this is the probability for the most unlikely tail
  end
end

% In case of a two tailed test, the type-I error rate (alpha) refers to
% both tails of the distribution, whereas the stat.prob value computed sofar
% corresponds with one tail, i.e. the probability, under the assumption of
% no effect or no difference (the null hypothesis), of obtaining a result
% equal to or more extreme than what was actually observed. The decision
% rule whether the null-hopothesis should be rejected given the observed
% probability therefore should consider alpha divided by two, to correspond
% with the probability in one of the tails (the most extreme tail). This
% is conceptually equivalent to performing a Bonferroni correction for the
% two tails.
%
% An alternative solution to distribute the alpha level over both tails is
% achieved by multiplying the probability with a factor of two, prior to
% thresholding it wich cfg.alpha.  The advantage of this solution is that
% it results in a p-value that corresponds with a parametric probability.
% Below both options are realized
if strcmp(cfg.correcttail, 'prob') && cfg.tail==0
  stat.prob = stat.prob .* 2;
  stat.prob(stat.prob>1) = 1; % clip at p=1
  % also correct the probabilities in the pos/negcluster fields
  if isfield(stat, 'posclusters')
    for i=1:length(stat.posclusters)
      stat.posclusters(i).prob = stat.posclusters(i).prob*2;
      if stat.posclusters(i).prob>1; stat.posclusters(i).prob = 1; end
    end
  end
  if isfield(stat, 'negclusters')
    for i=1:length(stat.negclusters)
      stat.negclusters(i).prob = stat.negclusters(i).prob*2;
      if stat.negclusters(i).prob>1; stat.negclusters(i).prob = 1; end
    end
  end
elseif strcmp(cfg.correcttail, 'alpha') && cfg.tail==0
  cfg.alpha = cfg.alpha / 2;
end

% compute range of confidence interval p ? 1.96(sqrt(var(p))), with var(p) = var(x/n) = p*(1-p)/N
stddev = sqrt(stat.prob.*(1-stat.prob)/Nrand);
stat.cirange = 1.96*stddev;

if isfield(stat, 'posclusters')
  for i=1:length(stat.posclusters)
    stat.posclusters(i).stddev  = sqrt(stat.posclusters(i).prob.*(1-stat.posclusters(i).prob)/Nrand);
    stat.posclusters(i).cirange =  1.96*stat.posclusters(i).stddev;
    if i==1 && stat.posclusters(i).prob<cfg.alpha && stat.posclusters(i).prob+stat.posclusters(i).cirange>=cfg.alpha
      ft_warning('FieldTrip:posCluster_exceeds_alpha', sprintf('The p-value confidence interval of positive cluster #%i includes %.3f - consider increasing the number of permutations!', i, cfg.alpha));
    end
  end
end
if isfield(stat, 'negclusters')
  for i=1:length(stat.negclusters)
    stat.negclusters(i).stddev  = sqrt(stat.negclusters(i).prob.*(1-stat.negclusters(i).prob)/Nrand);
    stat.negclusters(i).cirange =  1.96*stat.negclusters(i).stddev;
    if i==1 && stat.negclusters(i).prob<cfg.alpha && stat.negclusters(i).prob+stat.negclusters(i).cirange>=cfg.alpha
      ft_warning('FieldTrip:negCluster_exceeds_alpha', sprintf('The p-value confidence interval of negative cluster #%i includes %.3f - consider increasing the number of permutations!', i, cfg.alpha));
    end
  end
end

if ~isfield(stat, 'prob')
  ft_warning('probability was not computed');
else
  switch lower(cfg.correctm)
    case 'max'
      % the correction is implicit in the method
      ft_notice('using a maximum-statistic based method for multiple comparison correction\n');
      ft_notice('the returned probabilities and the thresholded mask are corrected for multiple comparisons\n');
      stat.mask = stat.prob<=cfg.alpha;
      stat.posdistribution = posdistribution;
      stat.negdistribution = negdistribution;
    case 'tfce'
      ft_notice('using a threshold free cluster enhancement based method for multiple comparison correction\n');
      ft_notice('the returned probabilities and the thresholded mask are corrected for multiple comparisons\n');
      stat.mask = stat.prob<=cfg.alpha;
    case 'cluster'
      % the correction is implicit in the method
      ft_notice('using a cluster-based method for multiple comparison correction\n');
      ft_notice('the returned probabilities and the thresholded mask are corrected for multiple comparisons\n');
      stat.mask = stat.prob<=cfg.alpha;
    case 'bonferroni'
      ft_notice('performing Bonferroni correction for multiple comparisons\n');
      ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
      stat.mask = stat.prob<=(cfg.alpha ./ numel(stat.prob));
    case 'holm'
      % test the most significatt significance probability against alpha/N, the second largest against alpha/(N-1), etc.
      ft_notice('performing Holm-Bonferroni correction for multiple comparisons\n');
      ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
      [pvals, indx] = sort(stat.prob(:));                                   % this sorts the significance probabilities from smallest to largest
      k = find(pvals > (cfg.alpha ./ ((length(pvals):-1:1)')), 1, 'first'); % compare each significance probability against its individual threshold
      mask = (1:length(pvals))'<k;
      stat.mask = zeros(size(stat.prob));
      stat.mask(indx) = mask;
    case 'hochberg'
      % test the most significant significance probability against alpha/N, the second largest against alpha/(N-1), etc.
      ft_notice('performing Hochberg''s correction for multiple comparisons (this is *not* the Benjamini-Hochberg FDR procedure!)\n');
      ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
      [pvals, indx] = sort(stat.prob(:));                     % this sorts the significance probabilities from smallest to largest
      k = find(pvals <= (cfg.alpha ./ ((length(pvals):-1:1)')), 1, 'last'); % compare each significance probability against its individual threshold
      mask = (1:length(pvals))'<=k;
      stat.mask = zeros(size(stat.prob));
      stat.mask(indx) = mask;
    case 'fdr'
      ft_notice('performing FDR correction for multiple comparisons\n');
      ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
      stat.mask = fdr(stat.prob, cfg.alpha);
    otherwise
      ft_notice('not performing a correction for multiple comparisons\n');
      stat.mask = stat.prob<=cfg.alpha;
  end
end

% return the observed statistic
if ~isfield(stat, 'stat')
  stat.stat = statobs;
end

if exist('statrand', 'var')
  stat.ref = mean(statrand,2);
end

% return optional other details that were returned by the statfun
stat = copyfields(statfull, stat, fieldnames(statfull));

ft_warning(ws); % revert to original state

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear all;
%%%====== Settings ======%%%
model = 'Multi-purpose CNN'; % 'Deep SR-ITM' (ICCV'19) or 'Multi-purpose CNN' (ACCV'18)
SDR_file = './data/test/testset_SDR.mat'; % input .mat file
HDR_file = './data/test/testset_HDR.mat'; % GT .mat file
scale = 2; % scale factor for SR
pred_file = sprintf('./pred/pred_x%d.mat', scale); % result .mat file
%%%======================%%%
addpath('utils');
% load data
disp(['Testing for scale ', num2str(scale), '...'])
disp('Loading file...')
SDR = load(SDR_file);
HDR = load(HDR_file);
data = SDR.SDR;
label = HDR.HDR;

% initialize
psnr_all = zeros(1, size(data, 4));
ssim_all = zeros(1, size(data, 4));
mpsnr_all = zeros(1, size(data, 4));
msssim_all = zeros(1, size(data, 4));
pred = single(zeros(size(data)));

% load net
disp('Loading net...')
if strcmp(model, 'Deep SR-ITM')
    netstruct = load(sprintf('./net/x%d.mat', scale));
elseif strcmp(model, 'Multi-purpose CNN')
    netstruct = load(sprintf('./net/Multi-purpose_CNN_x%d.mat', scale));
end
if strcmp(model, 'Deep SR-ITM')
    netstruct.net.layers(172).type='dagnn.SubPixel_Conv_cpu';
elseif strcmp(model, 'Multi-purpose CNN')
    netstruct.net.layers(37).type='dagnn.SubPixel_Conv_cpu';
end
net = dagnn.DagNN.loadobj(netstruct.net);
move(net,'cpu');
net.mode = 'test' ;
pred_index = net.getVarIndex('pred');
net.conserveMemory = true;

% test
disp('Testing starts...')
for fr = 1:size(data, 4)
    % read frames
    SDR_YUV = single(data(:, :, :, fr));
    HDR_YUV = single(label(:, :, :, fr));
    % normalize
    SDR_YUV = SDR_YUV/255;
    HDR_YUV = HDR_YUV/1023;
    % create LR data
    SDR_LR_YUV = imresize(SDR_YUV, 1/scale);
    % prediction
    net.eval({'input', SDR_LR_YUV});
    pred_fr = net.vars(pred_index).value;
    pred(:, :, :, fr) = min(max(pred_fr, 0), 1);

    %%% Evaluation %%% (comment or uncomment appropriate lines)
    % *some metrics may be slow*
    psnr_all(fr) = psnr(HDR_YUV, pred(:, :, :, fr), 1); % PSNR
%     ssim_all(fr) = ssim(HDR_YUV, pred(:, :, :, fr)); % SSIM
%     mpsnr_all(fr) = mPSNR_HDR(HDR_YUV, pred(:, :, :, fr), -3, 3); % mPSNR
%     msssim_all(fr) = msssim(HDR_YUV*1023, pred(:, :, :, fr)*1023, 1023);

    disp(['#', num2str(fr), ' PSNR: ', num2str(psnr_all(fr)), ' dB'])
%     disp(['#', num2str(fr), ' PSNR: ', num2str(psnr_all(fr)), ' dB', ' SSIM: ', num2str(ssim_all(fr)),...
%         ' mPSNR: ', num2str(mpsnr_all(fr)), ' dB', ' MS-SSIM: ', num2str(msssim_all(fr))]);
end
disp(['Average PSNR: ', num2str(mean(psnr_all)), ' dB'])
% disp(['Avg PSNR: ', num2str(mean(psnr_all)), ' dB', ' Avg SSIM: ', num2str(mean(ssim_all)),...
%     ' Avg mPSNR: ', num2str(mean(mpsnr_all)), ' dB', ' Avg MS-SSIM: ', num2str(mean(msssim_all))]);
% REFERENCE x2: Avg PSNR: 35.5966 dB Avg SSIM: 0.9833 Avg mPSNR: 38.0449 dB Avg MS-SSIM: 0.9843
% REFERENCE x4: Avg PSNR: 33.69 dB Avg SSIM: 0.9657 Avg mPSNR: 36.13 Avg MS-SSIM: 0.9754

% save as .mat file
disp('Saving...')
save(pred_file, 'pred', '-v7.3');
disp('Done!')

author: sooyekim

repo: Deep-SR-ITM


text:
%  INTERNAL FUNCTION: doubling_solve solves the linear equation X=A*X*B+C using doubling algorithm
%
%  ::
%
%     [P,retcode]=doubling(A,B,C);
%
%

author: jmaih


text:
classdef circularGraph < handle
% CIRCULARGRAPH Plot an interactive circular graph to illustrate connections in a network.
%
%% Syntax
% circularGraph(X)
% circularGraph(X,'PropertyName',propertyvalue,...)
% h = circularGraph(...)
%
%% Description
% A 'circular graph' is a visualization of a network of nodes and their
% connections. The nodes are laid out along a circle, and the connections
% are drawn within the circle. Click on a node to make the connections that
% emanate from it more visible or less visible. Click on the 'Show All'
% button to make all nodes and their connections visible. Click on the
% 'Hide All' button to make all nodes and their connections less visible.
%
% Required input arguments.
% X : A symmetric matrix of numeric or logical values.
%
% Optional properties.
% Colormap : A N by 3 matrix of [r g b] triples, where N is the
%            length(adjacenyMatrix).
% Label    : A cell array of N strings.
%%
% Copyright 2014 The MathWorks, Inc.
  properties
    Node = node(0,0); % Array of nodes
    ColorMap;         % Colormap
    Label;            % Cell array of strings
    ShowButton;       % Turn all nodes on
    HideButton;       % Turn all nodes off
  end

  methods
    function this = circularGraph(adjacencyMatrix,varargin)
      % Constructor
      p = inputParser;

      defaultColorMap = parula(length(adjacencyMatrix));
      defaultLabel = cell(length(adjacencyMatrix));
      for i = 1:length(defaultLabel)
        defaultLabel{i} = num2str(i);
      end

      addRequired(p,'adjacencyMatrix',@(x)(isnumeric(x) || islogical(x)));
      addParameter(p,'ColorMap',defaultColorMap,@(colormap)length(colormap) == length(adjacencyMatrix));
      addParameter(p,'Label'   ,defaultLabel   ,@iscell);

      parse(p,adjacencyMatrix,varargin{:});
      this.ColorMap = p.Results.ColorMap;
      this.Label    = p.Results.Label;

      this.ShowButton = uicontrol(...
        'Style','pushbutton',...
        'Position',[0 40 80 40],...
        'String','Show All',...
        'Callback',@circularGraph.showNodes,...
        'UserData',this);

      this.HideButton = uicontrol(...
        'Style','pushbutton',...
        'Position',[0 0 80 40],...
        'String','Hide All',...
        'Callback',@circularGraph.hideNodes,...
        'UserData',this);

      fig = gcf;
      set(fig,...
        'UserData',this,...
        'CloseRequestFcn',@circularGraph.CloseRequestFcn);

      % Draw the nodes
      delete(this.Node);
      t = linspace(-pi,pi,length(adjacencyMatrix) + 1).'; % theta for each node
      extent = zeros(length(adjacencyMatrix),1);
      for i = 1:length(adjacencyMatrix)
        this.Node(i) = node(cos(t(i)),sin(t(i)));
        this.Node(i).Color = this.ColorMap(i,:);
        this.Node(i).Label = this.Label{i};
      end

      % Find non-zero values of s and their indices
      [row,col,v] = find(adjacencyMatrix);

      % Calculate line widths based on values of s (stored in v).
      minLineWidth  = 0.5;
      lineWidthCoef = 5;
      lineWidth = v./max(v);
      if sum(lineWidth) == numel(lineWidth) % all lines are the same width.
        lineWidth = repmat(minLineWidth,numel(lineWidth),1);
      else % lines of variable width.
        lineWidth = lineWidthCoef*lineWidth + minLineWidth;
      end

      % Draw connections on the Poincare hyperbolic disk.
      %
      % Equation of the circles on the disk:
      % x^2 + y^2
      % + 2*(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1))*x
      % - 2*(u(1)-v(1))/(u(1)*v(2)-u(2)*v(1))*y + 1 = 0,
      % where u and v are points on the boundary.
      %
      % Standard form of equation of a circle
      % (x - x0)^2 + (y - y0)^2 = r^2
      %
      % Therefore we can identify
      % x0 = -(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1));
      % y0 = (u(1)-v(1))/(u(1)*v(2)-u(2)*v(1));
      % r^2 = x0^2 + y0^2 - 1

      for i = 1:length(v)
        if row(i) ~= col(i)
          if abs(row(i) - col(i)) - length(adjacencyMatrix)/2 == 0
            % points are diametric, so draw a straight line
            u = [cos(t(row(i)));sin(t(row(i)))];
            v = [cos(t(col(i)));sin(t(col(i)))];
            this.Node(row(i)).Connection(end+1) = line(...
              [u(1);v(1)],...
              [u(2);v(2)],...
              'LineWidth', lineWidth(i),...
              'Color', this.ColorMap(row(i),:),...
              'PickableParts','none');
          else % points are not diametric, so draw an arc
            u  = [cos(t(row(i)));sin(t(row(i)))];
            v  = [cos(t(col(i)));sin(t(col(i)))];
            x0 = -(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1));
            y0 =  (u(1)-v(1))/(u(1)*v(2)-u(2)*v(1));
            r  = sqrt(x0^2 + y0^2 - 1);
            thetaLim(1) = atan2(u(2)-y0,u(1)-x0);
            thetaLim(2) = atan2(v(2)-y0,v(1)-x0);

            if u(1) >= 0 && v(1) >= 0
              % ensure the arc is within the unit disk
              theta = [linspace(max(thetaLim),pi,50),...
                       linspace(-pi,min(thetaLim),50)].';
            else
              theta = linspace(thetaLim(1),thetaLim(2)).';
            end

            this.Node(row(i)).Connection(end+1) = line(...
              r*cos(theta)+x0,...
              r*sin(theta)+y0,...
              'LineWidth', lineWidth(i),...
              'Color', this.ColorMap(row(i),:),...
              'PickableParts','none');
          end
        end
      end

      axis image;
      ax = gca;
      for i = 1:length(adjacencyMatrix)
        extent(i) = this.Node(i).Extent;
      end
      extent = max(extent(:));
      ax.XLim = ax.XLim + extent*[-1 1];
      fudgeFactor = 1.75; % Not sure why this is necessary. Eyeballed it.
      ax.YLim = ax.YLim + fudgeFactor*extent*[-1 1];
      ax.Visible = 'off';
      ax.SortMethod = 'depth';

      fig = gcf;
      fig.Color = [1 1 1];
    end

  end

  methods (Static = true)
    function showNodes(this,~)
      % Callback for 'Show All' button
      n = this.UserData.Node;
      for i = 1:length(n)
        n(i).Visible = true;
      end
    end

    function hideNodes(this,~)
      % Callback for 'Hide All' button
      n = this.UserData.Node;
      for i = 1:length(n)
        n(i).Visible = false;
      end
    end

    function CloseRequestFcn(this,~)
      % Callback for figure CloseRequestFcn
      c = this.UserData;
      for i = 1:length(c.Node)
        delete(c.Node(i));
      end
      delete(gcf);
    end

  end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [output] = volumefillholes(input, along)

% VOLUMEFILLHOLES is a helper function for segmentations
%
% See also VOLUMETHRESHOLD, VOLUMESMOOTH

% ensure that SPM is available, needed for spm_bwlabel
ft_hastoolbox('spm8up', 3) || ft_hastoolbox('spm2', 1);

if nargin<2
  inflate = false(size(input)+2);                   % grow the edges along each dimension
  inflate(2:end-1, 2:end-1, 2:end-1) = (input~=0);  % insert the original volume
  [lab, num] = spm_bwlabel(double(~inflate), 18);   % note that 18 is consistent with imfill, 26 is not
  if num>1
    inflate(lab~=lab(1)) = true;
    output  = inflate(2:end-1, 2:end-1, 2:end-1);   % trim the edges
  else
    output = input;
  end

else
  output = input;
  dim    = size(input);
  switch along
    case 1
      for i=1:dim(1)
        slice = reshape(input(i,:,:),dim([2 3]));
        im = imfill(slice,8,'holes');
        output(i,:,:) = im;
      end

    case 2
      for i=1:dim(2)
        slice = reshape(input(:,i,:),dim([1 3]));
        im = imfill(slice,8,'holes');
        output(:,i,:) = im;
      end

    case 3
      for i=1:dim(3)
        slice = reshape(input(:,:,i),dim([1 2]));
        im = imfill(slice,8,'holes');
        output(:,:,i) = im;
      end

    otherwise
      ft_error('invalid dimension %d to slice the volume', along);
  end % switch along
end % if nargin


text:
function [U, Theta, V, numiter ] = OR1MP(m, n, r, Known, data, opts )
%ORTHOGONAL-RANK-1-MATRIX-PURSUIT Infinit dimension matching pursuit for
%low rank matrix
%
%For the problem           min    L(X) = ||X - Y||^2
%                          s.t.   rank(X) <= r
%                          X = Theta * M = sum_i theta_i * M_i = U Theta V'
%   Detail variable
%   Input:
%         m ---- row number
%         n ---- column number
%         r ---- number of basis
%         Known ---- index of the known spot in the matrix
%         data ---- the content of the known spot in the matrix
%         opts ---- parameters for the algorithm
%   Output:
%         U ---- output matrix: U
%         Theta ---- output matrix: Theta
%         V ---- output matrix: V
%         numiter ---- number of iterations
%
%   Copyright Zheng Wang @ Arizona State University
%   $Date: 2013/01/29$

fn = mfilename;
error(nargchk(5, 6, nargin));
if nargin < 6
    epsilon = 1e-4; % convergence threshold
else
    epsilon = opts.epsilon;
end

%addpath('PROPACK/');
%addpath('largescale_ops/');
%addpath('SLEP_package_4.1/');

% initialization,
[indm, indn] = ind2sub([m, n], Known);
data( data == 0 )= eps;
res = sparse(indm, indn, data, m, n);
[indm, indn, data] = find(res);
U = [];
V = [];
Msup = [];

verbosity = 1;
printsyb = ['-', 'X', '|'];
if verbosity == 1
    fprintf('\nIteration:        ');
end

i = 0;
W = 0;
oldresnorm = 0;
gresnorm = 1;
yy = [];
nnorm = norm(data, 'fro');
% main iteration
% In OR1MP, the stop criterion is small gradient of residual
while (i < r) && (gresnorm > epsilon )
    % 1. find the top singular pair of the residual and update the gresnorm
    resvec = data - W;
    sparse_update(res, resvec); % sparse update the res using resvec

    [u, ~, v] = topsvd(res, 20); % run our power method for 10 iterations
    %[u, ~, v] = topsvd(res, 1);
    %[u, ~, v] = lansvd(res, 1, 'L'); % fast sparse svd using PROPACK
    %[u, s, v] = svds(res, 1); % use matlab sparse top svd

    resnorm = normest(res, 'fro')/nnorm;
    gresnorm = abs(resnorm - oldresnorm);
    oldresnorm = resnorm;

    % 2. update the weight Theta, the pursuit basis is uv', its weight is s.
    Mi = sparse_inp(u', v', indm, indn)';

    % b) use incremental inverse to solve the least sqare problem
    if i~=0
        Minv = inverse_incremental(Minv, Msup'*Mi, Mi'*Mi);
    else
        Minv = 1/(Mi'*Mi);
    end
    yy = [yy; Mi'*data];
    Theta = Minv*yy;
    Msup = [Msup Mi];

    U = [U u];
    V = [V v];

    % 3. update the learned matrix W = U' * diag(Theta) * V;
    W = Msup * Theta;

    if verbosity == 1
        fprintf('\b\b\b\b\b\b\b\b  %c:%4d', printsyb(1+mod(i,length(printsyb))), i);
    end
    i = i + 1;
end
% V = diag(Theta)*V;
numiter = i;
% fprintf( '\n OR1MP run %d rounds! \n', numiter);

% %* Sparse selection: for sharpe low rank problem, we may need the lasso
% %fine selection. Use lasso to learn a more sparse weights.
% % Starting point
% opts.init=2;        % starting from a zero point
% % termination criterion
% opts.tFlag=5;       % run .maxIter iterations
% opts.maxIter=100;   % maximum number of iterations
% % normalization
% opts.nFlag=0;       % without normalization
% % regularization
% opts.rFlag=1;       % the input parameter 'rho' is a ratio in (0, 1)
% opts.mFlag=0;       % treating it as compositive function
% opts.lFlag=0;       % Nemirovski's line search
% % get the final sparse Theta by lasso in SLEP toolbox
% [Theta, ~, ~] = LeastR(Msup, data, 0.00001, opts); %  [W, funVal, ValueL] 00001

function Ninv = inverse_incremental(Minv, MMi, d)
% calculate the inverse of the blocked matrix of
P = MMi' * Minv; % vector
q = 1/(d - P*MMi); % scaler
y = q*P;
Ninv = [Minv+P'*y, -y'; -y, q];

function [u, s, v] = topsvd(A, round)
% calculate the top svd of matrix A using round iterations
stopeps = 1e-3;
[m,n] = size(A);
u = ones(m,1); % this is the sigma*u
vo = 0;
for i=1:round
    v = u'*A/(norm(u))^2;
    u = A*v'/(norm(v))^2;
    if norm(v-vo) < stopeps
        break
    end
    vo = v;
end
nu = norm(u);
nv = norm(v);
u = u/nu;
v = v'/nv;
% v = v/nv;
s = nu*nv;

author: andrewssobral

repo: lrslibrary


text:
classdef spinpreference
%SPINPREFERENCE   Abstract class for managing preferences for SPIN, SPIN2 and
%SPIN3.
%   SPINPREFERENCE is an abstract class for managing preferences when solving
%   a time-dependent PDE defined by a SPINOPERATOR. SPINPREF (for SPIN in 1D),
%   SPINPREF2 (for SPIN2 in 2D), SPINPREF3 (for SPIN3 in 3D) and SPINPREFSPHERE
%   (for SPINSPHERE on the sphere) are full implementations.
%
% See also SPINPREF, SPINPREF2, SPINPREF3, SPINPREFSPHERE.

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )
        dataplot             % Which data to plot when complex values (STRING)
        dealias              % To use dealiasing with 2/3-rule (STRING)
        iterplot             % Plot every ITERPLOT iterations (1x1 INT)
        Nplot                % Number of points for movie grid (1x1 INT)
        plot                 % Plot options (STRING)
        scheme               % Time-stepping scheme (STRING)
    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef GemButton < GemUserInterfaceObject
    % GemButton GEM class for a button control
    %
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        % The button width and height may need to be set if an image is to be applied to the button
        % before it is resized
        ButtonWidth
        ButtonHeight

        FontSize = 11
        BorderSize = 1
        FontAngle = 'normal'
        Text
        ToolTip
        RGBImage
        OriginalImage
        ImageHighlight
        BackgroundColour = [0, 0, 0]

        Selected = false
        Highlighted = false;

        ShowTextOnButton

        SelectedColour
        UnSelectedColour
        HighlightColour
        HighlightSelectedColour

        AutoUpdateStatus = false
    end

    events
        ButtonClicked
    end

    properties (Access = protected)
        Tag
        Callback
    end

    methods
        function obj = GemButton(parent, text, tooltip, tag, callback)
            obj = obj@GemUserInterfaceObject(parent);
            obj.Text = text;
            obj.ToolTip = tooltip;
            obj.Tag = tag;
            obj.Callback = callback;
            obj.ShowTextOnButton = true;
            obj.SelectedColour = uint8(255*(obj.StyleSheet.SelectedBackgroundColour));
            obj.UnSelectedColour = uint8([150, 150, 150]);
            obj.HighlightColour = min(255, obj.UnSelectedColour + 100);
            obj.HighlightSelectedColour = min(255, obj.SelectedColour + 100);
        end

        function AddAndResizeImageWithHighlightMask(obj, new_image, mask_colour_rgb)

            if isempty(new_image)
                obj.OriginalImage = [];
            else
                new_image_size = [obj.ButtonHeight - 2*obj.BorderSize, obj.ButtonWidth - 2*obj.BorderSize];
                rgb_image = zeros([obj.ButtonHeight, obj.ButtonWidth, 3], class(new_image));
                rgb_image(1+obj.BorderSize:end-obj.BorderSize, 1+obj.BorderSize:end-obj.BorderSize, 1) = imresize(new_image(:, :, 1), new_image_size, 'bilinear');
                rgb_image(1+obj.BorderSize:end-obj.BorderSize, 1+obj.BorderSize:end-obj.BorderSize, 2) = imresize(new_image(:, :, 2), new_image_size, 'bilinear');
                rgb_image(1+obj.BorderSize:end-obj.BorderSize, 1+obj.BorderSize:end-obj.BorderSize, 3) = imresize(new_image(:, :, 3), new_image_size, 'bilinear');
                background = rgb_image(:, :, 1) == mask_colour_rgb(1) & rgb_image(:, :, 2) == mask_colour_rgb(2) & rgb_image(:, :, 3) == mask_colour_rgb(3);

                obj.ImageHighlight = GemUtilities.GetRGBImageHighlight(rgb_image, mask_colour_rgb);

                for index = 1 : 3
                    rgb_slice = rgb_image(:, :, index);
                    rgb_slice(background) = round(255*obj.StyleSheet.BackgroundColour(index));
                    rgb_image(:, :, index) = rgb_slice;
                end
                obj.OriginalImage = rgb_image;
            end

            obj.UpdateImageHighlight;
        end

        function CreateGuiComponent(obj, position)

            % If no image has been specified then create a blank image
            if isempty(obj.OriginalImage)
                obj.OriginalImage = GemUtilities.GetBlankButtonImage(position(3), position(4), obj.BorderSize, obj.BackgroundColour, obj.UnSelectedColour);
                button_size = size(obj.RGBImage);
                obj.ImageHighlight = GemUtilities.GetBorderImage(button_size(1:2), obj.BorderSize);
            end

            obj.UpdateRGBImageCache;

            if obj.ShowTextOnButton
                button_text = obj.Text;
            else
                button_text = '';
            end

            obj.GraphicalComponentHandle = uicontrol('Style', 'pushbutton', 'Parent', obj.Parent.GetContainerHandle, ...
                'String', button_text, 'Tag', obj.Tag, 'ToolTipString', obj.ToolTip, ...
                'FontAngle', obj.FontAngle, 'ForegroundColor', obj.StyleSheet.TextPrimaryColour, 'FontUnits', 'pixels', 'FontSize', obj.FontSize, ...
                'Callback', @obj.ButtonClickedCallback, 'Position', position, 'CData', obj.RGBImage);
        end

        function ChangeImage(obj, preview_image)
            obj.OriginalImage = preview_image;
            button_size = size(obj.RGBImage);
            obj.ImageHighlight = GemUtilities.GetBorderImage(button_size(1:2), obj.BorderSize);

            obj.UpdateImageHighlight
        end

        function Select(obj, is_selected)
            if obj.Selected ~= is_selected
                obj.Selected = is_selected;
                obj.UpdateImageHighlight;
            end
        end
    end

    methods (Access = protected)
        function ButtonClickedCallback(obj, ~, ~)

            if ~obj.IsCurrentlyRunning()
                obj.Select(true);
                notify(obj, 'ButtonClicked', CoreEventData(obj.Tag));
                obj.Callback(obj.Tag);

                % It is possible that the callback may lead to a rebuild of the interface and
                % thus deletion of the button that triggered the callback; in which case we
                % can't continue with the select
                if isvalid(obj) && ~obj.AutoUpdateStatus
                    obj.Select(false);
                end
            end
        end

        function input_has_been_processed = MouseHasMoved(obj, click_point, selection_type, src, eventdata)
            % This method is called when the mouse is moved

            obj.Highlight(true);
            input_has_been_processed = true;
        end

        function input_has_been_processed = MouseExit(obj, click_point, selection_type, src, eventdata)
            % This method is called when the mouse exits a control which previously
            % processed a MouseHasMoved event

            obj.Highlight(false);
            input_has_been_processed = true;
        end

    end

    methods (Access = private)
        function Highlight(obj, highlighted)
            if (highlighted ~= obj.Highlighted)
                obj.Highlighted = highlighted;
                obj.UpdateImageHighlight;
            end
        end

        function UpdateImageHighlight(obj)
            obj.UpdateRGBImageCache;

            if obj.ComponentHasBeenCreated && ishandle(obj.GraphicalComponentHandle)
                set(obj.GraphicalComponentHandle, 'CData', obj.RGBImage);
            end
        end

        function UpdateRGBImageCache(obj)
            % Take the button image from the OriginalImage cache property, and a highlight
            % of the appropriate colour, and then put into the RGBImage cache property

            if isempty(obj.OriginalImage)
                obj.RGBImage = [];
            else
                if obj.Selected
                    if obj.Highlighted
                        highlight_colour = obj.HighlightSelectedColour;
                    else
                        highlight_colour = obj.SelectedColour;
                    end
                else
                    if obj.Highlighted
                        highlight_colour = obj.HighlightColour;
                    else
                        highlight_colour = obj.UnSelectedColour;
                    end
                end
                obj.RGBImage = GemUtilities.AddHighlightToRGBImage(obj.OriginalImage, obj.ImageHighlight, highlight_colour);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function trans = rotationOz(varargin)
%ROTATIONOZ return 4x4 matrix of a rotation around z-axis
%
%   TRANS = rotationOz(THETA);
%   Returns the transform matrix corresponding to a rotation by the angle
%   THETA (in radians) around the Oz axis. A rotation by an angle of PI/2
%   would transform the vector [1 0 0] into the vector [0 1 0].
%
%   The returned matrix has the form:
%   [cos(THETA) -sin(THETA)  0  0]
%   [sin(THETA)  cos(THETA)  0  0]
%   [    0           0       1  0]
%   [    0           0       0  1]
%
%   TRANS = rotationOz(ORIGIN, THETA);
%   TRANS = rotationOz(X0, Y0, Z0, THETA);
%   Also specifies origin of rotation. The result is similar as performing
%   translation(-X0, -Y0, -Z0), rotation, and translation(X0, Y0, Z0).
%
%
%   See also:
%   transforms3d, transformPoint3d, rotationOx, rotationOy
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/04/2004.
%

%   HISTORY
%   2008/11/24 changed convention for angle
%   30/04/2009 deprecate: use createRotationOz instead

% deprecation warning
warning('geom3d:deprecated', ...
    [mfilename ' is deprecated, use ''createRotationOz'' instead']);

% call current implementation
trans = createRotationOz(varargin{:});

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function download_algorithms(simulate_no_download)
% Downloads and compiles/preps different external codes needed for RIGOR
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

clc

% if you don't want to download the files (for testing this script)
if ~exist('simulate_no_download', 'var')
    simulate_no_download = false;
end
no_download_dir = 'pre_downloaded';

fprintf(2, '---------------------------------------------------------------------------\n                              Download Notice                              \n---------------------------------------------------------------------------\nBy running this script, you accept all licensing agreements accompanied\nwith the 3rd party softwares that will now be downloaded and used later in\nour scripts.\n\nPress ''y'' to accept (any other key to stop the script): ');
user_inp = input('','s');
if isempty(user_inp) || ~strncmpi(user_inp(1), 'y', 1)
    return;
end

% destination algorithms directory
curr_dir = pwd;
code_root_dir = fullfile(fileparts(which(mfilename)), '..');
extern_src_rel = 'extern_src';
algos_dir = fullfile(code_root_dir, extern_src_rel);
dont_delete = {'DataHash', 'extra_gb_code', 'para_pseudoflow', 'fuxin_lib_src', 'stein_boundaryprocessing'};
if exist(algos_dir,'dir')
    % only delete the dir/files not in dont_delete
    d = dir(algos_dir);
    for idx = 1:length(d)
        if strcmp(d(idx).name,'.') || strcmp(d(idx).name,'..')
            continue;
        end
        if ~any(strcmp(d(idx).name, dont_delete))
            if d(idx).isdir
                rmdir(fullfile(algos_dir,d(idx).name), 's');
            else
                delete(fullfile(algos_dir,d(idx).name));
            end
        end
    end
else
    mkdir(algos_dir);
end

% create temp download directory
if ~simulate_no_download
    temp_dir = fullfile(code_root_dir, sprintf('temp%d', randi(1e8)));
    mkdir(temp_dir);
else
    temp_dir = no_download_dir;
end

try
    % Point Fuxin library to the right directory
    fuxin_lib_rel = fullfile(extern_src_rel, 'fuxin_lib_src');
    replaceInTextFile(fullfile(code_root_dir, 'internal_params.m'), '''/home/ahumayun/videovolumes/fuxin_lib_src''', ['fullfile(fp.code_root_dir, ''', fuxin_lib_rel, ''')'], true);
    classregtree_dir = fullfile(fuxin_lib_rel, '@classregtree_fuxin', 'private');
    fprintf(1, 'mex''ing fuxin library code\n');
    eval(sprintf('mex -O %s/regtreeEval.cpp -output %s/regtreeEval', classregtree_dir, classregtree_dir));
    eval(sprintf('mex -O %s/regtree_findbestsplit.cpp -output %s/regtree_findbestsplit', classregtree_dir, classregtree_dir));
    fprintf(1, 'Done mex''ing\n');

    % move stein boundary processing code to the segmentation folder
    mkdir(fullfile(algos_dir, 'segmentation'));
    copyfile(fullfile(algos_dir, 'stein_boundaryprocessing'), fullfile(algos_dir, 'segmentation', 'stein_boundaryprocessing'));

    % download Leordeanu''s GB code
    success = download_code('http://109.101.234.42/documente/code/doc_8.zip', ...
        fullfile(temp_dir, 'doc_8.zip'), 'Leordeanu''s GB flow code', simulate_no_download);

    if success
        mkdir(fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb'));
        unzip(fullfile(temp_dir, 'doc_8.zip'), fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb'));
        movefile(fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb', 'Gb_Code_Oct2012', '*'), fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb'));
        rmdir(fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb', 'Gb_Code_Oct2012'), 's');
        replaceInTextFile(fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb', 'Gb_data_lambda.m'), 'gb\(f\) = T/2 \+ sqrt\(\(T \.\^ 2 \) /4 \- D\);', 'gb(f) = T(f)/2 + sqrt((T(f) .^ 2 ) /4 - D(f));');
        replaceInTextFile(fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb', 'Gb_data_lambda.m'), 'or_C\(f, 1\) = -Ms\(:, 2\);', 'or_C(f, 1) = -Ms(f, 2);');
        replaceInTextFile(fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb', 'Gb_data_lambda.m'), 'or_C\(f, 2\) = gb\(:\) - Ms\(:, 1\);', 'or_C(f, 2) = gb(f) - Ms(f, 1);');
        adjustAttributes(fullfile(algos_dir, 'segmentation', 'boundaries--leordeanu_ECCV_2012_gb'));
    end

    % download Joseph Lim's SketchTokens code
    success = download_code('https://github.com/joelimlimit/SketchTokens/archive/master.zip', ...
        fullfile(temp_dir, 'master.zip'), 'Joseph Lim''s SketchTokens code', simulate_no_download);
    success = success & download_code('http://people.csail.mit.edu/lim/lzd_cvpr2013/st_data.tgz', ...
        fullfile(temp_dir, 'st_data.tgz'), 'Joseph Lim''s SketchTokens data', simulate_no_download);

    if success
        mkdir(fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens'));
        unzip(fullfile(temp_dir, 'master.zip'), fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens'));
        untar(fullfile(temp_dir, 'st_data.tgz'), fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens'));
        rmdir(fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens', 'SketchTokens-master', 'models'), 's');
        movefile(fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens', 'models'), fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens', 'SketchTokens-master'));
        movefile(fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens', 'SketchTokens-master', '*'), fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens'));
        rmdir(fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens', 'SketchTokens-master'), 's');
        cd(fullfile(algos_dir, 'segmentation', 'boundaries--lim_CVPR_2013_sketchtokens'));
        mex stDetectMex.cpp 'OPTIMFLAGS="$OPTIMFLAGS' '/openmp"'
        fprintf(1, 'Done mex''ing\n');
        cd(curr_dir);
    end

    % download Piotr's toolbox code
    success = download_code('http://vision.ucsd.edu/~pdollar/toolbox/piotr_toolbox.zip', ...
        fullfile(temp_dir, 'piotr_toolbox.zip'), 'Piotr''s toolbox code', simulate_no_download);

    if success
        mkdir(fullfile(algos_dir, 'toolboxes', 'piotr_toolbox'));
        unzip(fullfile(temp_dir, 'piotr_toolbox.zip'), fullfile(algos_dir, 'toolboxes', 'piotr_toolbox'));
        movefile(fullfile(algos_dir, 'toolboxes', 'piotr_toolbox', 'toolbox', '*'), fullfile(algos_dir, 'toolboxes', 'piotr_toolbox'));
        rmdir(fullfile(algos_dir, 'toolboxes', 'piotr_toolbox', 'toolbox'), 's');
        % do not compile for linux - there is some problem when running
        % struct edges with compiled code (rather than with the compiled
        % code that comes with the library)
        if ~isunix
            cd(fullfile(algos_dir, 'toolboxes', 'piotr_toolbox', 'external'));
            toolboxCompile;
            fprintf(1, 'Done mex''ing\n');
            cd(curr_dir);
        end
    end

    % download Piotr's Struct Edges code
    success = download_code('http://ftp.research.microsoft.com/downloads/389109f6-b4e8-404c-84bf-239f7cbf4e3d/releaseV3.zip', ...
        fullfile(temp_dir, 'release.zip'), 'Piotr''s Struct Edges code', simulate_no_download);

    if success
        mkdir(fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges'));
        unzip(fullfile(temp_dir, 'release.zip'), fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges'));
        movefile(fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges', 'releaseV3', '*'), fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges'));
        rmdir(fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges', 'releaseV3'), 's');
        cd(fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges'));
        if isunix
            mex private/edgesDetectMex.cpp -outdir private '-DUSEOMP' CFLAGS="\$CFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
            mex private/edgesNmsMex.cpp    -outdir private '-DUSEOMP' CFLAGS="\$CFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
        else
            % if windows
            mex private/edgesDetectMex.cpp -outdir private '-DUSEOMP' 'OPTIMFLAGS="$OPTIMFLAGS' '/openmp"'
            mex private/edgesNmsMex.cpp    -outdir private '-DUSEOMP' 'OPTIMFLAGS="$OPTIMFLAGS' '/openmp"'
        end
        fprintf(1, 'Done mex''ing\n');
        cd(curr_dir);
        replaceInTextFile(fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges', 'edgesDetect.m'), '\[E,O,inds,segs\]', '[E,T,O,inds,segs]');
        replaceInTextFile(fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges', 'edgesDetect.m'), '% perform nms', 'T = [];\n% perform nms');
        replaceInTextFile(fullfile(algos_dir, 'segmentation', 'boundaries--dollar_ICCV_2013_structedges', 'edgesDetect.m'), 'E\s*=\s*edgesNmsMex', 'T=edgesNmsMex');
    end

    % download Woodford's Image-Based Rendering and Stereo code
    success = download_code('http://www.robots.ox.ac.uk/~ojw/files/imrender_v2.4.zip', ...
        fullfile(temp_dir, 'imrender_v2.4.zip'), 'Woodford''s Image-Based Rendering and Stereo code', simulate_no_download);

    if success
        unzip(fullfile(temp_dir, 'imrender_v2.4.zip'), fullfile(algos_dir, 'segmentation'));
        cd(fullfile(algos_dir, 'segmentation', 'imrender', 'vgg'));
        % segment as a test - which also compiles the file
        temp = vgg_segment_gb(imread('peppers.png'), 0.5, 10, 10);
        fprintf(1, 'Done mex''ing\n');
        cd(curr_dir);
    end

    % download export_fig code
    success = download_code('https://codeload.github.com/ojwoodford/export_fig/legacy.zip/master', ...
        fullfile(temp_dir, 'ojwoodford-export_fig.zip'), 'Woodford''s matlab export fig code', simulate_no_download);

    if success
        mkdir(fullfile(algos_dir, 'utils'));
        unzip(fullfile(temp_dir, 'ojwoodford-export_fig.zip'), fullfile(algos_dir, 'utils'));
        d = dir(fullfile(algos_dir, 'utils', '*export_fig*'));
        movefile(fullfile(algos_dir, 'utils', d(1).name), fullfile(algos_dir, 'utils', 'export_fig'));
        replaceInTextFile(fullfile(code_root_dir, 'utils', 'drawFigFrames.m'), '~/videovolumes/extern_src', algos_dir, true);
    end

    % download Vedaldi's vlfeat code
    success = download_code('https://github.com/vlfeat/vlfeat/archive/master.zip', ...
        fullfile(temp_dir, 'vlfeat.zip'), 'Vedaldi''s vlfeat code', simulate_no_download);

    if success
        unzip(fullfile(temp_dir, 'vlfeat.zip'), fullfile(algos_dir, 'toolboxes'));
        d = dir(fullfile(algos_dir, 'toolboxes', '*vlfeat*'));
        movefile(fullfile(algos_dir, 'toolboxes', d(1).name), fullfile(algos_dir, 'toolboxes', 'vlfeat'));
    end

    % change the extern code path in the internal parameters script
    replaceInTextFile(fullfile(code_root_dir, 'internal_params.m'), '''/home/ahumayun/videovolumes/extern_src''', ['fullfile(fp.code_root_dir, ''', extern_src_rel, ''')'], true);

    % delete temp directory
    if ~simulate_no_download
        rmdir(temp_dir, 's');
    end

catch exception
    % remove temp dir
    rmdir(temp_dir, 's');
    rethrow(exception)
end

end

function success = download_code(url, filepath, desc, simulate_no_download)
try
    fprintf(1, 'Downloading %s ...\n', desc);
    if ~simulate_no_download
        [f, status] = urlwrite(url, filepath);
        success = status == 1;
    else
        success = true;
    end
catch exception
    success = false;
end

if success
    if ~simulate_no_download
        fprintf(1, 'Done downloading\n');
    else
        fprintf(1, 'Simulated download call\n');
    end
else
    fprintf(2, 'Downloading failure\n');
end
end

function transferMatlabFromGNU(makefilepath, dest_dir)
fd = fopen(makefilepath, 'r');
str = fread(fd);
str = char(str');
fclose(fd);

files = {};
re_matches = regexp(str, '(?:matlab)\s+:=\s+(\S+.(?:(mat)|m)(?:\s*\\?\s*((\$\(wildcard)|(\)))?\s+))*', 'tokens');
for idx = 1:length(re_matches)
    f = regexp(re_matches{idx}{1}, '(\S+.(?:(mat)|m))', 'tokens');
    files = [files cellfun(@(x)x, f)];
end
for idx = 1:length(files)
    movefile(fullfile(fileparts(makefilepath), files{idx}), dest_dir);
end
end

function adjustAttributes(folder_path)
% dont need to change file permissions on windows
if ispc == 1
    return;
end

d = dir(folder_path);
for idx = 1:length(d)
    if strcmp(d(idx).name,'.') || strcmp(d(idx).name,'..')
        continue;
    end

    curr_path = fullfile(folder_path,d(idx).name);
    if d(idx).isdir == 1
        adjustAttributes(curr_path);
    else
        unix(['chmod 0644 "' curr_path '"']);
    end
end
end

author: Cloud-CV

repo: object-proposals


text:
function new_plugin = MimParseGuiPluginClass(plugin_name, plugin_class, suggested_category, default_category, default_mode, reporting)
    % MimParseGuiPluginClass. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     Fetches information about any GUI plugins for the TD MIM Toolkit which
    %     are available.
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    new_plugin = [];
    new_plugin.PluginName = plugin_name;

    new_plugin.ToolTip = plugin_class.ToolTip;
    new_plugin.ButtonText = plugin_class.ButtonText;
    if ~isempty(plugin_class.Category)
        new_plugin.Category = plugin_class.Category;
    else
        if ~isempty(suggested_category)
            new_plugin.Category = suggested_category;
        else
            new_plugin.Category = default_category;
        end
    end
    new_plugin.HidePluginInDisplay = plugin_class.HidePluginInDisplay;
    new_plugin.ButtonWidth = plugin_class.ButtonWidth;
    new_plugin.ButtonHeight = plugin_class.ButtonHeight;
    new_plugin.AlwaysRunPlugin = false; % Ensures plugin name is not in italics

    property_list = properties(plugin_class);

    if ismember('PluginInterfaceVersion', property_list)
        new_plugin.PluginInterfaceVersion = plugin_class.PluginInterfaceVersion;
    elseif ismember('PTKVersion', property_list)
        new_plugin.PluginInterfaceVersion = plugin_class.PTKVersion;
    else
        new_plugin.PluginInterfaceVersion = 1;
    end

    if ismember('Visibility', property_list)
        new_plugin.Visibility = plugin_class.Visibility;
    end

    if ismember('Mode', property_list) && ~isempty(plugin_class.Mode)
        new_plugin.Mode = plugin_class.Mode;
    else
        new_plugin.Mode = default_mode;
    end

    if ismember('SubMode', property_list)
        new_plugin.SubMode = plugin_class.SubMode;
    else
        new_plugin.SubMode = [];
    end

    if ismember('Location', property_list)
        new_plugin.Location = plugin_class.Location;
    else
        new_plugin.Location = 100;
    end
end


text:
 function [hash_table, hash_coarse] = at_dense_hashtable(cnnfeat1,cnnfeat1fine)

x_coarse_size = size(cnnfeat1,1);
y_coarse_size = size(cnnfeat1,2);

x_fine_size = size(cnnfeat1fine,1);
y_fine_size = size(cnnfeat1fine,2);

% scale = x_fine_size/x_coarse_size;
% if scale ~= whos
%   error('aspect ratio should be preserved');
% end

% x_coarse_size = 5;
% y_coarse_size = 4;
% scale = 2;
% x_fine_size = x_coarse_size * scale;
% y_fine_size = y_coarse_size * scale;

% [x_coarse,y_coarse] = meshgrid(1:x_coarse_size,1:y_coarse_size);
hash_coarse = reshape(1:(x_coarse_size*y_coarse_size),x_coarse_size,y_coarse_size);

hash_fine = imresize(hash_coarse,[x_fine_size y_fine_size],'nearest');
[x_fine,y_fine] = meshgrid(1:y_fine_size,1:x_fine_size);

Nhash = max(hash_coarse(:));

hash_table = cell(1,Nhash);
hash_fine = hash_fine(:);
x_fine = x_fine(:);
y_fine = y_fine(:);
for ii=1:Nhash
  a = find(hash_fine == ii);
  hash_table{ii} = [x_fine(a)'; y_fine(a)'];
end

author: HajimeTaira


text:
function [week, sow] = time2weektow(time)

% SYNTAX:
%   [week, sow] = time2weektow(time);
%
% INPUT:
%   time = GPS time (continuous since 6-1-1980)
%
% OUTPUT:
%   week = GPS week
%   sow  = GPS seconds-of-week
%
% DESCRIPTION:
%   Conversion from GPS time in continuous format (similar to datenum) to
%   GPS time in week, seconds-of-week.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

sec_in_week = 7*86400;

sow  = rem(time, sec_in_week);
week = (time - sow) / sec_in_week;

author: goGPS-Project


text:
% Copyright 2017 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     https://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Returns set of parameters
% set light_mode = true to run the code in a sub optimal but faster mode
% set light_mode = false to obtain the results reported in the RED paper

function params = GetSuperResADMMParams(light_mode)

params = GetSuperResFPParams(light_mode);

% admm parameter
params.beta = 1e-3;

% number of denoising applications
params.inner_denoiser_iters = 1;

% relaxation parameter of ADMM
params.alpha = 1;

return

author: google

repo: RED

save_path: github-repos/MATLAB/google-RED


text:
%kkbitset 'Set All Bits to One'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kbitset.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkbitset(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kbitset - Set All Bits to One
%
%  DESCRIPTION
% The "SET" operator sets all bits to one for each data
% point in the input data object, \fBInput".
%
% Executing "SET" runs the program \fIkarith1\fP with the -set flag.
%
%
%  "\fBData Type"" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "\fBMap Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "\fBValidity Mask"" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "\fBExplicit Location and Time Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "\fBFailure Modes""
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith1
%
%  RESTRICTIONS
% Operations on complex data are not supported at this time.
% For the bitwise operations (not, set, clear), if the source data is
% floating point, it is cast to long integer before processing.  This
% may result in loss of precision.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkbitset(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkbitset(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith1"  -set'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%scf - Set current figure, without making it visible.
%
%  This function is useful in batches where figures must remain invisible. It
%  is equivalent to figure(h) but does not make h visible.
%
%  USAGE
%
%    scf(h)
%
%    h              figure handle
%
%  SEE ALSO
%
%    See also Hide, sca, StartBatch.
%

% Copyright (C) 2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function scf(h)

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

% Check parameters
if ~ishandle(h) || ~strcmp(get(h,'type'),'figure'),
  error('Incorrect figure handle (type ''help scf'' for details).');
end

set(0,'CurrentFigure',h);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [X, z, mu] = mixBernRnd(d, k, n)
% Generate samples from a Bernoulli mixture distribution.
% Input:
%   d: dimension of data
%   k: number of components
%   n: number of data
% Output:
%   X: d x n data matrix
%   z: 1 x n response variable
%   mu: d x k parameters of each Bernoulli component
% Written by Mo Chen (sth4nth@gmail.com).

% w = dirichletRnd(1,ones(1,k)/k);
w = ones(1,k)/k;
z = discreteRnd(w,n);
mu = rand(d,k);
X = zeros(d,n);
for i = 1:k
    idx = z==i;
    X(:,idx) = bsxfun(@le,rand(d,sum(idx)), mu(:,i));
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%  OVERLAPPING BLOCK-BASED SIGNAL SUBSPACE PROCESSING  %%
   %%                  FOR MOTION TRACKING                 %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nx=2;                     % filter size in x domain is 2*nx+1
ny=2;                     % filter size in y domain is 2*ny+1
N=(2*nx+1)*(2*ny+1);      % Number of basis functions
%
% ref_image:              reference image or f1(x,y) (INPUT)
% test_image:             test image or f2(x,y)      (INPUT)

% Example:
%
      NX=48; NY=44;
      square=randn(10:10); rect=randn(20,14);
      ix=-10:10; ix=ix(:)*ones(1,17);
      jy=-8:8; jy=ones(21,1)*jy;
      ell=(sqrt((ix/10).^2+(jy/8).^2) <= 1).*randn(21,17);
      ref_image=zeros(NX,NY);
      ref_image(5:14,17:26)=square;
      ref_image(19:38,5:18)=rect;
      ref_image(15:35,20:36)=ell;
      ref_image=ref_image+.05*randn(NX,NY);
      test_image=zeros(NX,NY);
      test_image(5+1:14+1,17-1:26-1)=2*square;
      test_image(19:38,5+2:18+2)=-1.4*rect;
      test_image(15-1:35-1,20-2:36-2)=.9*ell;
      test_image=test_image+.05*randn(NX,NY);

%
[NX,NY]=size(ref_image);  % size of reference or test image
%
nxb=5;                    % size of block in x domain
nxb2=(nxb-1)/2;
nyb=5;                   % size of block in y domain
nyb2=(nyb-1)/2;
nb=nxb*nyb;               % number of pixels in a block
%
MX=NX+2*(nx+nxb2);        % size of zero-padded image in x domain
MY=NY+2*(ny+nyb2);        % size of zero-padded image in y domain
%
f1=zeros(MX,MY);
f1(nx+nxb2+1:nx+nxb2+NX,ny+nyb2+1:ny+nyb2+NY)=ref_image;
                                     % zero-padded reference image
%
f2=zeros(MX,MY);
f2(nx+nxb2+1:nx+nxb2+NX,ny+nyb2+1:ny+nyb2+NY)=test_image;
                                     % zero-padded test image
%
% for each block perform signal subspace processing
%
hx=zeros(MX,MY);
hy=hx;
ix=-nx:nx; jy=-ny:ny;

for i=nx+nxb2+1:MX-nx-nxb2; i
 IX=i-nxb2:i+nxb2;         % x domain block indices
 for j=ny+nyb2+1:MY-ny-nyb2;
  JY=j-nyb2:j+nyb2;        % y domain block indices
%
  g=zeros(nb,N);           % array containing block reference image
                           % and its shifted versions
  icount=0;
  for ii=-nx:nx;
   for jj=-ny:ny;
    icount=icount+1;
    g(:,icount)=reshape(f1(IX+ii,JY+jj),nb,1);
   end;
  end;
  [U,D,V]=svd(g);               % perform svd
  psi=U(:,1:N).';               % orthogonal basis functions
  g2=reshape(f2(IX,JY),nb,1);   % test image block
  F2=conj(psi)*g2;              % projection coefficients of
                                % test image
  F1=conj(psi)*g;          % projection coefficients of
                           % reference image and its shifted versions
  H=pinv(F1)*F2;                % estimate of filter
  h=zeros(2*nx+1,2*ny+1);
  icount=0;
  for ii=1:2*nx+1;
     for jj=1:2*ny+1;
        icount=icount+1;
        h(ii,jj)=H(icount);
     end;
  end;

  % estimate shifts
  %
  Eh=sum(sum(h.^2));
  hx(i,j)=sum(sum((ix(:)*ones(1,2*ny+1)).*(h.^2)))/Eh;  % x domain shift
  hy(i,j)=sum(sum((ones(2*nx+1,1)*jy).*(h.^2)))/Eh;     % y domain shift
 end;
end;
hx=hx.*(abs(hx) <= nx);
hy=hy.*(abs(hy) <= ny);

% lowpass filter estimate of motion vectors
%
ix=-MX/2:MX/2-1; jy=-MY/2:MY/2-1;
winx=(.54+.46*cos(pi*(ix/(.4*MX)))).*(abs(ix) <= .4*MX);
winy=(.54+.46*cos(pi*(jy/(.4*MY)))).*(abs(jy) <= .4*MY);
W=winx(:)*winy;
HX=real(iftx(ifty(ftx(fty(hx)).*W)));
HY=real(iftx(ifty(ftx(fty(hy)).*W)));

hx=hx(nx+nxb2+1:nx+nxb2+NX,ny+nyb2+1:ny+nyb2+NY);
hy=hy(nx+nxb2+1:nx+nxb2+NX,ny+nyb2+1:ny+nyb2+NY);
HX=HX(nx+nxb2+1:nx+nxb2+NX,ny+nyb2+1:ny+nyb2+NY);
HY=HY(nx+nxb2+1:nx+nxb2+NX,ny+nyb2+1:ny+nyb2+NY);

% display
%
dx=1;              % x domain sample spacing
dy=1;              % y domain sample spacing
x=-NX/2:NX/2-1; y=-NY/2:NY/2-1;
%
G=abs(ref_image)';
xg=max(max(G)); ng=min(min(G)); cg=256/(xg-ng);
image(x,y,256-cg*(G-ng)); axis image;
xlabel('Spatial Domain X')
ylabel('Spatial Domain Y')
title('Reference Image f_1 (x,y)')
print P8.7.ps
pause(1)
%
G=abs(test_image)';
xg=max(max(G)); ng=min(min(G)); cg=256/(xg-ng);
image(x,y,256-cg*(G-ng)); axis image;
xlabel('Spatial Domain X')
ylabel('Spatial Domain Y')
title('Test Image f_2 (x,y)')
print P8.8.ps
pause(1)
%
xgg=max(max(abs(ref_image)));
I=(abs(ref_image) > .1*xgg);
gx=(HX.*I); gy=(HY.*I);
G=ref_image';
xg=max(max(G)); ng=min(min(G)); cg=256/(xg-ng);
image(x,y,256-cg*(G-ng)); axis image; axis xy
hold on
quiver(x(:)*ones(1,NY),ones(NX,1)*y,gx,gy,2)
xlabel('Spatial Domain X')
ylabel('Spatial Domain Y')
title('Motion Vector Image')
print P8.9.ps
hold off

author: Sable

repo: mcbench-benchmarks


text:
%> @brief Catkin-style basic dependencies from matlab_utilities
%> @param vararigin List of subdirectories to add to the Matlab path
%> Each dependency must be a subdirectory under this one.
function [] = addpath_thirdparty_packages(varargin)

% Get path of this file
base_path = fileparts(mfilename('fullpath'));

for i = 1:nargin
    dep = varargin{i};

%     if strcmp(dep,'yaml')
%         jcp = fullfile(prefdir,'javaclasspath.txt');
%         jcpID = fopen(jcp,'a');
%         fprintf(jcpID,...
%             [fullfile(pwd,'third','yaml','snakeyaml-1.10.jar')]...
%             );
%         fclose(jcpID);
%     end
    % Check then add path
    dep_path = fullfile(base_path, dep);
    assert(isfolder(dep_path), 'Invalid dependency (directory does not exist): %s', dep_path);
    addpath(dep_path);

    % Check for specific initialization
    %     init_path = fullfile(dep_path, [dep, '_init.m']);
    %     if exist(init_path, 'file')
    %         % Invoke initialization function
    %         feval([dep, '_init']);
    %     end
end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function piecewise_linear_product_integral_test ( )

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

  piecewise_linear_product_integral_test01 ( );
  piecewise_linear_product_integral_test02 ( );
  piecewise_linear_product_integral_test03 ( );
  piecewise_linear_product_integral_test04 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PIECEWISE_LINEAR_PRODUCT_INTEGRAL_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This test example is taken from demo_gabmulappr

% Setup parameters for the Gabor system and length of the signal
L=576; % Length of the signal
a=32;   % Time shift
M=72;  % Number of modulations
N=L/a;
fs=44100; % assumed sampling rate
SNRtv=63; % signal to noise ratio of change rate of time-variant system

% construction of slowly time variant system
% take an initial vector and multiply by random vector close to one
A = [];
c1=(1:L/2); c2=(L/2:-1:1); c=[c1 c2].^(-1); % weight of decay x^(-1)
A(1,:)=(tester_rand(1,L)-0.5).*c;  % convolution kernel
Nlvl = exp(-SNRtv/10);
Slvl = 1-Nlvl;
for ii=2:L;
     A(ii,:)=(Slvl*circshift(A(ii-1,:),[0 1]))+(Nlvl*(tester_rand(1,L)-0.5));
end;
A = A/norm(A)*0.99; % normalize matrix

% perform best approximation by gabor multiplier
g=gabtight(a,M,L);
sym1=gabmulappr(A,g,a,M);

% Now do the same using the general frame algorithm.

F=frame('dgt',g,a,M);

sym2=framemulappr(F,F,A);

norm(sym1-reshape(sym2,M,N))

% Test for exactness

testsym=tester_crand(M,N);
FT=frsynmatrix(F,L);

T=FT*diag(testsym(:))*FT';

sym1b=gabmulappr(T,g,a,M);
sym2b=framemulappr(F,F,T);

norm(testsym-sym1b)
norm(testsym-reshape(sym2b,M,N))

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function cdf = semicircular_cdf ( x, a, b )

%*****************************************************************************80
%
%% SEMICIRCULAR_CDF evaluates the Semicircular CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameter of the PDF.
%    0.0 < B.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= a - b )

    cdf = 0.0;

  elseif ( x <= a + b )

    y = ( x - a ) / b;

    cdf = 0.5 + ( y * sqrt ( 1.0 - y * y ) + asin ( y ) ) / pi;

  elseif ( a + b < x )

    cdf = 1.0;

  end

  return
end


text:
function [Rot] = AxisAngle2Rot(axisAngle)

    theta = norm(axisAngle, 2);

    nx = axisAngle / theta;

    nx = [    0  -nx(3)  nx(2);
           nx(3)     0  -nx(1);
          -nx(2)  nx(1)     0 ];

    Rot = eye(3) + sin(theta) * nx + (1-cos(theta))*nx^2;

author: TadasBaltrusaitis

repo: OpenFace


text:
function varargout = size(X,d)
% SIZE (overloaded)

if nargin == 1 && nargout > 1 && nargout == length(X.dim)
    for i = 1:length(X.dim);
        varargout{i} = X.dim(i);
    end
    return
elseif nargin == 1 && nargout > 1 && nargout < length(X.dim)
    % Weird but this is MATLAB standard
    for i = 1:nargout-1
        varargout{i} = X.dim(i);
    end
    varargout{end+1} = prod(X.dim(nargout:end));
    return
end

if nargin == 1
    varargout{1} = X.dim;
else
    if (d < 1) || d~=round(d)
        error('Dimension argument must be a positive integer scalar in the range 1 to 2^31.');
    end
    if d <= length(X.dim)
        varargout{1} = X.dim(d);
    else
        varargout{1} = 1;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% op_CSIRemoveLipids.m
%
% Removes lipids from CSI data using L2 regularization. This minimizes the
% equation norm(x - x_0, 2) + beta * norm(W'x, 2)
%
% INPUT:
% MRSIStruct        = MRSI structure used in FID-A
% lipidComponents   = number of lipid spectra in the lipid basis
% lineWidthRange    = range of linewidth used for building lipid basis
% ppmRange          = ppm range used in building lipid basis
% beta              = regularization term
% plotBasis         = plot basis spectra
%
%
% OUTPUT:
% MRSIStruct        = MRSI structure with lipids removed

function [MRSIStruct] = op_CSIRemoveLipids(MRSIStruct, basisArguments, plottingArguments)
    arguments
        MRSIStruct (1, 1) struct
        basisArguments.lipidComponenets (1, 1) double = 1000
        basisArguments.lineWidthRange (1, 2) double = [1 80]
        basisArguments.lipidPPMRange (1, 2) double = [0.3 1.9000]
        basisArguments.beta (1, 1) double = 1e-4
        plottingArguments.plotBasis (1, 1) logical = false
    end
    % extract arguments from name value pairs
    lipidComponents = basisArguments.lipidComponenets;
    lineWidthRange = basisArguments.lineWidthRange;
    lipidPPMRange = basisArguments.lipidPPMRange;
    beta = basisArguments.beta;

    spectraSize = getSizeFromDimensions(MRSIStruct, {'t'});

    % calculate lipid basis used for L2 regularization
    lipidBasis = createLipipBasis(MRSIStruct, lipidComponents, lineWidthRange, lipidPPMRange);
    if(plottingArguments.plotBasis)
        figure
        plot(MRSIStruct.ppm, flip(real(lipidBasis),1));
    end
    % calculate solution from the basis
    L2Solution = inv(eye(spectraSize) + beta * (lipidBasis * lipidBasis'));

    [MRSIStruct, prevPermute, prevShape] = reshapeDimensions(MRSIStruct, {'t'});
    data = getData(MRSIStruct);

    data = L2Solution * data;
    MRSIStruct = setData(MRSIStruct, data);
    MRSIStruct = reshapeBack(MRSIStruct, prevPermute, prevShape);
end

% calculate lipid basis
function lipidBasis = createLipipBasis(MRSIStruct, lipidComponents, lineWidthRange, lipidPPMRange)
    spectralWidth = getSpectralWidth(MRSIStruct);

    spectralPoints = getSizeFromDimensions(MRSIStruct, {'t'});
    fidBasis = zeros(spectralPoints, lipidComponents);

    lipidStructure = load('Lip.mat', 'sysLip');
    lipidStructure = lipidStructure.sysLip;

    for iSpectra = 1:lipidComponents

        fidBasis(:, iSpectra) = getRandomLipidFids(spectralPoints, spectralWidth, ...
                                                  lineWidthRange, lipidPPMRange, ...
                                                  lipidStructure);
    end
    lipidBasis = fftshift(fft(fidBasis, [], 1), 1);
end

% calculate the single lipid spectra for the basis
function lipidFids = getRandomLipidFids(spectralPoints, spectralWidth, lineWidth, ...
                                        lipidPPMRange, lipidSystem)
    [randomLineWidth, randomPPM] = getRandomLineWidthandPPM(lineWidth, lipidPPMRange);
    lipidSystem.shifts = randomPPM;

    simulatedSignal = sim_onepulse(spectralPoints, spectralWidth, 3, randomLineWidth, lipidSystem);
    simulatedSignal = op_complexConj(simulatedSignal);
    simulatedSignal = addRandomPhase(simulatedSignal);
    %simulatedSignal = scaleSpectra(simulatedSignal, randomPPM, lipidPPMRange);
    lipidFids = simulatedSignal.fids;
end

% pick a random number from lower bounds and upper bounds
function randomNumber = randomNumberInRange(lowerBounds, upperBounds)
    difference = upperBounds - lowerBounds;
    randomNumber = lowerBounds + difference * rand(1);
end

function simulatedSignal = addRandomPhase(simulatedSignal)
    sepctralPhase = randomNumberInRange(-180, 180);
    simulatedSignal = op_addphase(simulatedSignal, sepctralPhase, 0, 4.65, 1);
end

% scale spectra based on a normal distribution. Signal near the center of lipid
% range will be scaled high and signal near the edges scaled lower.
function simulatedSignal = scaleSpectra(simulatedSignal, lipidPPM, lipidRange)
    normalProbability = normpdf(lipidPPM, mean(lipidRange), diff(lipidRange)/4);
    simulatedSignal = op_ampScale(simulatedSignal, normalProbability);
    simulatedSignal = op_ampScale(simulatedSignal, 10);

end

function [randomLineWidth, randomPPM] = getRandomLineWidthandPPM(lineWidth, lipidPPMRange)
    randomLineWidth = randomNumberInRange(lineWidth(1), lineWidth(2));
    randomPPM = randomNumberInRange(lipidPPMRange(1), lipidPPMRange(2));
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function nP=spacing2numVertices(F,V,pointSpacing)

A=sum(patchArea(F,V)); %Total area
At=(pointSpacing.^2*sqrt(3))/4; %Theoretical area of equilateral triangle
NF=(A/At);

E=patchEdges(F,1); %Edges
nE=size(E,1); %Number of edges
nF=size(F,1); %Number of faces
nV=size(V,1); %Number of vertices

%Compute Euler characteristic
X=size(V,1)-size(E,1)+size(F,1);

nRefScalar=(log(NF)-log(nF))/log(4);

if nRefScalar<0
    nRef=floor(nRefScalar);
    nRange=0:-1:nRef;
else
    nRef=ceil(nRefScalar);
    nRange=0:1:nRef;
end

nvR=nV.*ones(numel(nRange),1);
neR=nE.*ones(numel(nRange),1);

nfR=nF*4.^nRange';
for q=2:1:numel(nRange)
    if nRef>0
        nvR(q)=nvR(q-1)+neR(q-1);
    elseif nRef<0
        nvR(q)=(nvR(q-1)+X-nfR(q))/2;
    end
    neR(q)=-X+nfR(q)+nvR(q);
end

% [nvR -neR nfR]
% sum([nvR -neR nfR],2)

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON