text:
function geometry_test02015 ( )

%*****************************************************************************80
%
%% TEST02015 tests CYLINDER_POINT_DIST_SIGNED_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  test_num = 6;

  dist_test = [ 3.0, -0.5, 5.0, 8.0, -1.0, -0.25 ];
  p_test = [ ...
      4.0,   0.5,  0.0; ...
     -0.5,  -1.0,  0.0; ...
      4.0,   6.0,  0.0; ...
      0.75, -10.0, 0.0; ...
      0.0,   0.0,  0.0; ...
      0.25,  1.75, 0.0 ]';
  p1 = [ 0.0, -2.0, 0.0 ];
  p2 = [ 0.0,  2.0, 0.0 ];
  r = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02015\n' );
  fprintf ( 1, '  CYLINDER_POINT_DIST_SIGNED_3D computes the signed\n' );
  fprintf ( 1, '  distance to a cylinder.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Radius R = %f\n', r );
  fprintf ( 1, '  Center of bottom disk = %10f  %10f  %10f\n', p1(1:dim_num) );
  fprintf ( 1, '  Center of top disk =    %10f  %10f  %10f\n', p2(1:dim_num) );

  for test = 1 : test_num

    p(1:dim_num) = p_test(1:dim_num,test);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  P =    %10f  %10f  %10f\n', p(1:dim_num) );

    dist = cylinder_point_dist_signed_3d ( p1, p2, r, p );

    fprintf ( 1, '  Signed distance (computed) = %f\n', dist );
    fprintf ( 1, '  Signed distance (exact) =    %f\n', dist_test(test) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ n_data, x, fx ] = si_values ( n_data )

%*****************************************************************************80
%
%% SI_VALUES returns some values of the sine integral function.
%
%  Discussion:
%
%    SI(X) = integral ( 0 <= T <= X ) sin ( T ) / T dt
%
%    In Mathematica, the function can be evaluated by:
%
%      SinIntegral[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 16;

  fx_vec = [ ...
     0.4931074180430667E+00, ...
     0.5881288096080801E+00, ...
     0.6812222391166113E+00, ...
     0.7720957854819966E+00, ...
     0.8604707107452929E+00, ...
     0.9460830703671830E+00, ...
     0.1108047199013719E+01, ...
     0.1256226732779218E+01, ...
     0.1389180485870438E+01, ...
     0.1505816780255579E+01, ...
     0.1605412976802695E+01, ...
     0.1778520173443827E+01, ...
     0.1848652527999468E+01, ...
     0.1833125398665997E+01, ...
     0.1758203138949053E+01, ...
     0.1654140414379244E+01 ];

  x_vec = [ ...
      0.5E+00, ...
      0.6E+00, ...
      0.7E+00, ...
      0.8E+00, ...
      0.9E+00, ...
      1.0E+00, ...
      1.2E+00, ...
      1.4E+00, ...
      1.6E+00, ...
      1.8E+00, ...
      2.0E+00, ...
      2.5E+00, ...
      3.0E+00, ...
      3.5E+00, ...
      4.0E+00, ...
      4.5E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function Population = EnvironmentSelection(Population,N)
% The environmental selection of NSGA-II

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Non-dominated sorting
    PopObj = Population.objs;
    [FrontNo,MaxFNo] = NDSort(PopObj,N);
    Next   = false(1,length(FrontNo));
    Next(FrontNo<MaxFNo) = true;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(PopObj,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [A_lh, A_rh]=cortical_label_adf(subject, p_val)
% Computes the area of the different cortical labels
% and compare them to the normal range
% Uses p_value to detect the abnormal areas
% Uses the lh/rh.parc.txt files
%

%
% cortical_labeling_afd_txt.m
%
% Original Author: Laurence Wastiaux
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

if (nargin<2 | nargin>2)
    msg=sprintf('USAGE: [L, R]=cortical_label_adf(subject, p_val)');
    disp(msg)
end

A_lh=zeros(1,84);
A_rh=zeros(1,84);
namerh=[];
namelh=[];

labmapfile=('/space/lyon/1/fsdev/freesurfer_dev/Simple_surface_labels2002.txt');
[label name val1 val2 val3 val4]=textread(labmapfile,'%d %s %d %d %d %d',85);

%% a- Check if the left and right stats files exist %%
LabelDir=strcat(subject, '/label');
dlabel=dir(LabelDir);
if(length(dlabel)<=2)
    mess=sprintf('Directory %s is empty', LabelDir);
    error(mess)
end
if(length(dlabel)>2)
    if ((exist(strcat(LabelDir, '/rh.parc.txt'))==0) | (exist(strcat(LabelDir, '/lh.parc.txt'))==0))
        mess=sprintf('Cannot find the stats files for subject %s', subject);
        error(mess)
    end
end

%% b- Read areas for each label & for both hemispheres using the lh/rh.parc.txt files %%
sn=strread(subject, '%s', 'delimiter', '/');
subject_name=char(sn(length(sn)));
rh_LabelFile=strcat(subject, '/label/rh.parc.txt');
lh_LabelFile=strcat(subject, '/label/lh.parc.txt');

lh_fid=fopen(lh_LabelFile);
rh_fid=fopen(rh_LabelFile);
if(lh_fid == -1)
    msg=sprintf('Cannot open file %s', lh_LabelFile);
    error(msg)
elseif(rh_fid ==-1)
    msg=sprintf('Cannot open file %s', rh_LabelFile);
    error(msg)
else
    while(feof(lh_fid)==0)
        s=fgetl(lh_fid);
        if(strfind(s,'structure name')>1)
            s=fgetl(lh_fid);
            lhpt=ftell(lh_fid);
            break;
        end
    end
    fseek(lh_fid, lhpt,'bof');
    while(feof(lh_fid)==0)
        s=fgetl(lh_fid);
        a=strread(s,'%d',2);  %area
        aa=strread(s,'%s',11); %name
        llab=0;
        for test=1:length(name)
            if(strfind(char(name(test)),char(aa(11)))>=1)
                llab=label(test);
                break;
            end
        end
        if(llab==0)
            disp(sprintf('impossible to find the label'));
        elseif(A_lh(llab)~=0)
            disp(sprintf('already attributed'));
        else
            A_lh(llab)=a(2);
        end
    end
    while(feof(rh_fid)==0)
        s=fgetl(rh_fid);
        if(strfind(s,'structure name')>1)
            s=fgetl(rh_fid);
            rhpt=ftell(rh_fid);
            break;
        end
    end
    fseek(rh_fid, rhpt,'bof');
    while(feof(rh_fid)==0)
        s=fgetl(rh_fid);
        a=strread(s,'%d',2);
        aa=strread(s,'%s',11);
        llab=0;
        for test=1:length(name)
            if(strfind(char(name(test)),char(aa(11)))>=1)
                llab=label(test);
                break;
            end
        end
        if(llab==0)
            disp(sprintf('impossible to find the label'));
        elseif(A_rh(llab)~=0)
            disp(sprintf('already attributed'));
        else
            A_rh(llab)=a(2);
        end
    end
end
nl=length(name)-1;
if(length(A_lh)< nl)
    for lleft=(length(A_lh)+1):nl
        A_lh(lleft)=0;
    end
end
if(length(A_rh)< nl)
    for lleft=(length(A_rh)+1):nl
        A_rh(lleft)=0;
    end
end
if (sum(A_lh)==0 | sum(A_rh)==0)
    msg=sprintf('Total area is null');
    error(msg)
else
    total_lharea=sum(A_lh);
    total_rharea=sum(A_rh);
end
for i=1:nl
    A_lh(i)= 100 * A_lh(i) / total_lharea;
    A_rh(i)= 100 * A_rh(i) / total_rharea;
end
% Compute p_values %

count=0;
for k =1:length(A_lh)
    %disp(k)

    [Rpval_inf,Rpval_sup, Lpval_inf, Lpval_sup]=compute_pval(k,A_rh(k),A_lh(k));
    if ((Lpval_inf<p_val | Lpval_sup <p_val) & (A_lh(k)~=0))
        fprintf('%s, left hemisphere: percent area of %s (%d) = %g (pval_inf=%g, pval_sup=%g)\n' ,subject_name, char(name(k+1)), k , A_lh(k),Lpval_inf, Lpval_sup);
        count=1;
    end
    if ((Rpval_inf<p_val | Rpval_sup <p_val)& (A_rh(k)~=0))
        fprintf('%s, right hemisphere: percent area of %s (%d) = %g (pval_inf=%g, pval_sup=%g)\n' ,subject_name, char(name(k+1)),k, A_rh(k),Rpval_inf,Rpval_sup);
        count=1;
    end

end
if(count==0)
    fprintf('%s: subject is normal\n', subject_name);
end

% subfunction compute_pval() %
function [pinf_rh, psup_rh, pinf_lh, psup_lh]=compute_pval(llabel,R,L)
%load('/space/okapi/3/data/laurence/ADF/cortical_labeling/PercentArea_labels.mat'); % loads D_lh and D_rh
% load('/space/okapi/3/data/laurence/ADF/cortical_labeling/PercentArea_labels_parctxt.mat');
% D_lh=D2_lh;
% D_rh=D2_rh;
%rh_stat_file='/space/okapi/3/data/laurence/ADF/cortical_labeling/rh.CorticalLabelingPercentArea_txt.adf';
%lh_stat_file='/space/okapi/3/data/laurence/ADF/cortical_labeling/lh.CorticalLabelingPercentArea_txt.adf';
%%% Get the table's directory %%%
if(getenv('FREESURFER_HOME'))
    fsh=getenv('FREESURFER_HOME');
    fsafdDir=strcat(fsh, '/fsafd');
else
    error(sprintf('Impossible to find FREESURFER_HOME\n'));
end
rh_stat_file=strcat(fsafdDir, '/rh.CorticalLabelingPercentArea_txt.adf');
lh_stat_file=strcat(fsafdDir, '/lh.CorticalLabelingPercentArea_txt.adf');
fidrh=fopen(rh_stat_file);
fidlh=fopen(lh_stat_file);
if(fidrh==-1 | fidlh==-1)
    mess=sprintf('Could not find %s or %s', rh_stat_file, lh_stat_file);
    error(mess)
end
while(strfind(fgetl(fidrh), '#'))
    pos=ftell(fidrh);
end
fseek(fidrh, pos, 'bof');
D_rhtmp=fscanf(fidrh, '%g');
nrowr=length(D_rhtmp)/84;
D_rh=(reshape(D_rhtmp, [84, nrowr]))';
while(strfind(fgetl(fidlh), '#'))
    pos=ftell(fidlh);
end
fseek(fidlh, pos, 'bof');
D_lhtmp=fscanf(fidlh, '%g');
nrow=length(D_lhtmp)/84;
D_lh=(reshape(D_lhtmp, [84,nrow]))';
%distribution%
pas=0.01;
x=0:pas:18;
[h1] = hist(D_lh(:,llabel),x);
[h2] = hist(D_rh(:,llabel),x);
%[h1] = hist(D2_lh(:,llabel),x);
%[h2] = hist(D2_rh(:,llabel),x);
p1 = h1/sum(h1);
p2 = h2/sum(h2);
d1inf=find(x<=L);
d2inf=find(x<=R);
d1sup=find(x>L);
d2sup=find(x>R);
x1inf=x(d1inf);
x1sup=x(d1sup);
x2inf=x(d2inf);
x2sup=x(d2sup);
p1inf=p1(1:length(x1inf));
p2inf=p2(1:length(x2inf));
p1sup=p1(length(x)-length(x1sup)+1:end);
p2sup=p2(length(x)-length(x2sup)+1:end);
if( L>=0 & length(x1inf)>1 & length(x1sup)>1)
    pinf_lh=trapz(x1inf,p1inf)/pas;
    psup_lh=trapz(x1sup,p1sup)/pas;
else   % No subdivisions of the last intervals
    pinf_lh=0;
    psup_lh=0;
end
if(R>=0 & length(x2inf)>1 & length(x2sup)>1)
    pinf_rh=trapz(x2inf,p2inf)/pas;
    psup_rh=trapz(x2sup,p2sup)/pas;
else
    pinf_rh=0;
    psup_rh=0;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = boundedCentroidalVoronoi2d(germs, box, varargin)
%BOUNDEDCENTROIDALVORONOI2D Create a 2D Centroidal Voronoi Tesselation in a box.
%
%   [N, E, F] = boundedCentroidalVoronoi2d(GERMS, BOX)
%   GERMS are N-by-2 point array, BOX is given as [xmin xmax ymin ymax].
%   Algorithm is an iteration of voronoi diagram computations, using at
%   each steps the centroids of previous diagram as germs for the new
%   diagram.
%
%   [N, E, F] = boundedCentroidalVoronoi2d(GERMS, BOX, NITER)
%   Specifies the number of iterations. Default is 10.
%
%   [N, E, F, G] = boundedCentroidalVoronoi2d(...)
%   also returns the positions of germs/centroids for each face. If the
%   number of iteration was sufficient, location of germs should correspond
%   to centroids of faces 'fc' computed using:
%   fc(i,:) = polygonCentroid(n(f{i}, :));
%
%   Example
%   [n, e, f] = boundedCentroidalVoronoi2d(rand(20, 2)*100, [0 100 0 100]);
%   drawGraph(n, e, f);
%
%   See also
%     graphs, boundedVoronoi2d, centroidalVoronoi2d, clipGraph
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2007-01-12
% Copyright 2007-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

% number of iteration
nIter = 10;
if ~isempty(varargin)
    nIter = varargin{1};
end

% limits and size of the box
x0 = box(1); x1 = box(2);
y0 = box(3); y1 = box(4);
dx = x1 - x0;  dy = y1 - y0;

% far points to bound the voronoi diagram
farPoints = [...
    x1 + 10 * dx, y1 + 10 * dy;...
    x0 - 10 * dx, y1 + 10 * dy;...
    x0 - 10 * dx, y0 - 10 * dy;...
    x1 + 10 * dx, y0 - 10 * dy];

% iterate bounded voronoi tesselation
for i = 1:nIter
    % generate Voronoi diagram, and clip with the box
    [n, e, f] = voronoi2d([germs ; farPoints]);
    [n, e, f] = clipGraph(n, e, f, box);

    % centroid of each face will be used as germs for next iteration
    for j = 1:length(f)
        face = n(f{j}, :);
        germs(j, 1:2) = polygonCentroid(face);
    end
end

% result is given in n, e, and f, eventually germs
varargout{1} = n;
varargout{2} = e;
varargout{3} = f;
if nargout > 3
    varargout{4} = germs;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function computeClinicalResults_RF(pathRF,nBoot,nSplit,testSplit,testCost,seeds,matlabPATH)

startpath = pwd;

cd(pathRF), load('training'), load('testing')
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
comb =  {{[1,2,3],...      % Loco        : Age, Subtype, T_Stage
        [1,2,4],...        % Loco        : Age, Subtype, N_Stage
        [1,2,5],...        % Loco        : Age, Subtype, TNM_Stage
        [1,2,3,4]},...     % Loco        : Age, Subtype, T_Stage, N_Stage
        {[1,2,3],...       % Distant     : Age, Subtype, T_Stage
        [1,2,4],...        % Distant     : Age, Subtype, N_Stage
        [1,2,5],...        % Distant     : Age, Subtype, TNM_Stage
        [1,2,3,4]},...     % Distant     : Age, Subtype, T_Stage, N_Stage
        {[1,2,3],...       % Death       : Age, Subtype, T_Stage
        [1,2,4],...        % Death       : Age, Subtype, N_Stage
        [1,2,5],...        % Death       : Age, Subtype, TNM_Stage
        [1,2,3,4]}};       % Death       : Age, Subtype, T_Stage, N_Stage
mkdir('batchLog_RFclinical'), cd('batchLog_RFclinical'), pathBatch = pwd;
time = 60; % Time in the wainting loop

% PRODUCE BATCH COMPUTATIONS TO FIND THE BEST CLINICAL VARIABLES
seed = seeds(1); batch = 0;
save('workspace','pathRF','nameOutcomes','comb','testCost','nSplit','testSplit','nBoot','seed'), pause(2);
for o = 1:nOutcomes
    batch = batch + 1;
    nameScript = ['batch',num2str(batch),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    fprintf(fid,['[indClinic,bestCost] = estimateClinicalRF(pathRF,nameOutcomes{',num2str(o),'},comb{',num2str(o),'},testCost,nSplit,testSplit,nBoot,seed);\n']);
    fprintf(fid,['save(''result_batch',num2str(batch),''',''indClinic'',''bestCost'')\n']);
    fprintf(fid,['system(''touch batch',num2str(batch),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nOutcomes)
delete('workspace.mat')

% GROUP RESULTS
cd(pathBatch)
batch = 0; parameters = struct;
for o = 1:nOutcomes
    batch = batch + 1;
    load(['result_batch',num2str(batch)]) % Variables 'indClinic' and 'bestCost' gets out of there
    parameters.clinical.(nameOutcomes{o}) = indClinic;
    parameters.cost.(nameOutcomes{o}) = bestCost;
    delete(['result_batch',num2str(batch),'.mat'])
    clear indClinic bestCost
end
cd(pathRF), save('clinicalPARAM','parameters')

% COMPUTE THE RANDOM FORESTS
for o = 1:nOutcomes
    if strcmp(nameOutcomes{o},'DeathSign')
        outcome = training.outcomes.Death;
    else
        outcome = training.outcomes.(nameOutcomes{o});
    end
    indClinic = parameters.clinical.(nameOutcomes{o});
    cost = parameters.cost.(nameOutcomes{o});
    tableTrain = training.clinical.table(:,indClinic);
    cat = logical(training.clinical.categories(indClinic));
    rng(seeds(2)), [RF] = trainRF_table(tableTrain,outcome,cat,nBoot,cost);
    RF = compact(RF); % Compact version
    save(['RF_clinic_',nameOutcomes{o}],'RF')
end

% TEST THE RANDOM FORESTS
cd(pathRF), load('clinicalPARAM')
for o = 1:nOutcomes
    if strcmp(nameOutcomes{o},'DeathSign')
        outcome = testing.outcomes.Death;
        time = testing.timeToEvents.Death;
    else
        outcome = testing.outcomes.(nameOutcomes{o});
        time = testing.timeToEvents.(nameOutcomes{o});
    end
    censoring = 1 - outcome;
    results = struct;
    RF = load(['RF_clinic_',nameOutcomes{o}]); RF = struct2cell(RF); RF = RF{1};
    indClinic = parameters.clinical.(nameOutcomes{o});
    tableTest = testing.clinical.table(:,indClinic);
    [prob] = predictRF(tableTest,RF);
    results.probResponse = prob;
    [AUC,sensitivity,specificity,accuracy] = calcPerformMetrics(prob,outcome,0.5);
    CI = calcCI(prob,time,censoring);
    results.AUC = AUC; results.Sensitivity = sensitivity; results.Specificity = specificity; results.Accuracy = accuracy;
    results.CI = CI;
    save(['testResultsRF_clinic_',nameOutcomes{o}],'results'), clear results
end

cd(startpath)
end

author: mvallieres

repo: radiomics


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

% relationship of linear convolution with circular convolution

% 4-point circular convolution
x1=[1,2,3,4];
x2=[3,2,5,1];
N=length(x1);
for m=0:N-1
p(m+1)=mod(-m,N);
x2s(1+m)=x2(1+p(m+1));
end
x2s
for n=0:N-1
    x2sn=circshift(x2s',n);
    y2(n+1)=x1*x2sn;
end
y2

% 7-point circular convolution
x11=[ x1 0 0 0];
x22=[x2 0 0 0];
N=length(x11);
for m=0:N-1
p(m+1)=mod(-m,N);
end
for m=0:N-1
x22s(1+m)=x22(1+p(m+1));
end
for n=0:N-1
    x22sn=circshift(x22s',n);
    y22(n+1)=x11*x22sn;
end
y22

%linear convolution
y1=conv(x1,x2)

author: Sable

repo: mcbench-benchmarks


text:
function [W, stat] = gsp_learn_graph_log_degrees(Z, a, b, params)
%GSP_LEARN_GRAPH_LOG_DEGREES Learn graph from pairwise distances using negative log prior on nodes degrees
%   Usage:  [W, stat] = gsp_learn_graph_log_degrees(Z, a, b)
%           [W, stat] = gsp_learn_graph_log_degrees(Z, a, b, params)
%
%   Inputs:
%         Z         : Matrix with (squared) pairwise distances of nodes
%         a         : Log prior constant  (bigger a -> bigger weights in W)
%         b         : ||W||_F^2 prior constant  (bigger b -> more dense W)
%         params    : Optional parameters
%
%   Outputs:
%         W         : Weighted adjacency matrix
%         stat      : Optional output statistics (adds small overhead)
%
%
%   'W = gsp_learn_graph_log_degrees(Z, a, b, params)' computes a weighted
%   adjacency matrix $W$ from squared pairwise distances in $Z$, using the
%   smoothness assumption that $\text{trace}(X^TLX)$ is small, where $X$ is
%   the data (columns) changing smoothly from node to node on the graph and
%   $L = D-W$ is the combinatorial graph Laplacian. See the paper of the
%   references for the theory behind the algorithm.
%
%   Alternatively, Z can contain other types of distances and use the
%   smoothness assumption that
%
%   .. sum(sum(W .* Z))
%
%   .. math:: \sum_i\sum_j W_{ij}Z_{ij}
%
%   is small.
%
%   The minimization problem solved is
%
%   .. minimize_W sum(sum(W .* Z)) - a * sum(log(sum(W))) + b * ||W||_F^2/2 + c * ||W-W_0||_F^2/2
%
%   .. math:: \min_W \sum_i\sum_j W_{ij}Z_{ij} -\alpha\sum_i\log(\sum_jW_{ij}) + \frac{\beta}{2} \|W\|_F^2 + c/2 ||W - W0||_F^2
%
%   subject to $W$ being a valid weighted adjacency matrix (non-negative,
%   symmetric, with zero diagonal).
%
%   The algorithm used is forward-backward-forward (FBF) based primal dual
%   optimization (see references).
%
%   Example:::
%
%         G = gsp_sensor(256);
%         f1 = @(x,y) sin((2-x-y).^2);
%         f2 = @(x,y) cos((x+y).^2);
%         f3 = @(x,y) (x-.5).^2 + (y-.5).^3 + x - y;
%         f4 = @(x,y) sin(3*((x-.5).^2+(y-.5).^2));
%         X = [f1(G.coords(:,1), G.coords(:,2)), f2(G.coords(:,1), G.coords(:,2)), f3(G.coords(:,1), G.coords(:,2)), f4(G.coords(:,1), G.coords(:,2))];
%         figure; subplot(2,2,1); gsp_plot_signal(G, X(:,1)); title('1st smooth signal');
%         subplot(2,2,2); gsp_plot_signal(G, X(:,2)); title('2nd smooth signal');
%         subplot(2,2,3); gsp_plot_signal(G, X(:,3)); title('3rd smooth signal');
%         subplot(2,2,4); gsp_plot_signal(G, X(:,4)); title('4th smooth signal');
%         Z = gsp_distanz(X').^2;
%         % we can multiply the pairwise distances with a number to control sparsity
%         [W] = gsp_learn_graph_log_degrees(Z*25, 1, 1);
%         % clean up zeros
%         W(W<1e-5) = 0;
%         G2 = gsp_update_weights(G, W);
%         figure; gsp_plot_graph(G2); title('Graph with edges learned from above 4 signals');
%
%
%   Additional parameters
%   ---------------------
%
%   * *params.W_init*   : Initialization point. default: zeros(size(Z))
%   * *verbosity*       : Default = 1. Above 1 adds a small overhead
%   * *maxit*           : Maximum number of iterations. Default: 1000
%   * *tol*             : Tolerance for stopping criterion. Defaul: 1e-5
%   * *step_size*       : Step size from the interval (0,1). Default: 0.5
%   * *max_w*           : Maximum weight allowed for each edge (or inf)
%   * *w_0*             : Vector for adding prior c/2*||w - w_0||^2
%   * *c*               : multiplier for prior c/2*||w - w_0||^2 if w_0 given
%   * *fix_zeros*       : Fix a set of edges to zero (true/false)
%   * *edge_mask*       : Mask indicating the non zero edges if "fix_zeros"
%
%   If fix_zeros is set, an edge_mask is needed. Only the edges
%   corresponding to the non-zero values in edge_mask will be learnt. This
%   has two applications: (1) for large scale applications it is cheaper to
%   learn a subset of edges. (2) for some applications we don't want some
%   connections to be allowed, for example for locality on images.
%
%   The cost of each iteration is linear to the number of edges to be
%   learned, or the square of the number of nodes (numel(Z)) if fix_zeros
%   is not set.
%
%   The function is using the UNLocBoX functions sum_squareform and
%   squareform_sp.
%   The stopping criterion is whether both relative primal and dual
%   distance between two iterations are below a given tolerance.
%
%   To set the step size use the following rule of thumb: Set it so that
%   relative change of primal and dual converge with similar rates (use
%   verbosity > 1).
%
%   See also: gsp_learn_graph_l2_degrees gsp_distanz gsp_update_weights
%       squareform_sp sum_squareform gsp_compute_graph_learning_theta
%
%   References: kalofolias2016learn komodakis2015playing kalofolias2017large
%

% Author: Vassilis Kalofolias
% Testing: gsp_test_learn_graph
% Date: June 2015

%% Default parameters
if nargin < 4
    params = struct;
end

if not(isfield(params, 'verbosity')),   params.verbosity = 1;   end
if not(isfield(params, 'maxit')),       params.maxit = 1000;      end
if not(isfield(params, 'tol')),         params.tol = 1e-5;      end
if not(isfield(params, 'step_size')),   params.step_size = .5;      end     % from (0, 1)
if not(isfield(params, 'fix_zeros')),   params.fix_zeros = false;      end
if not(isfield(params, 'max_w')),       params.max_w = inf;         end

%% Fix parameter size and initialize
if isvector(Z)
    z = Z;  % lazy copying of matlab doesn't allocate new memory for z
else
    z = squareform_sp(Z);
end
% clear Z   % for large scale computation

z = z(:);
l = length(z);                      % number of edges
% n(n-1)/2 = l => n = (1 + sqrt(1+8*l))/ 2
n = round((1 + sqrt(1+8*l))/ 2);    % number of nodes

if isfield(params, 'w_0')
    if not(isfield(params, 'c'))
        error('When params.w_0 is specified, params.c should also be specified');
    else
        c = params.c;
    end
    if isvector(params.w_0)
        w_0 = params.w_0;
    else
        w_0 = squareform_sp(params.w_0);
    end
    w_0 = w_0(:);
else
    w_0 = 0;
end

% if sparsity pattern is fixed we optimize with respect to a smaller number
% of variables, all included in w
if params.fix_zeros
    if not(isvector(params.edge_mask))
        params.edge_mask = squareform_sp(params.edge_mask);
    end
    % use only the non-zero elements to optimize
    ind = find(params.edge_mask(:));
    z = full(z(ind));
    if not(isscalar(w_0))
        w_0 = full(w_0(ind));
    end
else
    z = full(z);
    w_0 = full(w_0);
end

w = zeros(size(z));

%% Needed operators
% S*w = sum(W)
if params.fix_zeros
    [S, St] = sum_squareform(n, params.edge_mask(:));
else
    [S, St] = sum_squareform(n);
end

% S: edges -> nodes
K_op = @(w) S*w;

% S': nodes -> edges
Kt_op = @(z) St*z;

if params.fix_zeros
    norm_K = normest(S);
    % approximation:
    % sqrt(2*(n-1)) * sqrt(nnz(params.edge_mask) / (n*(n+1)/2)) /sqrt(2)
else
    % the next is an upper bound if params.fix_zeros
    norm_K = sqrt(2*(n-1));
end

%% TODO: Rescaling??
% we want    h.beta == norm_K   (see definition of mu)
% we can multiply all terms by s = norm_K/2*b so new h.beta==2*b*s==norm_K

%% Learn the graph
% min_{W>=0}     tr(X'*L*X) - gc * sum(log(sum(W))) + gp * norm(W-W0,'fro')^2, where L = diag(sum(W))-W
% min_W       I{W>=0} + W(:)'*Dx(:)  - gc * sum(log(sum(W))) + gp * norm(W-W0,'fro')^2
% min_W                f(W)          +       g(L_op(W))      +   h(W)

% put proximal of trace plus positivity together
f.eval = @(w) 2*w'*z;    % half should be counted
%f.eval = @(W) 0;
f.prox = @(w, c) min(params.max_w, max(0, w - 2*c*z));  % all change the same

param_prox_log.verbose = params.verbosity - 3;
g.eval = @(z) -a * sum(log(z));
g.prox = @(z, c) prox_sum_log(z, c*a, param_prox_log);
% proximal of conjugate of g: z-c*g.prox(z/c, 1/c)
g_star_prox = @(z, c) z - c*a * prox_sum_log(z/(c*a), 1/(c*a), param_prox_log);

if w_0 == 0
    % "if" not needed, for c = 0 both are the same but gains speed
    h.eval = @(w) b * norm(w)^2;
    h.grad = @(w) 2 * b * w;
    h.beta = 2 * b;
else
    h.eval = @(w) b * norm(w)^2 + c * norm(w - w_0,'fro')^2;
    h.grad = @(w) 2 * ((b+c) * w - c * w_0);
    h.beta = 2 * (b+c);
end

%% My custom FBF based primal dual (see [1] = [Komodakis, Pesquet])
% parameters mu, epsilon for convergence (see [1])
mu = h.beta + norm_K;     %TODO: is it squared or not??
epsilon = lin_map(0.0, [0, 1/(1+mu)], [0,1]);   % in (0, 1/(1+mu) )

% INITIALIZATION
% primal variable ALREADY INITIALIZED
%w = params.w_init;
% dual variable
v_n = K_op(w);
if nargout > 1 || params.verbosity > 1
    stat.f_eval = nan(params.maxit, 1);
    stat.g_eval = nan(params.maxit, 1);
    stat.h_eval = nan(params.maxit, 1);
    stat.fgh_eval = nan(params.maxit, 1);
    stat.pos_violation = nan(params.maxit, 1);
end
if params.verbosity > 1
    fprintf('Relative change of primal, dual variables, and objective fun\n');
end

tic
gn = lin_map(params.step_size, [epsilon, (1-epsilon)/mu], [0,1]);              % in [epsilon, (1-epsilon)/mu]
for i = 1:params.maxit
    Y_n = w - gn * (h.grad(w) + Kt_op(v_n));
    y_n = v_n + gn * (K_op(w));
    P_n = f.prox(Y_n, gn);
    p_n = g_star_prox(y_n, gn); % = y_n - gn*g_prox(y_n/gn, 1/gn)
    Q_n = P_n - gn * (h.grad(P_n) + Kt_op(p_n));
    q_n = p_n + gn * (K_op(P_n));

    if nargout > 1 || params.verbosity > 2
        stat.f_eval(i) = f.eval(w);
        stat.g_eval(i) = g.eval(K_op(w));
        stat.h_eval(i) = h.eval(w);
        stat.fgh_eval(i) = stat.f_eval(i) + stat.g_eval(i) + stat.h_eval(i);
        stat.pos_violation(i) = -sum(min(0,w));
    end
    rel_norm_primal = norm(- Y_n + Q_n, 'fro')/norm(w, 'fro');
    rel_norm_dual = norm(- y_n + q_n)/norm(v_n);

    if params.verbosity > 3
        fprintf('iter %4d: %6.4e   %6.4e   %6.3e', i, rel_norm_primal, rel_norm_dual, stat.fgh_eval(i));
    elseif params.verbosity > 2
        fprintf('iter %4d: %6.4e   %6.4e   %6.3e\n', i, rel_norm_primal, rel_norm_dual, stat.fgh_eval(i));
    elseif params.verbosity > 1
        fprintf('iter %4d: %6.4e   %6.4e\n', i, rel_norm_primal, rel_norm_dual);
    end

    w = w - Y_n + Q_n;
    v_n = v_n - y_n + q_n;

    if rel_norm_primal < params.tol && rel_norm_dual < params.tol
        break
    end
end
stat.time = toc;
if params.verbosity > 0
    fprintf('# iters: %4d. Rel primal: %6.4e Rel dual: %6.4e  OBJ %6.3e\n', i, rel_norm_primal, rel_norm_dual, f.eval(w) + g.eval(K_op(w)) + h.eval(w));
    fprintf('Time needed is %f seconds\n', stat.time);
end

% use the following for testing:
% g.L = K_op;
% g.Lt = Kt_op;
% g.norm_L = norm_K;
% [w, info] = fbf_primal_dual(w, f, g, h, params);
% %[w, info] = fb_based_primal_dual(w, f, g, h, params);

%%

if params.verbosity > 3
    figure; plot(real([stat.f_eval, stat.g_eval, stat.h_eval])); hold all; plot(real(stat.fgh_eval), '.'); legend('f', 'g', 'h', 'f+g+h');
    figure; plot(stat.pos_violation); title('sum of negative (invalid) values per iteration')
    figure; semilogy(max(0,-diff(real(stat.fgh_eval'))),'b.-'); hold on; semilogy(max(0,diff(real(stat.fgh_eval'))),'ro-'); title('|f(i)-f(i-1)|'); legend('going down','going up');
end

if params.fix_zeros
    w = sparse(ind, ones(size(ind)), w, l, 1);
end

if isvector(Z)
    W = w;
else
    W = squareform_sp(w);
end

repo: gspbox


text:
function drawStopActionAcrobot(soln,p)

clf; hold on;

length = p.l1+p.l2;
axis equal; axis(length*[-1,1,-1,1]); axis off;

% Start by drawing the trajectory of the tip
nCurve = 250;
t = linspace(soln.grid.time(1), soln.grid.time(end), nCurve);
Z = soln.interp.state(t);
[~,p2] = acrobotKinematics(Z,p);
plot(p2(1,:),p2(2,:),'LineWidth',2,'Color',[0.2,0.2,0.8]);

% Now draw the stop-action frames
nFrame = 10;
t = linspace(soln.grid.time(1), soln.grid.time(end), nFrame);
Z = soln.interp.state(t);

for i=1:nFrame
    z = Z(:,i);

    [p1,p2] = acrobotKinematics(z,p);
    pos = [[0;0],p1,p2];

    val = 0.3 + 0.7*(i-1)/(nFrame-1);
    color = hsv2rgb([0.0, 0.8, val]);

    plot(0,0,'ks','MarkerSize',25,'LineWidth',3)
    plot(pos(1,:),pos(2,:),'Color',color,'LineWidth',3)
    plot(pos(1,1),pos(2,1),'k.','MarkerSize',20)
    plot(pos(1,2),pos(2,2),'k.','MarkerSize',20)
    plot(pos(1,3),pos(2,3),'k.','MarkerSize',40)

end

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function sL = Laue(s)
% return the corresponding Laue group

% maybe we have alread computed the Laue group
if ~isempty(s.LaueRef)

  sL = s.LaueRef;

elseif s.isLaue % if it is already a Laue group then there is nothing to do

  sL = s;

else % in the meantime isLaue has computed the Laue group :)

  sL = s.LaueRef;

end

author: mtex-toolbox

repo: mtex

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


text:
%EVALINIT Calculates initial solution vector from given expressions.
%
%   [ U0 ] = EVALINIT( PROB, C_EXPR ) Given a valid PROB struct
%   evaluates the inital value expressions in C_EXPR (or PHYS struct
%   if C_EXPR is not given), and returns a solution vector U0.
%
%   See also EVALEXPR

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function tests = RegionFeatureTest
    tests = functiontests(localfunctions);
    clc
end

function teardownOnce(tc)
    close all
end

function thresh_test(tc)
    castle = iread('castle.png', 'double', 'grey');

    ithresh(castle);
    t = otsu(castle);
    t = niblack(castle, -0.2, 35);
    tc.verifyEqual(size(t), size(castle));

    castle = iread('castle.png', 'grey');
    ithresh(castle);
    t = otsu(castle);
    tc.verifyTrue( t > 0 && t < 255);
    t = niblack(castle, -0.2, 35);
    tc.verifyEqual(size(t), size(castle));
end

function colorseg_test(tc)

    im = iread('yellowtargets.png', 'gamma', 'sRGB', 'double');
    % kmeans clustering
    K = 4;
    [cls, cxy,resid] = colorkmeans(im, K);

    tc.verifyEqual(size(cls), size(im(:,:,1)));
    tc.verifySize(cxy, [2 K]);
    tc.verifyClass(resid, 'double');

    % assign to clusters
    cls = colorkmeans(im, cxy);
    tc.verifyEqual(size(cls), size(im(:,:,1)));
    tc.verifyEqual( min(cls(:)), 1);
    tc.verifyEqual( max(cls(:)), K);
end

function iblobs_test(tc)
    im = zeros(20, 20);
    im(4:16,1:6) = 1;
    im(8:12, 12:18) = 1;

    A2 = (16-4+1)*(6-1+1);
    A3 = (12-8+1)*(18-12+1);
    A1 = 20*20-A2-A3;

    f = iblobs(im);

    %tc.verifyClass(f, 'RegionFeature');
    tc.verifySize(f, [1 3]);

    s = char(f);
    tc.verifyEqual(size(s,1), 3);

    tc.verifyEqual(f(1).area, A1);
    tc.verifyEqual(f(1).umin, 1);
    tc.verifyEqual(f(1).vmin, 1);
    tc.verifyEqual(f(1).umax, 20);
    tc.verifyEqual(f(1).vmax, 20);
    tc.verifyEqual(f(1).class, 0);
    tc.verifyEqual(f(1).label, 1);
    tc.verifyEqual(f(1).touch, true);
    tc.verifyEqual(f(1).parent, uint32(0));

    % im(4:16,1:6) = 1;
    tc.verifyEqual(f(2).area, A2);
    tc.verifyEqual(f(2).umin, 1);
    tc.verifyEqual(f(2).vmin, 4);
    tc.verifyEqual(f(2).umax, 6);
    tc.verifyEqual(f(2).vmax, 16);
    tc.verifyEqual(f(2).class, 1);
    tc.verifyEqual(f(2).touch, true);
    tc.verifyEqual(f(2).parent, uint32(0));
    tc.assertEqual(f(2).uc, 3.5, 'absTol', 1e-3);
    tc.assertEqual(f(2).vc, 10, 'absTol', 1e-3);
    tc.assertEqual(f(2).aspect, 0.456, 'absTol', 1e-3);
    tc.assertEqual(f(2).theta, pi/2, 'absTol', 1e-3);

    % im(8:12, 12:18) = 1;
    tc.verifyEqual(f(3).area, A3);
    tc.verifyEqual(f(3).umin, 12);
    tc.verifyEqual(f(3).vmin, 8);
    tc.verifyEqual(f(3).umax, 18);
    tc.verifyEqual(f(3).vmax, 12);
    tc.verifyEqual(f(3).class, 1);
    tc.verifyEqual(f(3).label, 3);
    tc.verifyEqual(f(3).touch, false);
    tc.verifyEqual(f(3).parent, uint32(1));
    tc.assertEqual(f(3).uc, 15, 'absTol', 1e-3);
    tc.assertEqual(f(3).vc, 10, 'absTol', 1e-3);

    tc.assertEqual(f(3).aspect, 1/sqrt(2), 'absTol', 1e-3);
    tc.assertEqual(f(3).theta, 0, 'absTol', 1e-3);

    f.plot_box();

    % test filters
    f = iblobs(im, 'class', 1);
    tc.verifyEqual(length(f), 2);
    tc.verifyEqual(f.class, [1 1]);
    f = iblobs(im, 'class', 0);
    tc.verifyEqual(f.class, [0]);

    f = iblobs(im, 'touch', 0);
    tc.verifyEqual(f.touch, [false]);

    f = iblobs(im, 'touch', 1)
    tc.verifyEqual(f.touch, [true true]);

    f = iblobs(im, 'area', [1 Inf]);
    tc.verifyEqual(length(f), 3);

    f = iblobs(im, 'area', [30 40]);
    tc.verifyEqual(length(f), 1);
    tc.verifyEqual(f.area, [A3]);

    f = iblobs(im, 'area', [50 100]);
    tc.verifyEqual(length(f), 1);
    tc.verifyEqual(f.area, [A2]);

    f = iblobs(im, 'area', [20 100], 'class', 1);
    tc.verifyEqual(length(f), 2);
    tc.verifyTrue(all(f.area >= 20) && all(f.area <=100));
    tc.verifyEqual(f.class, [1 1]);

    f = iblobs(im, 'area', [20 100], 'class', 1, 'touch', 1);
    tc.verifyEqual(length(f), 1);
    tc.verifyTrue(f.area >= 20 && f.area <=100);
    tc.verifyEqual(f.class, 1);
    tc.verifyEqual(f.touch, true);

    % test boundary stuff
    f = iblobs(im, 'boundary');
    s = char(f);
    tc.verifyEqual(size(s,1), 3);

    tc.assertEqual(f(2).perimeter, 34.0, 'absTol', 1e-4);
    tc.assertEqual(f(3).perimeter, 20.0, 'absTol', 1e-4);

    % test circularity on a big circle for accuracy
    im = kcircle(100); n = size(im,1);
    im = [ zeros(10,n+20); zeros(n,10) im zeros(n,10); zeros(10,n+20)];
    f = iblobs(im, 'boundary');
    tc.assertEqual(length(f), 2);
    tc.assertEqual(f(2).circularity, 1, 'absTol', 5e-3);

    f = iblobs(im, 'boundary');
    tc.assertEqual(length(f), 2);
    tc.assertEqual(f(2).circularity, 1, 'absTol', 5e-3);
    tc.assertEqual(f(2).aspect, 1, 1e-3);
    tc.assertEqual(f.class, [0 1]);

    % test a circularity filter, now that we have a circle
    f = iblobs(im, 'boundary', 'circularity', [0.9 2]);
    tc.assertEqual(length(f), 1);
    tc.assertEqual(f.circularity, 1, 'absTol', 5e-3);
    tc.assertEqual(f.aspect, 1, 'absTol', 5e-3);
    tc.assertEqual(f.class, 1);

end

function MSER_test(tc)
    tc.assumeTrue(exist('vl_mser') > 0)
    castle = iread('castle2.png', 'double', 'grey');
    [mser,nsets] = imser(castle, 'area', [100 20000]);
    tc.verifyEqual(nsets, 71);
    tc.verifyEqual(size(castle), size(mser));
    tc.verifyEqual( min(mser(:)), 0);
    tc.verifyEqual( max(mser(:)), nsets-1);
end

function graphseg_test(tc)
    tc.assumeTrue(exist('graphseg') > 0)
    im = iread('58060.jpg');
    [label, m] = igraphseg(im, 1500, 100, 0.5);
    tc.verifyEqual(m, 28);
    tc.verifyEqual(size(im(:,:,1)), size(label));
    tc.verifyEqual( min(label(:)), 1);
    tc.verifyEqual( max(label(:)), m);
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function Volume=polygon2voxel(FV,VolumeSize,mode,Yxz,hollow)
% This function POLYGON2VOXEL will convert a Triangulated Mesh into a
% Voxel Volume which will contain the discretized mesh. Discretization of a
% polygon is done by splitting/refining the face, until the longest edge
% is smaller than 0.5 voxels. Followed by setting the voxel beneath the vertice
% coordinates of that small triangle to one.
%
% Volume=polygon2voxel(FV,VolumeSize,Mode,Yxz);
%
% Inputs,
%   FV : A struct containing FV.faces with a facelist Nx3 and FV.vertices
%        with a Nx3 vertice list. Such a structure is created by Matlab
%        Patch function
%   VolumeSize : The size of the output volume, example [100 100 100]
%   Mode : (optional) if set to:
%               'none', The vertices data is directly used as coordinates
%                       in the voxel volume.
%               'wrap', The vertices data is directly used as coordinates
%                       in the voxel volume, coordinates outside are
%                       circular wrapped to the inside.
%               'auto', The vertices data is translated and
%                       scaled with a scalar to fit inside the new volume.
%               'center', coordinate 0,0,0 is set as the center of the volume
%                       instead of the corner of the voxel volume.
%               'clamp', The vertices data is directly used as coordinates
%                       in the voxel volume, coordinates outside are
%                       clamped to the inside.
%   (Optional)
%   Yxz : If true (default) use Matlab convention 1th dimension Y,
%         2th dimension X, and last dimension Z. Otherwise 1th
%         dimension is X, 2th Y and last Z.
%
%
% Outputs,
%   Volume : The 3D logical volume, with all voxels part of the discretized
%           mesh one, and all other voxels zero.
%
% Example,
%   % Compile the c-coded function
%   mex polygon2voxel_double.c -v
%
%   % Load a triangulated mesh of a sphere
%   load sphere;
%
%   % Show the mesh
%   figure, patch(FV,'FaceColor',[1 0 0]); axis square;
%
%   % Convert the mesh to a voxelvolume
%   Volume=polygon2voxel(FV,[50 50 50],'auto');
%
%   % Show x,y,z slices
%   figure,
%   subplot(1,3,1), imshow(squeeze(Volume(25,:,:)));
%   subplot(1,3,2), imshow(squeeze(Volume(:,25,:)));
%   subplot(1,3,3), imshow(squeeze(Volume(:,:,25)));
%
%   %  Show iso surface of result
%   figure, patch(isosurface(Volume,0.1), 'Facecolor', [1 0 0]);
%
% Example2,
%   % Compile the c-coded function
%   mex polygon2voxel_double.c -v
%
%   % Make A Volume with a few blocks
%   I = false(120,120,120);
%   I(40:60,50:70,60:80)=1; I(60:90,45:75,60:90)=1;
%   I(20:60,40:80,20:60)=1; I(60:110,35:85,10:60)=1;
%
%   % Convert the volume to a triangulated mesh
%   FV = isosurface(I,0.8);
%
%   % Convert the triangulated mesh back to a surface in a volume
%   J = polygon2voxel(FV,[120, 120, 120],'none');
%   % Fill the volume
%   J=imfill(J,'holes');
%
%   % Differences between original and reconstructed
%   VD = abs(J-I);
%
%   % Show the original Mesh and Mesh of new volume
%   figure,
%   subplot(1,3,1),  title('original')
%     patch(FV,'facecolor',[1 0 0],'edgecolor','none'), camlight;view(3);
%   subplot(1,3,2), title('converted');
%     patch(isosurface(J,0.8),'facecolor',[0 0 1],'edgecolor','none'), camlight;view(3);
%   subplot(1,3,3), title('difference');
%     patch(isosurface(VD,0.8),'facecolor',[0 0 1],'edgecolor','none'), camlight;view(3);
%
% Function is written by D.Kroon University of Twente (May 2009)
%
% Simple modification by Jianxiong to fill in the holes

if ~exist('hollow','var')
    hollow = false;
end

if(nargin<4), Yxz=true; end

% Check VolumeSize size
if(length(VolumeSize)==1)
    VolumeSize=[VolumeSize VolumeSize VolumeSize];
end
if(length(VolumeSize)~=3)
    error('polygon2voxel:inputs','VolumeSize must be a array of 3 elements ')
end

% Volume Size must always be an integer value
VolumeSize=round(VolumeSize);

sizev=size(FV.vertices);
% Check size of vertice array
if((sizev(2)~=3)||(length(sizev)~=2))
    error('polygon2voxel:inputs','The vertice list is not a m x 3 array')
end

sizef=size(FV.faces);
% Check size of vertice array
if((sizef(2)~=3)||(length(sizef)~=2))
    error('polygon2voxel:inputs','The vertice list is not a m x 3 array')
end

% Check if vertice indices exist
if(max(FV.faces(:))>size(FV.vertices,1))
    error('polygon2voxel:inputs','The face list contains an undefined vertex index')
end

% Check if vertice indices exist
if(min(FV.faces(:))<1)
    error('polygon2voxel:inputs','The face list contains an vertex index smaller then 1')
end

% Matlab dimension convention YXZ
if(Yxz)
    FV.vertices=FV.vertices(:,[2 1 3]);
end

switch(lower(mode(1:2)))
    case {'au'} % auto
        % Make all vertices-coordinates positive
        FV.vertices=FV.vertices-min(FV.vertices(:));
        scaling=min((VolumeSize-1)./(max(FV.vertices(:))));
        % Make the vertices-coordinates to range from 0 to 100
        FV.vertices=FV.vertices*scaling+1;
        Wrap=0;
    case {'ce'} % center
        % Center the vertices
        FV.vertices=FV.vertices+repmat((VolumeSize/2),size(FV.vertices,1),1);
        Wrap=0;
    case {'wr'} %wrap
        Wrap=1;
    case{'cl'} % clamp
        Wrap=2;
    otherwise
        Wrap=0;
end

% Separate the columns;
FacesA=double(FV.faces(:,1));
FacesB=double(FV.faces(:,2));
FacesC=double(FV.faces(:,3));
VerticesX=double(FV.vertices(:,1));
VerticesY=double(FV.vertices(:,2));
VerticesZ=double(FV.vertices(:,3));

% Volume size to double
VolumeSize=double(VolumeSize);

% Call the mex function
Volume=polygon2voxel_double(FacesA,FacesB,FacesC,VerticesX,VerticesY,VerticesZ,VolumeSize,Wrap);

if ~hollow
    for i=1:size(Volume,3)
        Volume(:,:,i) = imfill(Volume(:,:,i),'holes');
    end
end

author: zouchuhang

repo: LayoutNet


text:
%  When you load any ANALYZE or NIfTI file with 'load_nii.m', and view
%  it with 'view_nii.m', you may find that the image is L-R flipped.
%  This is because of the confusion of radiological and neurological
%  convention in the medical image before NIfTI format is adopted. You
%  can find more details from:
%
%  http://www.rotman-baycrest.on.ca/~jimmy/UseANALYZE.htm
%
%  Sometime, people even want to convert RAS (standard orientation) back
%  to LAS orientation to satisfy the legend programs or processes. This
%  program is only written for those purpose. So PLEASE BE VERY CAUTIOUS
%  WHEN USING THIS 'FLIP_LR.M' PROGRAM.
%
%  With 'flip_lr.m', you can convert any ANALYZE or NIfTI (no matter
%  3D or 4D) file to a flipped NIfTI file. This is implemented simply
%  by flipping the affine matrix in the NIfTI header. Since the L-R
%  orientation is determined there, so the image will be flipped.
%
%  Usage: flip_lr(original_fn, flipped_fn, [old_RGB],[tolerance],[preferredForm])
%
%  original_fn  -  filename of the original ANALYZE or NIfTI (3D or 4D) file
%
%  flipped_fn  -  filename of the L-R flipped NIfTI file
%
%  old_RGB (optional)  -  a scale number to tell difference of new RGB24
%	from old RGB24. New RGB24 uses RGB triple sequentially for each
%	voxel, like [R1 G1 B1 R2 G2 B2 ...]. Analyze 6.0 from AnalyzeDirect
%	uses old RGB24, in a way like [R1 R2 ... G1 G2 ... B1 B2 ...] for
%	each slices. If the image that you view is garbled, try to set
%	old_RGB variable to 1 and try again, because it could be in
%	old RGB24. It will be set to 0, if it is default or empty.
%
%  tolerance (optional) - distortion allowed for non-orthogonal rotation
%	or shearing in NIfTI affine matrix. It will be set to 0.1 (10%),
%	if it is default or empty.
%
%  preferredForm (optional)  -  selects which transformation from voxels
%	to RAS coordinates; values are s,q,S,Q.  Lower case s,q indicate
%	"prefer sform or qform, but use others if preferred not present".
%	Upper case indicate the program is forced to use the specificied
%	tranform or fail loading.  'preferredForm' will be 's', if it is
%	default or empty.	- Jeff Gunter
%
%  Example: flip_lr('avg152T1_LR_nifti.nii', 'flipped_lr.nii');
%           flip_lr('avg152T1_RL_nifti.nii', 'flipped_rl.nii');
%
%  You will find that 'avg152T1_LR_nifti.nii' and 'avg152T1_RL_nifti.nii'
%  are the same, and 'flipped_lr.nii' and 'flipped_rl.nii' are also the
%  the same, but they are L-R flipped from 'avg152T1_*'.
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function flip_lr(original_fn, flipped_fn, old_RGB, tolerance, preferredForm)

   if ~exist('original_fn','var') | ~exist('flipped_fn','var')
      error('Usage: flip_lr(original_fn, flipped_fn, [old_RGB],[tolerance])');
   end

   if ~exist('old_RGB','var') | isempty(old_RGB)
      old_RGB = 0;
   end

   if ~exist('tolerance','var') | isempty(tolerance)
      tolerance = 0.1;
   end

   if ~exist('preferredForm','var') | isempty(preferredForm)
      preferredForm= 's';				% Jeff
   end

   nii = load_nii(original_fn, [], [], [], [], old_RGB, tolerance, preferredForm);
   M = diag(nii.hdr.dime.pixdim(2:5));
   M(1:3,4) = -M(1:3,1:3)*(nii.hdr.hist.originator(1:3)-1)';
   M(1,:) = -1*M(1,:);
   nii.hdr.hist.sform_code = 1;
   nii.hdr.hist.srow_x = M(1,:);
   nii.hdr.hist.srow_y = M(2,:);
   nii.hdr.hist.srow_z = M(3,:);
   save_nii(nii, flipped_fn);

   return;					% flip_lr

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function ft_preamble(cmd, varargin)

% FT_PREAMBLE is a helper function that is included in many of the FieldTrip
% functions and which takes care of some general settings and operations at the
% begin of the function.
%
% This ft_preamble m-file is a function, but internally it executes a
% number of private scripts in the callers workspace. This allows the
% private script to access the variables in the callers workspace and
% behave as if the script were included as a header file in C-code.
%
% See also FT_POSTAMBLE

% Copyright (C) 2011-2012, Robert Oostenveld, DCCN
%
% 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$

% ideally this would be a script, because the local variables would then be
% shared with the calling function. Instead, this is a function which then
% passes the variables explicitely to another script which is eval'ed.

% the following section ensures that these scripts are included as
% dependencies when using the MATLAB compiler
%
%#function ft_preamble_init
%#function ft_preamble_debug
%#function ft_preamble_trackconfig
%#function ft_preamble_provenance
%#function ft_preamble_loadvar
%#function ft_preamble_randomseed

% this is a trick to pass the input arguments into the ft_preamble_xxx script
assignin('caller', 'iW1aenge_preamble', varargin);

full_cmd = ['ft_preamble_' cmd];
cmd_exists = false;

if exist(full_cmd, 'file')
  % Matlab can find commands in a private subdirectory; Octave cannot.
  % If pwd is already the private directory, or if using Matlab then
  % the command can be evaluated directly
  cmd_exists = true;

elseif ~ft_platform_supports('exists-in-private-directory')
  % Octave does not find files by name in a private directory, so the full
  % filename must be specified.
  private_dir=fullfile(fileparts(which(mfilename)),'private');
  full_path=fullfile(private_dir,[full_cmd '.m']);

  cmd_exists=exist(full_path,'file');
  full_cmd_parts = {'ft_tmp_orig_pwd=pwd();',...
    'ft_tmp_orig_pwd_cleaner='...
    'onCleanup(@()cd(ft_tmp_orig_pwd));',...
    sprintf('cd(''%s'');',private_dir),...
    [full_cmd ';'],...
    'clear ft_tmp_orig_pwd_cleaner;'};
  full_cmd = sprintf('%s',full_cmd_parts{:});
end

if ~cmd_exists
  % XXX earlier versions would not do anything if ~cmd_exists,
  % but fail silently (without raising an error or warning).
  % Should that behavior be kept?
  ft_error('Could not run %s - does not seem to exist', full_cmd);
end

evalin('caller', full_cmd);

author: spm


text:
function [feature] = TDD(inf,tra,cnn_feature,scale_x,scale_y,num_cell)
% TDD: perform trajectory pooling over convolutional feature maps.
% Input:
%       inf: information of trajectories from iDTs (10*N)
%       traj: extracted trajectories (2L*N)
%       cnn_feature: cnn feature maps (convlutional layers: W*H*C*L)
%       scale_x: width ratio
%       scale_y: height ratio
%       num_cell: the number of cell in temporal dimension
% Output:
%       feature: trajectory pooled descriptors ((C*NUM_CELL) *N)

if ~isempty(inf)
	ind = inf(7,:)==1;
	inf = inf(:,ind);
	tra = tra(:,ind);
end

if ~isempty(inf)
	NUM_DIM = size(cnn_feature,3);
	NUM_DES = size(inf,2);
	TRA_LEN = size(tra,1)/2;

	num_fea = TRA_LEN / num_cell;

	pos = reshape(tra,2,[])-1;
	pos = round(bsxfun(@rdivide,pos,[scale_x;scale_y]) + 1);
	pos = bsxfun(@max,pos,[1;1]);
	pos = bsxfun(@min,pos,[size(cnn_feature,2);size(cnn_feature,1)]);
	pos = reshape(pos,TRA_LEN*2,[]);

	cnn_feature = permute(cnn_feature,[1,2,4,3]);
	offset = [TRA_LEN-1:-1:0];
	size_mat = [size(cnn_feature,1),size(cnn_feature,2),size(cnn_feature,3)];
	cnn_feature = reshape(cnn_feature,[],NUM_DIM);

	cur_x = pos(1:2:end,:);
	cur_y = pos(2:2:end,:);
	cur_t = bsxfun(@minus,inf(1,:),offset');

	tmp = cnn_feature(sub2ind(size_mat,cur_y,cur_x,cur_t),:)';
	tmp = reshape(tmp,NUM_DIM,num_fea,[]);
	feature = reshape(sum(tmp,2),[],NUM_DES);
else
	feature = [];
end

end

author: wanglimin

repo: TDD

save_path: github-repos/MATLAB/wanglimin-TDD


text:
function [varargout]=patchCleanUnused(F,V)

% function [Fc,Vc,indFix2,logicValid,indUni]=patchCleanUnused(F,V)
%-------------------------------------------------------------------------
%
%
% 2010 Created
% 2019/03/04 Added basic description header to function
% 2019/04/22 Added handling of cell input
%-------------------------------------------------------------------------

%%

if isa(F,'cell')
    indUni=[];
    for q=1:1:numel(F)
        f=F{q};
        indUni=[indUni;f(:)];
    end
else
    indUni=unique(F(:));
end
indUni=unique(indUni(~isnan(indUni)));

%%

numPointsOriginal=size(V,1); %Number of original points
numPointsNew=numel(indUni); %Number of points in new set
Vc=V(indUni,:); %Select relevant points

%Fix indices in faces matrix
indFix1=1:1:numPointsNew;
indFix2=zeros(numPointsOriginal,1);
indFix2(indUni)=indFix1;

Fc=F;
if isa(F,'cell')
    for q=1:1:numel(F)
        fc=F{q};
        logicValid=~isnan(fc);
        fc(logicValid)=indFix2(fc(logicValid));
        Fc{q}=fc;
    end
else
    logicValid=~isnan(F);
    Fc(logicValid)=indFix2(F(logicValid));
end

%Output
varargout{1}=Fc;
varargout{2}=Vc;
varargout{3}=indFix2;
varargout{4}=logicValid;
varargout{5}=indUni;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function alpha = p15_alpha ( )

%*****************************************************************************80
%
%% P15_ALPHA returns ALPHA for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real ALPHA, the value of ALPHA.
%
  alpha = 0.0;

  return
end


text:
function p03_demo ( iteration_max, h, fh )

%*****************************************************************************80
%
%% P03_DEMO runs the 2D demo problem #3, with mesh size H.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    06 February 2006
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer ITERATION_MAX, the maximum number of iterations that DISTMESH
%    should take.  (The program might take fewer iterations if it detects convergence.)
%
%    Input, real H, the mesh spacing parameter.
%
%    Input, external FH, the mesh density function.
%
  if ( nargin < 1 )
    iteration_max = 300;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P03_DEMO - Note:\n' );
    fprintf ( 1, '  No value of ITERATION_MAX was supplied.\n' );
    fprintf ( 1, '  The default value ITERATION_MAX = %d will be used.\n', ...
      iteration_max );
  end

  if ( nargin < 2 )
    h = 0.15;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P03_DEMO - Note:\n' );
    fprintf ( 1, '  No value of H was supplied.\n' );
    fprintf ( 1, '  The default value H = %f will be used.\n', h );
  end

  if ( nargin < 3 )
    fh = @p03_fh;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P03_DEMO - Note:\n' );
    fprintf ( 1, '  No value of FH was supplied.\n' );
    fprintf ( 1, '  The default uniform density function will be used.\n', h );
  end
%
%  Put the random number generator into a fixed initial state.
%
  rand ( 'state', 111 );
%
%  Set the rendering method for the current figure to Z-buffering.
%
  set ( gcf, 'rend', 'z' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 3:\n' );
  fprintf ( 1, '  The square with a circular hole, h = %f\n', h )

  fd = @p03_fd;
  box = [ -1.0,-1.0; 1.0,1.0 ];
  fixed = [ -1.0,-1.0; -1.0,1.0; 1.0,-1.0; 1.0,1.0 ];

  [ p, t ] = distmesh_2d ( fd, fh, h, box, iteration_max, fixed );

  post_2d ( p, t, fh )
%
%  Write a PostScript image of the triangulation.
%
  [ node_num, junk ] = size ( p );
  [ tri_num, junk ] = size ( t );
  p = p';
  t = t';
  node_show = 0;
  triangle_show = 1;

  triangulation_order3_plot ( 'p03_mesh.eps', node_num, p, tri_num, ...
    t, node_show, triangle_show );
%
%  Write a text file containing the nodes.
%
  r8mat_write ( 'p03_nodes.txt', 2, node_num, p );
%
%  Write a text file containing the triangles.
%
  i4mat_write ( 'p03_elements.txt', 3, tri_num, t );

  return
end


text:
function [h,g,a,info] = wfilt_optsymb(N)
%WFILT_OPTSYMB  Optimizatized Symmetric Self-Hilbertian Filters
%
%   Usage: [h,g,a] = wfilt_optsymb(N);
%
%   `[h,g,a]=wfiltdt_optsymb(N)` with $N \in {1,2,3}$ returns filters
%   suitable with optimized symmetry suitable for for dual-tree complex
%   wavelet transform tree B.
%
%   Examples:
%   ---------
%   :::
%     wfiltinfo('optsymb3');
%
%   References: dubase08
%

% AUTHOR: Zdenek Prusa

[ha,~,a,info] = wfilt_optsyma(N);

hlp = ha{1}.h;
offset = -(numel(hlp)/2);
range = (0:numel(hlp)-1) + offset;

% Create the filters according to the reference paper.
%
% REMARK: The phase of the alternating +1 and -1 is crucial here.
%
    harr = [...
            flipud(hlp),...
            (-1).^(range).'.*hlp,...
            ];

htmp=mat2cell(harr,size(harr,1),ones(1,size(harr,2)));

h(1:2,1) = cellfun(@(hEl)struct('h',hEl,'offset',offset),htmp(1:2),...
                   'UniformOutput',0);
g = h;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function result = AngleAxisRotatePoint(angle_axis, pt)

angle_axis = angle_axis(1:3);

theta2 = dot(angle_axis,angle_axis);

if (theta2 > 0.0)
    % Away from zero, use the rodriguez formula
    %
    %   result = pt costheta + (w x pt) * sintheta + w (w . pt) (1 - costheta)
    %
    % We want to be careful to only evaluate the square root if the
    % norm of the angle_axis vector is greater than zero. Otherwise
    % we get a division by zero.

    theta = sqrt(theta2);
    w = angle_axis / theta;

    costheta = cos(theta);
    sintheta = sin(theta);

    % w_cross_pt = cross(w, pt);
    w_cross_pt = xprodmat(w) * pt;

    %w_dot_pt = dot(w, pt);
    w_dot_pt = w * pt;

    result = pt * costheta + w_cross_pt * sintheta + (w' * (1 - costheta)) * w_dot_pt;

else
    % Near zero, the first order Taylor approximation of the rotation
    % matrix R corresponding to a vector w and angle w is
    %
    %   R = I + hat(w) * sin(theta)
    %
    % But sintheta ~ theta and theta * w = angle_axis, which gives us
    %
    %  R = I + hat(w)
    %
    % and actually performing multiplication with the point pt, gives us
    % R * pt = pt + w x pt.
    %
    % Switching to the Taylor expansion at zero helps avoid all sorts
    % of numerical nastiness.

    %w_cross_pt = cross(angle_axis, pt);
    w_cross_pt = xprodmat(angle_axis) * pt; % vectorize version

    result = pt + w_cross_pt;

end

author: jianxiongxiao

repo: ProfXkit


text:
function [data,units] = compute_dangle_smallest_wing(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  danglel = diff(-trx(fly).wing_anglel);
  dangler = diff(trx(fly).wing_angler);
  data{i} = danglel;
  idx = trx(fly).wing_arear_mm(1:end-1) <= trx(fly).wing_areal_mm(1:end-1);
  data{i}(idx) = dangler(idx);
  data{i} = data{i} ./ trx(fly).dt;

end
units = parseunits('rad/s');

author: kristinbranson

repo: JAABA


text:
%klibrary ' '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros library.pane file
%
% Parameters:
% MultiChoice: lang 'Language Dependency?', default: 1: 'what language is this library written in?'
%    Choices are:
%   1: 'C'
%   2: 'C++'
%   3: 'Fortran'
%
% Example: klibrary( {'lang',1})
%
% Khoros helpfile follows below:

function varargout = klibrary(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = klibrary(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'lang', 1};
maxval={0};
minval={0};
istoggle=[0];
was_set=istoggle * 0;
paramtype={'MultiChoice'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; 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(0);
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
callKhoros([w 'library" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function kern = gaussianwhiteKernExpandParam(kern, params)

% GAUSSIANWHITEKERNEXPANDPARAM Create kernel structure from gaussian white
%                              kernel's parameters.
% FORMAT
% DESC returns a gaussian white kernel structure filled with the parameters in the given
%	vector. This is used as a helper function to enable parameters to be
%	optimised in, for example, the NETLAB optimisation functions.
% RETURN kern : kernel structure with the given parameters in the relevant
%	   locations.
% ARG kern : the kernel structure in which the parameters are to be
% ARG param : vector of parameters which are to be placed in the kernel
%	   structure.
%
% SEEALSO : gaussianwhiteKernParamInit, gaussianwhiteKernExtractParam, kernExpandParam
%
% COPYRIGHT : Mauricio Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009

% KERN

if kern.isArd
    if kern.nIndFunct == 1,
        kern.precisionT =  params(1:end-1)';
    else
        kern.precisionT = reshape(params(1:end-1), size(kern.precisionT));
    end
else
    kern.precisionT =  params(1:end-1);
end
kern.sigma2Noise = params(end);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function y = downsample2(x,K)

for i = 1:size(x,3)
y(:,:,i) = downsample(downsample(x(:,:,i),K)',K)';
end

end

author: cszn

repo: IRCNN

save_path: github-repos/MATLAB/cszn-IRCNN


text:
function mv = mv_helpMenu(mv, hfig);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (5): help menu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mv.ui.helpMenu = uimenu(hfig, 'ForegroundColor', [.6 .2 0], 'Label', 'Help', ...
                              'Separator', 'on');

% web page callback:
% web http://white.stanford.edu/newlm/index.php/Multi_Voxel_UI
cb = 'web http://white.stanford.edu/newlm/index.php/Multi_Voxel_UI';
mv.ui.webHelp = uimenu(mv.ui.helpMenu, 'Label', 'Time Course UI page', 'Separator', 'off', ...
    'CallBack', cb);

% web page (external browser) callback:
% web http://white.stanford.edu/newlm/index.php/Multi_Voxel_UI -browser
cb = 'web http://white.stanford.edu/newlm/index.php/Multi_Voxel_UI -browser';
cb = 'web http://white/newlm/index.php/Main_Page -browser';
mv.ui.webHelp2 = uimenu(mv.ui.helpMenu, 'Label', 'Time Course UI page (external browser)', 'Separator', 'off', ...
    'CallBack', cb);

% mrVista web page callback:
% web web http://white.stanford.edu/newlm/index.php/Main_Page
cb = 'web http://white.stanford.edu/newlm/index.php/Main_Page';
mv.ui.webWiki = uimenu(mv.ui.helpMenu,  'Label',  'mrVista wiki',  ...
    'Separator',  'off',  'Callback',  cb);

% mrVista web page (external browser) callback:
% web web http://white.stanford.edu/newlm/index.php/Main_Page
% -browser
mv.ui.webWiki2 = uimenu(mv.ui.helpMenu,  'Label',  'mrVista wiki (external browser)',  ...
    'Separator',  'off',  'Callback',  cb);

% identify callback for menu item
cb = 'helpFindCallback; ';
uimenu(mv.ui.helpMenu,  'Label',  'Identify Callback for a menu item',  ...
        'Separator',  'on',  'Callback',  cb);

return

author: vistalab

repo: vistasoft


text:
function [dateNum, dateStr] = datestr2datenum(dateStr)
dateNum = 0;

c = str2cell(dateStr, ' ');
if length(c)<2
    return
end
dateStr = c{1};
timeStr = c{2};

dateStr(dateStr=='-')='';
timeStr(timeStr==':')='';

timeNum = str2num(timeStr);
dateNum = str2num(dateStr) + timeNum/1e6;

dateStr = [c{1}, ' ', c{2}];

author: BUNPC


text:
function plotbyfactor(x,y,f)
%
% very simplistic version of plotbyfactor.
% no symbols, will crash if more than 8 levels.
%

z = unique(f);
cols = ['b','r','g','m','c','y','k','w'];

for (i = 1:length(z))
  u = find(f==z(i));
  plot(x(u),y(u),'.','color',cols(i));
  hold on;
end;
hold off;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function y = perform_waveatoms_transform(x,dir, options)

% perform_waveatoms_transform - interface to WaveAtom transform
%
%   y = perform_waveatoms_transform(x,dir, options);
%
%   The waveatom toolbox can be downloaded from
%       http://www.waveatom.org/
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
issym = getoptions(options, 'issym', 0);
force_real = getoptions(options, 'force_real', 1);

if dir==1
    d = nb_dims(x);
else
    if iscell(x)
        d = nb_dims(x{1});
    else
        d = nb_dims(x);
    end
end

pat = 'p'; tp = 'ortho';
if d==1
    if dir==1
        if issym
            y = fatom1sym(x,pat,[1 1]);
        else
            y = fwa1(x,pat,tp);
        end
    else
        if issym
            y = iatom1sym(x,pat,[1 1]);
        else
            y = iwa1(x,pat,tp);
        end
    end
else
    if dir==1
        if issym
            y = fwa2sym(x,pat,[1 1]);
        else
            y = fwa2(x,pat,tp);
        end
    else
        if issym
            y = iwa2sym(x,pat,[1 1]);
        else
            y = iwa2(x,pat,tp);
        end
    end
end

if force_real && dir==-1
    y = real(y);
end

function c = fwa2(x,pat,tp)
% fwa2 - 2D forward wave atom transform
% -----------------
% INPUT
% --
% x is a real N-by-N matrix. N is a power of 2.
% --
% pat specifies the type of frequency partition which satsifies
% parabolic scaling relationship. pat can either be 'p' or 'q'.
% --
% tp is the type of tranform.
% 	'ortho': orthobasis
% 	'directional': real-valued frame with single oscillation direction
% 	'complex': complex-valued frame
% -----------------
% OUTPUT
% --
% c is a cell array which contains the wave atom coefficients. If
% tp=='ortho', then c{j}{m1,m2}(n1,n2) is the coefficient at scale j,
% frequency index (m1,m2) and spatial index (n1,n2). If
% tp=='directional', then c{j,d}{m1,m2}(n1,n2) with d=1,2 are the
% coefficients at scale j, frequency index (m1,m2) and spatial index
% (n1,n2). If tp=='complex', then c{j,d}{m1,m2)(n1,n2) with d=1,2,3,4
% are the coefficients at scale j, frequency index (m1,m2) and spatial
% index (n1,n2).
% -----------------
% Written by Lexing Ying and Laurent Demanet, 2007

if( ismember(tp, {'ortho','directional','complex'})==0 | ismember(pat, {'p','q','u'})==0 )    error('wrong');  end

if(strcmp(tp, 'ortho')==1)
    %---------------------------------------------------------
    N = size(x,1);
    H = N/2;
    lst = freq_pat(H,pat);
    %------------------
    f = fft2(x) / sqrt(prod(size(x)));
    A = N;
    c = cell(length(lst),1);
    %------------------
    for s=1:length(lst)
        nw = length(lst{s});
        c{s} = cell(nw,nw);
        for I=0:nw-1
            for J=0:nw-1
                if(lst{s}(I+1)==0 & lst{s}(J+1)==0)
                    c{s}{I+1,J+1} = [];
                else
                    B = 2^(s-1);
                    D = 2*B;
                    Ict = I*B;      Jct = J*B; %starting position in freq
                    if(mod(I,2)==0)
                        Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
                    else
                        Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
                    end
                    if(mod(J,2)==0)
                        Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
                    else
                        Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
                    end
                    res = zeros(D,D);
                    for id=0:1
                        if(id==0)
                            Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                        else
                            Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                        end
                        for jd=0:1
                            if(jd==0)
                                Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                            else
                                Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                            end
                            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                        end
                    end
                    c{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
                end
            end
        end
    end

elseif(strcmp(tp, 'directional')==1)
    %---------------------------------------------------------
    N = size(x,1);
    H = N/2;
    lst = freq_pat(H,pat);
    %------------------
    f = fft2(x) / sqrt(prod(size(x)));
    A = N;
    c1 = cell(length(lst),1);
    c2 = cell(length(lst),1);
    %------------------
    for s=1:length(lst)
        nw = length(lst{s});
        c1{s} = cell(nw,nw);
        c2{s} = cell(nw,nw);
        for I=0:nw-1
            for J=0:nw-1
                if(lst{s}(I+1)==0 & lst{s}(J+1)==0)
                    c1{s}{I+1,J+1} = [];
                    c2{s}{I+1,J+1} = [];
                else
                    B = 2^(s-1);
                    D = 2*B;
                    res = zeros(D,D);
                    Ict = I*B;      Jct = J*B; %starting position in freq
                    if(mod(I,2)==0)
                        Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
                    else
                        Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
                    end
                    if(mod(J,2)==0)
                        Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
                    else
                        Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
                    end

                    res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    c1{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));

                    res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    c2{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
                end
            end
        end
    end
    c = [c1 c2];

elseif(strcmp(tp, 'complex')==1)
    %---------------------------------------------------------
    N = size(x,1);
    H = N/2;
    lst = freq_pat(H,pat);
    %------------------
    f = fft2(x) / sqrt(prod(size(x)));
    A = N;
    c1 = cell(length(lst),1);
    c2 = cell(length(lst),1);
    c3 = cell(length(lst),1);
    c4 = cell(length(lst),1);
    %------------------
    for s=1:length(lst)
        nw = length(lst{s});
        c1{s} = cell(nw,nw);
        c2{s} = cell(nw,nw);
        c3{s} = cell(nw,nw);
        c4{s} = cell(nw,nw);
        for I=0:nw-1
            for J=0:nw-1
                if(lst{s}(I+1)==0 & lst{s}(J+1)==0)
                    c1{s}{I+1,J+1} = [];
                    c2{s}{I+1,J+1} = [];
                    c3{s}{I+1,J+1} = [];
                    c4{s}{I+1,J+1} = [];
                else
                    B = 2^(s-1);
                    D = 2*B;
                    res = zeros(D,D);
                    Ict = I*B;      Jct = J*B; %starting position in freq
                    if(mod(I,2)==0)
                        Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
                    else
                        Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
                    end
                    if(mod(J,2)==0)
                        Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
                    else
                        Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
                    end

                    res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    c1{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));

                    res = zeros(D,D);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    c2{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));

                    res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    c3{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));

                    res = zeros(D,D);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
                    c4{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
                end
            end
        end
    end
    c = [c1 c2 c3 c4];
end

function x = iwa2(c,pat,tp)
% iwa2 - 2D inverse wave atom transform
% -----------------
% INPUT
% --
% c is a cell array which contains the wave atom coefficients. If
% tp=='ortho', then c{j}{m1,m2}(n1,n2) is the coefficient at scale j,
% frequency index (m1,m2) and spatial index (n1,n2). If
% tp=='directional', then c{j,d}{m1,m2}(n1,n2) with d=1,2 are the
% coefficients at scale j, frequency index (m1,m2) and spatial index
% (n1,n2). If tp=='complex', then c{j,d}{m1,m2)(n1,n2) with d=1,2,3,4
% are the coefficients at scale j, frequency index (m1,m2) and spatial
% index (n1,n2).
% --
% pat specifies the type of frequency partition which satsifies
% parabolic scaling relationship. pat can either be 'p' or 'q'.
% --
% tp is the type of tranform.
% 	'ortho': orthobasis
% 	'directional': real-valued frame with single oscillation direction
% 	'complex': complex-valued frame
% -----------------
% OUTPUT
% --
% x is a real N-by-N matrix. N is a power of 2.
% -----------------
% Written by Lexing Ying and Laurent Demanet, 2007

if( ismember(tp, {'ortho','directional','complex'})==0 | ismember(pat, {'p','q','u'})==0 )    error('wrong');  end

if(strcmp(tp, 'ortho')==1)
    %---------------------------------------------------------
    T = 0;
    for s=1:length(c)
        nw = length(c{s});
        for I=1:nw
            for J=1:nw
                T = T + prod(size(c{s}{I,J}));
            end
        end
    end
    N = sqrt(T);
    H = N/2;
    lst = freq_pat(H,pat);
    A = N;
    f = zeros(A,A);
    %------------------
    for s=1:length(lst)
        nw = length(lst{s});
        for I=0:nw-1
            for J=0:nw-1
                if(~isempty(c{s}{I+1,J+1}))
                    B = 2^(s-1);
                    D = 2*B;
                    Ict = I*B;      Jct = J*B; %starting position in freq
                    if(mod(I,2)==0)
                        Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
                    else
                        Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
                    end
                    if(mod(J,2)==0)
                        Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
                    else
                        Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
                    end
                    res = fft2(c{s}{I+1,J+1}) / sqrt(prod(size(c{s}{I+1,J+1}))); %res = zeros(D,D);
                    for id=0:1
                        if(id==0)
                            Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                        else
                            Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                        end
                        for jd=0:1
                            if(jd==0)
                                Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                            else
                                Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                            end
                            f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);
                        end
                    end
                end
            end
        end
    end
    %------------------
    x = ifft2(f) * sqrt(prod(size(f)));

elseif(strcmp(tp, 'directional')==1)
    %---------------------------------------------------------
    c1 = c(:,1);
    c2 = c(:,2);

    T = 0;
    for s=1:length(c1)
        nw = length(c1{s});
        for I=1:nw
            for J=1:nw
                T = T + prod(size(c1{s}{I,J}));
            end
        end
    end
    N = sqrt(T);
    H = N/2;
    lst = freq_pat(H,pat);
    A = N;
    f = zeros(A,A);
    %------------------
    for s=1:length(lst)
        nw = length(lst{s});
        for I=0:nw-1
            for J=0:nw-1
                if(~isempty(c1{s}{I+1,J+1}))
                    B = 2^(s-1);
                    D = 2*B;
                    Ict = I*B;      Jct = J*B; %starting position in freq
                    if(mod(I,2)==0)
                        Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
                    else
                        Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
                    end
                    if(mod(J,2)==0)
                        Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
                    else
                        Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
                    end

                    res = fft2(c1{s}{I+1,J+1}) / sqrt(prod(size(c1{s}{I+1,J+1}))); %res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);

                    res = fft2(c2{s}{I+1,J+1}) / sqrt(prod(size(c2{s}{I+1,J+1}))); %res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);
                end
            end
        end
    end
    x = ifft2(f) * sqrt(prod(size(f)));

elseif(strcmp(tp, 'complex')==1)
    %---------------------------------------------------------
    c1 = c(:,1);
    c2 = c(:,2);
    c3 = c(:,3);
    c4 = c(:,4);

    T = 0;
    for s=1:length(c1)
        nw = length(c1{s});
        for I=1:nw
            for J=1:nw
                T = T + prod(size(c1{s}{I,J}));
            end
        end
    end
    N = sqrt(T);
    H = N/2;
    lst = freq_pat(H,pat);
    A = N;
    f = zeros(A,A);
    %------------------
    for s=1:length(lst)
        nw = length(lst{s});
        for I=0:nw-1
            for J=0:nw-1
                if(~isempty(c1{s}{I+1,J+1}))
                    B = 2^(s-1);
                    D = 2*B;
                    Ict = I*B;      Jct = J*B; %starting position in freq
                    if(mod(I,2)==0)
                        Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
                    else
                        Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
                    end
                    if(mod(J,2)==0)
                        Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
                    else
                        Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
                    end

                    res = fft2(c1{s}{I+1,J+1}) / sqrt(prod(size(c1{s}{I+1,J+1}))); %res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);

                    res = fft2(c2{s}{I+1,J+1}) / sqrt(prod(size(c2{s}{I+1,J+1}))); %res = zeros(D,D);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);

                    res = fft2(c3{s}{I+1,J+1}) / sqrt(prod(size(c3{s}{I+1,J+1}))); %res = zeros(D,D);
                    Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
                    Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);

                    res = fft2(c4{s}{I+1,J+1}) / sqrt(prod(size(c4{s}{I+1,J+1}))); %res = zeros(D,D);
                    Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
                    Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                    f(mod(Idx,A)+1,mod(Jdx,A)+1) = f(mod(Idx,A)+1,mod(Jdx,A)+1) + ( Icf.'*Jcf ) .* res(mod(Idx,D)+1,mod(Jdx,D)+1);
                end
            end
        end
    end
    x = ifft2(f) * sqrt(prod(size(f)));
end

function lst = freq_pat(H,pat)
% freq_pat.m - generate frequency partition
%
% input:
%   H:    half length of the frequency span, H needs to be greater than 16
%   pat:  type of frequency partition which satsifies parabolic scaling relationship
%         equal to 'p' or 'q'
%
% output:
%   lst:  data representing the partition of frequency
%
% Written by Lexing Ying and Laurent Demanet, 2006

if(H<16)    error('H needs to be at least 16');  end

if(    pat=='p') %parabolic scaling by lexing
    len = ceil(log2(H)/2)+1;
    lst = cell(len,1);
    lst{1} = [1 1];    lst{2} = [0 1];    lst{3} = [0 1 1 1];
    cnt = 16;
    idx = 3;    rad = 4;
    while(cnt<H)
        old = cnt;
        lst{idx} = [lst{idx} 1 1];
        cnt = cnt + 2*rad;
        idx = idx+1;
        rad = 2*rad;
        trg = min(4*old,H);
        lst{idx} = [zeros(1,cnt/rad), ones(1,(trg-cnt)/rad)];
        cnt = trg;
    end
elseif(pat=='q') %parabolic scaling by laurent
    len = floor(log2(H)/2)+1;
    lst = cell(len,1);
    lst{1} = [1 1];    lst{2} = [0 1 1 1];
    cnt = 8;
    idx = 2;      rad = 2;
    while(cnt<H)
        old = cnt;
        lst{idx} = [lst{idx} 1 1];
        cnt = cnt + 2*rad;
        idx = idx+1;
        rad = 2*rad;
        trg = min(4*old,H);
        lst{idx} = [zeros(1,cnt/rad), ones(1,(trg-cnt)/rad)];
        cnt = trg;
    end
elseif(pat=='u')
    len = ceil(log2(H)/2)+1; %uniform partitioning
    lst = cell(len,1);
    B = 2^(len-1);
    lst{end} = ones(1,H/B);
else
    error('wrong pat');
end

function r = kf_rt(w,n)
% kf_rt.m - right bump
%
% Written by Lexing Ying and Laurent Demanet, 2006

r = zeros(size(w));
an = pi/2*(n+1/2);
en = (-1)^n;
en1 = (-1)^(n+1);

r = exp(-i*w/2) .* ( exp(i*an)*g(en*(w-pi*(n+1/2))) );

function r = g(w)
% g.m - 'g' function Villemoes's construction
%
% Written by Lexing Ying and Laurent Demanet, 2006
r = zeros(size(w));
gd = w<5*pi/6 & w>-7*pi/6;
r(gd) = abs(sf(w(gd)-3*pi/2));

%----------------------------------------------------------------------
function r = sf(w)
r = zeros(size(w));

aw = abs(w);
r(aw<=2*pi/3) = 0;

gd = aw>=2*pi/3 & aw<=4*pi/3;
r(gd) = 1/sqrt(2) * hf(w(gd)/2+pi);

gd = aw>=4*pi/3 & aw<=8*pi/3;
r(gd) = 1/sqrt(2) * hf(w(gd)/4);

r(aw>8*pi/2) = 0;

function r = hf(w)
w = mod(w+pi,2*pi) - pi;
r = zeros(size(w));
w = abs(w);

r = sqrt(2) * cos(pi/2 * beta(3*w/pi-1));
r(w<=pi/3) = sqrt(2);
r(w>=2*pi/3) = 0;

function r = beta(x)
r = x.^4 .*(35-84*x + 70*x.^2 -20*x.^3);

function r = kf_lf(w,n)
% kf_lf.m - left bump
%
% Written by Lexing Ying and Laurent Demanet, 2006

r = zeros(size(w));
an = pi/2*(n+1/2);
en = (-1)^n;
en1 = (-1)^(n+1);

r = exp(-i*w/2) .* ( exp(-i*an)*g(en1*(w+pi*(n+1/2))) );

function res = fwa2sym(x,pat,tp)
% fwa2sym - 2D forward wave atom transform (symmetric version)
% -----------------
% INPUT
% --
% x is a real N-by-N matrix. N is a power of 2.
% --
% pat specifies the type of frequency partition which satsifies
% parabolic scaling relationship. pat can either be 'p' or 'q'.
% --
% tp is the type of tranform.
% 	'ortho': orthobasis
% 	'directional': real-valued frame with single oscillation direction
% 	'complex': complex-valued frame
% -----------------
% OUTPUT
% --
% res is an array containing all the wave atom coefficients. If
% tp=='ortho', then res is of size N-by-N. If tp=='
% If tp=='directional', then res is of size N-by-N-by-2. If tp=='complex',
% then res is of size N-by-N-by-4.
% -----------------
% Written by Lexing Ying and Laurent Demanet, 2007

if( ismember(tp, {'ortho','directional','complex'})==0 | ismember(pat, {'p','q','u'})==0 )    error('wrong');  end

call = fwa2(x,pat,tp);
N = size(x,1);
res = zeros(N,N,size(call,2));

for i=1:size(call,2)
    c = call(:,i);
    y = zeros(N,N);
    for s=1:length(c)
        D = 2^s;
        nw = length(c{s});
        for I=0:nw-1
            for J=0:nw-1
                if(~isempty(c{s}{I+1,J+1}))
                    y( I*D+[1:D], J*D+[1:D] ) = c{s}{I+1,J+1};
                end
            end
        end
    end
    res(:,:,i) = y;
end

function x = iwa2sym(res,pat,tp)
% iwa2sym - 2D inverse wave atom transform (symmetric version)
% -----------------
% INPUT
% --
% res is an array containing all the wave atom coefficients. If
% tp=='ortho', then res is of size N-by-N. If tp=='
% If tp=='directional', then res is of size N-by-N-by-2. If tp=='complex',
% then res is of size N-by-N-by-4.
% --
% pat specifies the type of frequency partition which satsifies
% parabolic scaling relationship. pat can either be 'p' or 'q'.
% --
% tp is the type of tranform.
% 	'ortho': orthobasis
% 	'directional': real-valued frame with single oscillation direction
% 	'complex': complex-valued frame
% -----------------
% OUTPUT
% --
% x is a real N-by-N matrix. N is a power of 2.
% -----------------
% Written by Lexing Ying and Laurent Demanet, 2007

if( ismember(tp, {'ortho','directional','complex'})==0 | ismember(pat, {'p','q','u'})==0 )    error('wrong');  end

N = size(res,1);
H = N/2;
lst = freq_pat(H,pat);
red = size(res,3);

call = cell(length(lst),red);
for i=1:red
    y = res(:,:,i);
    c = cell(length(lst),1);
    for s=1:length(lst)
        B = 2^(s-1);    D = 2*B;
        nw = length(lst{s});
        c{s} = cell(nw,nw);
        for I=0:nw-1
            for J=0:nw-1
                if(lst{s}(I+1)==0 & lst{s}(J+1)==0)
                    c{s}{I+1,J+1} = [];
                else
                    c{s}{I+1,J+1} = y( I*D+[1:D], J*D+[1:D] );
                end
            end
        end
    end
    call(:,i) = c;
end
x = iwa2(call,pat,tp);

author: gpeyre

repo: matlab-toolboxes


text:
function funUpload1(object_handle, event)

    global mode;
    global figure1;

    mode = 'Upload1';
    figure1 = true;

    set(object_handle, 'BackGroundColor', [0.4,0.4,0.4], 'Enable', 'off');

end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
function [result, result_raw] = par_evaluate(task__,noresult)
% Internal: Task processing function of the worker.
% Result = par_evaluate(Task)
%
% Both the task and the result are assumed to be in serialized form, and this function deserializes
% the task, evaluates it, and re-serializes the results.
%
% In:
%   Task : task record; A sequence of bytes which represents a cell array of the form:
%          {Number, Function, Arugment1, Arguments2, ...} where Number is an arbitrary task
%          number, Function is the function to evaluate, and ArgumentX are arguments to the function.
%
%   NoResult : Whether the task does not produce a result. (default: false)
%
% Out:
%   Result : result record; A sequence of bytes which represents a cell array of the form:
%            {Number, ResultData} where Number is the corresponding task number and ResultData is
%            the first output of the function.
%
%            If an exception happened during processing, result represents {Number, LastError} where
%            LastError is an error record as returned by lasterror.
%
% Notes:
%   Depending on the severity of the error (especially when the serialization itself fails), the
%   error record may contain an inaccurate number or just a string as error
%   record.
%
% See also:
%   par_worker
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-08-29

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

if nargin < 2
    noresult = false; end

result_raw = [];
try
    % deserialize the task (names mangled as the task may involve an "eval" call).
    task__ = hlp_deserialize(task__);
    % evaluate Function(Argument1,Argument2, ...) and construct the result record
    if noresult
        task__{2}(task__{3:end});
    else
        result_raw = {task__{1},task__{2}(task__{3:end})};
        % serialize the result record
        result = hlp_serialize(result_raw);
    end
catch e
    try
        % try to display and serialize the error
        fprintf('exception during task processing: %s\n',e.message);
        if exist('hlp_handleerror','file')
            hlp_handleerror(e);
        else
            for k = 1:length(e.stack)
                fprintf('  %s (%i)\n',e.stack(k).name,e.stack(k).line); end
        end
        if iscell(task__)
            result = hlp_serialize({task__{1},e});
        else
            result = hlp_serialize({NaN,e});
        end
    catch
        % fall back to minimal reporting
        try
            fprintf('exception during error serialization: %s\n',e.message);
            result = hlp_serialize({NaN,['Error during error serialization: ' e.message]});
        catch
            fprintf('serialization failed completely.\n');
            result = uint8(200);
        end
    end
end

repo: REST