text:
% 2021 EDIT: This function is now deprecated. The current version of the shaped PRESS
% simulation is run_simPressShaped.m. The current version employs coherence selection,
% resulting in a 4x simulation speed increase.
%
% run_simPressShaped_fast_phCyc.m
% Jamie Near, McGill University 2018.
%
% USAGE:
% out=run_simPressShaped_fast_phCyc(spinSys);
%
% DESCRIPTION:
% This script simulates a PRESS experiment with fully shaped refocusing
% pulses.  Phase cycling of refocusing pulses is performed.  Furthermore,
% simulations are run at various locations in space to account for the
% within-voxel spatial variation of the metabolite signal.  Summation
% across phase cycles and spatial positions is performed.  To achieve
% faster perfomance compared to the original 'run_simPressShaped.m' function,
% this code uses the method described by Yan Zhang et al. Med Phys 2017;44(8):
% 4169-78.  Some additional acceleration is currently performed using parfor
% loops in both x and y directions.  To enable the use of the MATLAB
% parallel computing toolbox, initialize the multiple worked nodes using
% "matlabpool size X" where "X" is the number of available processing
% nodes.  If the parallel processing toolbox is not available, then replace
% the "parfor" loop with a "for" loop.
%
% INPUTS:
% To run this script, there is technically only one input argument:
% spinSys           = spin system to simulate
%
% However, the user should also edit the following parameters as
% desired before running the function:
% refocWaveform     = name of refocusing pulse waveform.
% refTp             = duration of refocusing pulses[ms]
% Bfield            = Magnetic field strength in [T]
% Npts              = number of spectral points
% sw                = spectral width [Hz]
% Bfield            = magnetic field strength [Tesla]
% lw                = linewidth of the output spectrum [Hz]
% thkX              = slice thickness of x refocusing pulse [cm]
% thkY              = slice thickness of y refocusing pulse [cm]
% fovX              = full simulation FOV in the x direction [cm]
% fovY              = full simulation FOV in the y direction [cm]
% nX                = number of spatial grid points to simulate in x-direction
% nY                = number of spatial grid points to simulate in y-direction
% taus              = vector of pulse sequence timings  [ms]
% refPhCyc1         = vector of phase cycling steps for 1st refocusing pulse [degrees]
% refPhCyc2         = vector of phase cycling steps for 2nd refocusing pulse [degrees]
%
% OUTPUTS:
% out               = Simulation results, summed over all space.

function out=run_simPressShaped_fast_phCyc(spinSys)

% ************INPUT PARAMETERS**********************************
refocWaveform='sampleRefocPulse.pta'; %name of refocusing pulse waveform.
refTp=3.5; %duration of refocusing pulses[ms]
flipAngle=137;  %Flip Angle of the refocusing pulses [degrees] (e.g. Use 180 for Siemens PRESS.  Use 137 for GE PRESS).
Npts=2048; %number of spectral points
sw=2000; %spectral width [Hz]
Bfield=3; %magnetic field strength [Tesla]
lw=2; %linewidth of the output spectrum [Hz]
thkX=1.66; %slice thickness of x refocusing pulse [cm]
thkY=1.66; %slice thickness of y refocusing pulse [cm]
fovX=2.4; %size of the full simulation Field of View in the x-direction [cm]
fovY=2.4; %size of the full simulation Field of View in the y-direction [cm]
nX=32; %Number of grid points to simulate in the x-direction
nY=32; %Number of grid points to simulate in the y-direction
tau1=30; %TE1 for first spin echo [ms]
tau2=105; %TE2 for second spin echo [ms]
refPhCyc1=[0,90]; %phase cycling steps for 1st refocusing pulse [degrees]
refPhCyc2=[0,90]; %phase cycling steps for 2nd refocusing pulse [degrees]
% ************END OF INPUT PARAMETERS**********************************

%set up spatial grid
x=linspace(-fovX/2,fovX/2,nX); %X positions to simulate [cm]
y=linspace(-fovY/2,fovY/2,nY); %y positions to simulate [cm]

%Load RF waveform
refRF=io_loadRFwaveform(refocWaveform,'ref',0);

gamma=42577000; %gyromagnetic ratio

%Load spin systems
load spinSystems
sys=eval(['sys' spinSys]);

%Resample refocusing RF pulse from 400 pts to 100 pts to reduce
%computational workload
refRF=rf_resample(refRF,100);

Gx=(refRF.tbw/(refTp/1000))/(gamma*thkX/10000); %[G/cm]
Gy=(refRF.tbw/(refTp/1000))/(gamma*thkY/10000); %[G/cm]

%Initialize structures:
d_temp=cell(length(x),length(refPhCyc1));
d=cell(length(refPhCyc1));

%loop through space: If you are using the parfor loops below, and you are
%using an older version of MATLAB (e.g.R2012), don't forget to initialize
%the parallel processing toolbox workers using 'matlabpool open N' (for N
%workers, 12 max).  I don't think this is necessary for newer version of
%MATLAB.

%First loop through all x-positions, simulating only the first refocusing
%pulse.
%First loop through x direction (first refoc pulse only);

%for X=1:length(x)  %Use this if you don't have the MATLAB parallel processing toolbox
parfor X=1:length(x)  %Use this if you have the MATLAB parallel processing toolbox
    for RP1=1:length(refPhCyc1)
        disp(['Executing X-position ' num2str(X) ' of ' num2str(length(x)) ', '...
            'First Refoc phase cycle ' num2str(RP1) ' of ' num2str(length(refPhCyc1)) '!!!']);
        d_temp{X}{RP1}=sim_press_shaped_fastRef1(Bfield,sys,tau1,tau2,refRF,refTp,x(X),Gx,refPhCyc1(RP1),flipAngle);
    end
end

%calculate the average density matrix (Doing this inside a separate for
%loop because I couldn't figure out how to do this inside the parfor loop):
for X=1:length(x)
    for RP1=1:length(refPhCyc1)
        d{RP1}=sim_dAdd(d{RP1},d_temp{X}{RP1});
    end
end

% %Initialize structures:
out_temp=cell(length(y),length(refPhCyc1),length(refPhCyc2));
out=struct([]);

%Now loop through y direction (second refoc pulse only);
%for Y=1:length(y) %Use this if you don't have the MATLAB parallel processing toolbox
parfor Y=1:length(y) %Use this if you do have the MATLAB parallel processing toolbox
    for RP1=1:length(refPhCyc1)
        for RP2=1:length(refPhCyc2)
            disp(['Executing Y-position ' num2str(Y) ' of ' num2str(length(y)) ', '...
                'First Refoc phase cycle ' num2str(RP1) ' of ' num2str(length(refPhCyc1)) ', '...
                'Second Refoc phase cycle ' num2str(RP2) ' of ' num2str(length(refPhCyc2)) '!!!']);
            out_temp{Y}{RP1}{RP2}=sim_press_shaped_fastRef2(d{RP1},Npts,sw,Bfield,lw,sys,tau1,tau2,...
                refRF,refTp,y(Y),Gy,refPhCyc2(RP2),flipAngle);
        end
    end
end

%Now combine the outputs;  Again, doing this inside a separate for loop
%becuase I can't figure out how to do this inside the parfor loop:
for Y=1:length(y)
    for RP1=1:length(refPhCyc1)
        for RP2=1:length(refPhCyc2)
            out=op_addScans(out,out_temp{Y}{RP1}{RP2},xor(RP1-1,RP2-1));
        end
    end
end

%For consistent scaling across different shaped simulations, we need to :
%1.  Scale down by the total number of simulations run (since these were
%    all added together.
numSims=(nX*nY*length(refPhCyc1)*length(refPhCyc2));
out=op_ampScale(out,1/numSims);

%2.  Scale by the total size of the simulated region, relative to the size
%    of the voxel.
voxRatio=(thkX*thkY)/(fovX*fovY);
out=op_ampScale(out,1/voxRatio);

end

%Nested Function #1
function d = sim_press_shaped_fastRef1(Bfield,sys,tau1,tau2,RF,tp,dx,Gx,phCyc1,flipAngle)
%
% USAGE:
% d = sim_press_shaped_fastRef1(n,sw,Bfield,linewidth,sys,tau1,tau2,RF,tp,dx,Gx,phCyc1,flipAngle)
%
% DESCRIPTION:
% This function simulates only the first bit of the PRESS experiment, up to
% the beginning of the second refocusing pulse.  The excitation is
% simulated as an instantaneous rotation, and the refocusing pulse is
% simulated as a shaped rotation.
%
% This code is designed to be used in highly-accelerated shaped simulations,
% using the method described by Yan Zhang et al. Med Phys 2017;44(8):
% 4169-78.
%
% This code enables the choice of the phase of the refocusing pulse.  This
% enables phase cycling of the refocusing pulses by repeating simulations
% with different editing pulse phases, which is necessary to remove phase
% artefacts from the editing pulses.  A four step phase cycling scheme is typically
% sufficient, where both refocusing pulses are phase cycled by 0 and 90 degrees, and
% the phase are combined in the following way:
%
% signal = ([0 90] - [0 0]) + ([90 0] - [90 90]);
%
% where, in [X Y], X is the phase of the first refocusing pulse and Y is
% the phase of the second refocusing pulse
%
% Finally, this code simulates the spectrum at a given point in space (x),
% given the values of the slice selection gradient (Gx).  In order
% to fully simulate the MEGA-PRESS experiment, you have to run this
% simulation many times at various points in space (x), followed by
% sim_press_shaped_fastRef2.m, at all points in space (y).
%
% INPUTS:
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% sys       = spin system definition structure
% tau1      = echo time 1 in [ms].
% tau2      = echo time 2 in [ms].
% RF        = RF pulse definition structure for refoc pulses (obtain using 'io_loadRFwaveform.m')
% tp        = RF pulse duration in [ms]
% dx        = position offset in x-direction (corresponding to first refocusing pulse) [cm]
% dy        = position offset in y-direction (corresponding to second refocusing pulse) [cm]
% Gx        = gradient strength for first selective refocusing pulse [G/cm]
% Gy        = gradient strength for second selective refocusing pulse [G/cm]
% phCycl    = initial phase of the first refocusing pulse in [degrees];
% phCycl2   = initial phase of the second refocusing pulse in [degrees];
% flipAngle = flip angle of refocusing pulses [degrees] (Optional.  Default = 180 deg)
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using PRESS
%             sequence.

if nargin<10
    flipAngle=180;
end

if tau1<tp/1000
    error('ERROR:  Echo-time 1 cannot be less than duration of refocusing pulse! ABORTING!!');
end
if tau2<tp/1000
    error('ERROR:  Echo-time 2 cannot be less than duration of refocusing pulse! ABORTING!!');
end

%Set water to centre
centreFreq=2.3;
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H,d]=sim_Hamiltonian(sys,Bfield);

%Calculate new delays by subtracting the pulse duration from tau1 and tau2;
delays=zeros(2);
delays(1)=tau1-tp;
delays(2)=tau2-tp;
if sum(delays<0)
    error(['ERROR! The following taus are too short: ' num2str(find(delays<0)) '.']);
end

%BEGIN PULSE SEQUENCE************
d=sim_excite(d,H,'x');                                    %EXCITE
d=sim_evolve(d,H,delays(1)/2000);                            %Evolve by delays(1)/2
d=sim_shapedRF(d,H,RF,tp,flipAngle,90+phCyc1,dx,Gx);          %1st shaped 180 degree refocusing pulse
d=sim_evolve(d,H,(delays(1)+delays(2))/2000);                     %Evolve by (delays(1)+delays(2))/2
%END PULSE SEQUENCE**************

%After running this many times along x, the density matrices should be
%averaged, and then the average density matrix should be passed through
%'sim_press_shaped_fastRef2' at various different y-positions.

end

%Nested Function #2
function out = sim_press_shaped_fastRef2(d,n,sw,Bfield,linewidth,sys,tau1,tau2,RF,tp,dy,Gy,phCyc2,flipAngle)
%
% USAGE:
% out = sim_press_shaped_fastRef2(d,n,sw,Bfield,linewidth,sys,tau2,RF,tp,dy,Gy,phCyc2,flipAngle)
%
% DESCRIPTION:
% This function simulates only the last bit of the PRESS experiment, from the
% the beginning of the second refocusing pulse, to the end.  The refocusing
%pulse is simulated as a shaped rotation.
%
% This code is designed to be used in highly-accelerated shaped simulations,
% using the method described by Yan Zhang et al. Med Phys 2017;44(8):
% 4169-78.
%
% This code enables the choice of the phase of the refocusing pulse.  This
% enables phase cycling of the refocusing pulses by repeating simulations
% with different editing pulse phases, which is necessary to remove phase
% artefacts from the editing pulses.  A four step phase cycling scheme is typically
% sufficient, where both refocusing pulses are phase cycled by 0 and 90 degrees, and
% the phase are combined in the following way:
%
% signal = ([0 90] - [0 0]) + ([90 0] - [90 90]);
%
% where, in [X Y], X is the phase of the first refocusing pulse and Y is
% the phase of the second refocusing pulse
%
% Finally, this code simulates the spectrum at a given point in space (y),
% given the values of the slice selection gradient (Gy).  In order
% to fully simulate the MEGA-PRESS experiment, you have to first run
% sim_press_shaped_fastRef1.m at all points in space (x), followed by
% this code, at all points in space (y).
%
% INPUTS:
% d         = starting density matrix (obtained using 'sim_press_shaped_fastRef1.m')
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% sys       = spin system definition structure
% tau1      = echo time 1 in [ms].
% tau2      = echo time 2 in [ms].
% RF        = RF pulse definition structure for refoc pulses (obtain using 'io_loadRFwaveform.m')
% tp        = RF pulse duration in [ms]
% dx        = position offset in x-direction (corresponding to first refocusing pulse) [cm]
% dy        = position offset in y-direction (corresponding to second refocusing pulse) [cm]
% Gx        = gradient strength for first selective refocusing pulse [G/cm]
% Gy        = gradient strength for second selective refocusing pulse [G/cm]
% phCycl    = initial phase of the first refocusing pulse in [degrees];
% phCycl2   = initial phase of the second refocusing pulse in [degrees];
% flipAngle = flip angle of refocusing pulses [degrees] (Optional.  Default = 180 deg)
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using PRESS
%             sequence.

if nargin<14
    flipAngle=180;
end

if tau1<tp/1000
    error('ERROR:  Echo-time 1 cannot be less than duration of refocusing pulse! ABORTING!!');
end
if tau2<tp/1000
    error('ERROR:  Echo-time 2 cannot be less than duration of refocusing pulse! ABORTING!!');
end

%Set water to centre
centreFreq=2.3;
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H]=sim_Hamiltonian(sys,Bfield);

%Calculate new delays by subtracting the pulse duration from tau1 and tau2;
delays=zeros(2);
delays(1)=tau1-tp;
delays(2)=tau2-tp;
if sum(delays<0)
    error(['ERROR! The following taus are too short: ' num2str(find(delays<0)) '.']);
end

%BEGIN PULSE SEQUENCE************
d=sim_shapedRF(d,H,RF,tp,flipAngle,90+phCyc2,dy,Gy);          %2nd shaped 180 degree refocusing pulse
d=sim_evolve(d,H,delays(2)/2000);                            %Evolve by delays(2)/2
[out,~]=sim_readout(d,H,n,sw,linewidth,90);      %Readout along y (90 degree phase);
%END PULSE SEQUENCE**************

%Correct the ppm scale:
out.ppm=out.ppm-(4.65-centreFreq);

%Fill in structure header fields:
out.seq='press';
out.te=tau1+tau2;
out.sim='shaped';

%Additional fields for compatibility with FID-A processing tools.
out.sz=size(out.specs);
out.date=date;
out.dims.t=1;
out.dims.coils=0;
out.dims.averages=0;
out.dims.subSpecs=0;
out.dims.extras=0;
out.averages=1;
out.rawAverages=1;
out.subspecs=1;
out.rawSubspecs=1;
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=1;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;

end

author: CIC-methods

repo: FID-A

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


text:
function ctf = sortrtctf(ctfds)

RESPTHRESH = 10^7;

if nargin < 1
    ctfds = '';
end

ctf = ctf_read(ctfds);

stim_chan = find(strcmp(ctf.sensor.label,'STIM'));

if (length(stim_chan) == 1)
    stim_data = ctf.data(:,stim_chan,:);
else
    error('ERROR: could not properly locate STIM channel');
end

RTs = [];
for i = 1:size(stim_data,3)
    tmp = stim_data(:,1,i);
    RTs(end+1) = min(find(tmp > RESPTHRESH));
end

[RTs_sorted,RTsx] = sort(RTs);

ctf.data = ctf.data(:,:,RTsx);

author: PatternRecognition

repo: OpenBMI


text:
%______________________________________________________________________
function [y,C,R,c1]  = MissData(y,C,R,c1);
%______________________________________________________________________
% PROC missdata
% PURPOSE: eliminates the rows in y & matrices Z, G that correspond to
%          missing data (NaN) in y
% INPUT    y             vector of observations at time t  (n x 1 )
%          S             KF system matrices             (structure)
%                        must contain Z & G
% OUTPUT   y             vector of observations (reduced)   (# x 1)
%          Z G           KF system matrices     (reduced)   (# x ?)
%          L             To restore standard dimensions     (n x #)
%                        where # is the nr of available data in y
%______________________________________________________________________
  ix = ~isnan(y);
  e  = eye(size(y,1));
  L  = e(:,ix);

  y  =  y(ix);
  c1 =  c1(ix);
  C  =  C(ix,:);
  R  =  R(ix,ix);

author: european-central-bank

repo: BEAR-toolbox


text:
function result = SimulateSSCOMP_3Spaces(D, K, Ng, theta, SNR, shuffle)
% Simulates SSC-OMP algorithm
% D - Ambient space dimension
% K - Subspace dimension
% Ng - Number of vectors per subspace
% theta - angle between disjoint subspaces
% shuffle - whether to shuffle the vectors or not.

    if nargin < 6
        shuffle = false;
    end
    if nargin < 5
        SNR = Inf;
    end
    % number of subspaces = number of clusters
    ns = 3;
    % Let us form the subspaces
    [A, B, C] = spx.la.spaces.three_disjoint_spaces_at_angle(K, deg2rad(theta));
    % Put them together
    X = [A B C];
    % Put them to bigger dimension
    X = spx.la.spaces.k_dim_to_n_dim(X, D);
    % Perform a random orthonormal transformation
    O = orth(randn(D));
    X = O * X;
    % Split them again
    A = X(:, 1:K);
    B = X(:, K + (1:K));
    C = X(:, 2*K + (1:K));
    % coefficients for Ng vectors chosen randomly in subspace A
    coeffs_a = randn(K,Ng);
    % avoid small values
    coeffs_a = 2*sign(coeffs_a) + coeffs_a;
    % coefficients for Ng vectors chosen randomly in subspace B
    coeffs_b = randn(K,Ng);
    % coefficients for Ng vectors chosen randomly in subspace C
    coeffs_c = randn(K,Ng);
    % avoid small values
    coeffs_c = 2*sign(coeffs_c) + coeffs_c;
    % Actual vectors from the three subspaces
    XA = A * coeffs_a;
    XB = B * coeffs_b;
    XC = C * coeffs_c;
    % Prepare the overall set of signals
    X = [XA XB XC];
    % ground through clustering data
    true_labels = [1*ones(Ng,1) ; 2*ones(Ng,1) ; 3*ones(Ng,1)];
    % All signals are expected to  have a K-sparse representation
    if shuffle
        % We also shuffle the signals
        shuffled_indices = randperm(3*Ng);
        X = X(:, shuffled_indices);
        true_labels = true_labels(shuffled_indices);
    end
    if ~isinf(SNR)
        % We need to add some noise in the signals.
        noises = spx.data.noise.Basic.createNoise(X, SNR);
        % Preserve original data
        X0 = X;
        % Add noise
        X = X0 + noises;
    end
    % Keep data for reference.
    result.num_subspaces = 3;
    result.theta = theta;
    result.subspace_dimension = K;
    result.ambient_dimension = D;
    result.num_signals_per_subspace = Ng;
    result.num_total_signals = size(X, 2);
    result.true_labels = true_labels;
    result.X = X;

    tstart = tic;
    % Application of Sparse subspace clustering
    solver = spx.cluster.ssc.SSC_OMP(X, K, ns);
    solver.Quiet = true;
    ssc_result = solver.solve();
    elapsed_time = toc(tstart);
    result.elapsed_time = elapsed_time;
    result.C = solver.Representation;
    % fprintf('Sparse subspace clustering time spent: %.2f seconds\n', elapsed_time);

    cluster_labels = ssc_result.Labels;
    %combined_labels = [true_labels cluster_labels]';
    result.cluster_labels = cluster_labels;

    % Time to compare the clustering
    comparer = spx.cluster.ClusterComparison(true_labels, cluster_labels);
    result.comparison = comparer.fMeasure();
    % fprintf('Sparse subspace clustering results:\n');
    % comparer.printF1MeasureResult(result.comparison);
end

author: indigits

repo: sparse-plex


text:
function [curYawTheta, curHeight] = get_hdc_initial_value()
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

    % set the initial position in the hdcell network
%     global YAW_HEIGHT_HDC_Y_DIM;
%     global YAW_HEIGHT_HDC_H_DIM;
%     curYawTheta = floor(YAW_HEIGHT_HDC_Y_DIM / 2);  % in 1:36
    curYawTheta = 1;  % in 1:36

%     curHeight = floor(YAW_HEIGHT_HDC_H_DIM / 2);  % in 1:36
    curHeight = 1;  % in 1:36
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function boundaryBenchGraphs(pbDir)
% function boundaryBenchGraphs(pbDir)
%
% Create graphs, after boundaryBench(pbDir) has been run.
%
% See also boundaryBench.
%
% David Martin <dmartin@eecs.berkeley.edu>
% May 2003

fname = fullfile(pbDir,'scores.txt');
scores = dlmread(fname); % iid,thresh,r,p,f
fname = fullfile(pbDir,'score.txt');
score = dlmread(fname); % thresh,r,p,f

% create the overall PR graph
fname = fullfile(pbDir,'pr.txt');
pr = dlmread(fname); % thresh,r,p,f
h = figure(1); clf; hold on;
prplot(h,pr(:,2),pr(:,3),sprintf('F=%4.2f',score(4)));
ph=plot(score(2),score(3),'ko','MarkerFaceColor','k','MarkerSize',10);
lh=legend(ph,sprintf('F=%4.2f @(%4.2f,%4.2f) t=%4.2f',...
                     score(4),score(2),score(3),score(1)));
p=get(lh,'Position');
p(1)=0.25;
p(2)=0.15;
set(lh,'Position',p);
print(h,'-depsc2',fullfile(pbDir,'pr.eps'));
print(h,'-djpeg95','-r36',fullfile(pbDir,'pr_half.jpg'));
print(h,'-djpeg95','-r0',fullfile(pbDir,'pr_full.jpg'));
%close(h);

iids = imgList('test');
for i = 1:numel(iids),
  iid = iids(i);
  fprintf(2,'Processing image %d/%d (iid=%d)...\n',i,numel(iids),iid);

  % create PR graphs for this image
  fname = fullfile(pbDir,sprintf('%d_pr.txt',iid));
  pri = dlmread(fname);
  h = figure(1); clf; hold on;
  prplot(h,pri(:,2),pri(:,3),sprintf('%d F=%4.2f',iid,scores(i,5)));
  ph=plot(scores(i,3),scores(i,4),'ko','MarkerFaceColor','k','MarkerSize',10);
  lh=legend(ph,sprintf('F=%4.2f @(%4.2f,%4.2f) t=%4.2f',...
                       scores(i,5),scores(i,3),scores(i,4),scores(i,2)));
  p=get(lh,'Position');
  p(1)=0.25;
  p(2)=0.15;
  set(lh,'Position',p);
  print(h,'-depsc2',fullfile(pbDir,sprintf('%d_pr.eps',iid)));
  print(h,'-djpeg95','-r36',fullfile(pbDir,sprintf('%d_pr_half.jpg',iid)));
  print(h,'-djpeg95','-r0',fullfile(pbDir,sprintf('%d_pr_full.jpg',iid)));
  %close(h);
end

function prplot(h,r,p,ti)
figure(h);
plot(r,p,'ko-');
box on;
grid on;
set(gca,'Fontsize',12);
set(gca,'XTick',[0 .25 .5 .75 1]);
set(gca,'YTick',[0 .25 .5 .75 1]);
set(gca,'XGrid','on');
set(gca,'YGrid','on');
xlabel('Recall');
ylabel('Precision');
title(ti);
axis square;
axis([0 1 0 1]);

author: Cloud-CV

repo: object-proposals


text:
function z = searchMaxAcq(data,xi)
%SEARCHMAXACQ Wrapper of acquisition function for SEARCHMAX

if isstruct(xi) && isnumeric(data) % Arguments in wrong order
    temp = data;
    data = xi;
    xi = temp;
end

z = feval(data.acq{:},xi(:)',data.target,data.gpstruct,data.optimState,0);

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function y = prtRvUtilNormPdf(X,mu,sigma)
%Y = prtRvUtilNormPdf(X,mu,sigma)
%Y = prtRvUtilNormPdf(X)
%Y = prtRvUtilNormPdf(X,mu)

% Test Drop in replacement of normpdf from stats toolbox
% X = randn([2 3 3]); mu = 1; sigma = 2; prtUtilApproxEqual(prtRvUtilNormPdf(X,mu,sigma),normpdf(X,mu,sigma))

if nargin < 2 %|| isempty(mu) % normpdf returns [] if mu is []
    mu = 0;
end

if nargin < 3 %|| isempty(sigma) % normpdf returns [] if sigma is []
    sigma = 1;
end

xSize = size(X);

y = reshape(exp(prtRvUtilMvnLogPdf(X(:), mu, sigma.^2)),xSize);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function varargout = rank(varargin)
%RANK   Rank of a DISKFUN.
%   RANK(F) produces an estimate of the rank of F.
%
%   RANK(F, TOL) is the number of singular values of F greater than TOL*N,
%   where N is the first singular value of F.
%
% See also DISKFUN/LENGTH.

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

[varargout{1:nargout}] = rank@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% ODE_POOL uses the MATLABPOOL command to run the ODE code.
%
%  Discussion:
%
%    Output printed by the function appears directly on the screen.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2011
%
%  Author:
%
%    John Burkardt
%

%
%  Initialize the k and b ranges.
%
  bVals = 0.1 : 0.05 : 5;
  kVals = 1.5 : 0.05 : 5;
%
%  Begin the parameter sweep.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ODE_POOL\n' );
  fprintf ( 1, '  Sweep through sets of values of parameters B and K,\n' );
  fprintf ( 1, '  computing the solution of the ODE corresponding to each set.\n' );
  fprintf ( 1, '  For each solution X(T), determine the maximum value over time.\n' );
  fprintf ( 1, '  Construct a contour plot of XMAX(B,K).\n' );
  fprintf ( 1, '  Use the PARFOR command to carry out these computations in parallel.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of K values = %d\n', length ( kVals ) );
  fprintf ( 1, '  Number of B values = %d\n', length ( bVals ) );
  fprintf ( 1, '  Number of times the ODE must be solved = %d\n', ...
    length ( kVals ) * length ( bVals ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Begin computation\n' );

  matlabpool open local 4
%
%  Solve the ODE for every pair of K and B values and return the maximum
%  value over the time interval.
%
  tic
  peakVals = ode_fun ( bVals, kVals );
  toc

  matlabpool close
%
%  Now display am image of the data.
%
  ode_display ( bVals, kVals, peakVals );

author: johannesgerer

repo: jburkardt-m


text:
function s = ndf(f)
%NDF   Number of degrees of freedom (parameters) needed to represent a
%   CHEBFUN3 object.
%
% See also CHEBFUN3T/NDF.

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

% NDF(f) = sum of modal ranks multiplied with length of the corresponding
% factor quasimatrix plus number of entries in the core tensor.

if ( isempty(f) )
    s = 0;
else
    [r1, r2, r3] = rank(f);
    [m, n, p] = length(f);
    s = dot([r1, r2, r3], [m, n, p]) + numel(f.core);
end

end


text:
function bool = prtUtilIsLogicalScalar(value)
%bool = prtUtilIsLogicalScalar(value)

bool = islogical(value) && isscalar(value);


text:
function varargout = process_spike_triggered_average( varargin )
% PROCESS_SPIKE_TRIGGERED_AVERAGE: Computes the spike triggered average.
% Select a time window around the spikes of a specific neuron and average the LFPs of each electrode

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription()
    % Description the process
    sProcess.Comment     = 'Spike triggered average';
    sProcess.FileTag     = 'STA';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Electrophysiology';
    sProcess.Index       = 1230;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/e-phys/functions#Spike_triggered_average';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data'};
    sProcess.OutputTypes = {'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 2;
    % Options: Segment around spike
    sProcess.options.timewindow.Comment  = 'Spike time window: ';
    sProcess.options.timewindow.Type     = 'range';
    sProcess.options.timewindow.Value    = {[-0.150, 0.150],'ms',[]};
    % Options: Parallel Processing
    sProcess.options.parallel.Comment = 'Parallel processing';
    sProcess.options.parallel.Type    = 'checkbox';
    sProcess.options.parallel.Value   = 0;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess)
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs)
    % Initialize returned values
    OutputFiles = {};
    % Get options
    isParallel = sProcess.options.parallel.Value;
    TimeWindow = sProcess.options.timewindow.Value{1};

    % ===== LOAD INPUTS =====
    % Loads all the data outside of the parfor, so it doesn't fail
    nTrials = length(sInputs);
    DataMats = cell(1, nTrials);
    ChannelFlag = [];
    for iFile = 1:length(sInputs)
        DataMats{iFile} = in_bst_data(sInputs(iFile).FileName);
        if isempty(ChannelFlag)
            ChannelFlag = DataMats{iFile}.ChannelFlag;
        else
            ChannelFlag(DataMats{iFile}.ChannelFlag == -1) = -1;
        end
    end
    % Check time window
    if TimeWindow(1)>=0 || TimeWindow(2)<=0
        bst_report('Error', sProcess, sInputs, 'The time-selection must be around the spikes.');
        return;
    elseif (TimeWindow(1) <= DataMats{1}.Time(1)) && (TimeWindow(2) >= DataMats{1}.Time(end))
        bst_report('Error', sProcess, sInputs, 'The spike window has to be smaller than the trial window.');
        return;
    end
    % Sampling frequency
    sampling_rate = round(abs(1. / (DataMats{1}.Time(2) - DataMats{1}.Time(1))));
    % Load channel file
    ChannelMat = in_bst_channel(sInputs(1).ChannelFile);

    % === START COMPUTATION ===
    % Input time window
    time_segmentAroundSpikes = linspace(TimeWindow(1), TimeWindow(2), abs(TimeWindow(2))* sampling_rate + abs(TimeWindow(1))* sampling_rate + 1);
    % Get LPFs
    LFP_trials = cell(1, nTrials);
    if isParallel
        parfor iFile = 1:nTrials
            LFP_trials{iFile} = get_LFPs(DataMats{iFile}, ChannelMat, TimeWindow, time_segmentAroundSpikes, sampling_rate);
        end
    else
        for iFile = 1:nTrials
            LFP_trials{iFile} = get_LFPs(DataMats{iFile}, ChannelMat, TimeWindow, time_segmentAroundSpikes, sampling_rate);
        end
    end

    % ===== COMPUTE SPIKE TRIGGERED AVERAGE =====
    % The Spike Triggered Average should be a 3d matrix
    % Number of neurons x Frequencies x Electrodes
    % Ultimately the user will select the NEURON that wants to be displayed,
    % and a 2D image with the other two dimensions will appear, showing the
    % coherence of the spikes of that neuron with the LFPs on every
    % electrode on all frequencies.

    % Create a cell that holds all of the labels and one for the unique labels
    % This will be used to take the averages using the appropriate indices
    all_labels = {};
    labelsNeurons = {}; % Unique neuron labels (each trial might have different number of neurons). We need everything that appears.
    for iFile = 1:nTrials
        for iNeuron = 1:length(LFP_trials{iFile})
            all_labels{iNeuron,iFile} = LFP_trials{iFile}(iNeuron).label;
            labelsNeurons{end+1} = LFP_trials{iFile}(iNeuron).label;
        end
    end
    labelsNeurons = unique(labelsNeurons,'stable');

    % Compute STA per individual neuron
    for iNeuron = 1:length(labelsNeurons)
        % For each TRIAL, get the index of the label that corresponds to the appropriate neuron.
        for ii = 1:size(all_labels,1)
            for jj = 1:size(all_labels,2)
                logicalEvents(ii,jj) = strcmp(all_labels{ii,jj}, labelsNeurons{iNeuron});
            end
        end

        iEvents = zeros(size(all_labels,2),1);
        for iFile = 1:size(all_labels,2)
            temp = find(logicalEvents(:,iFile));
            if ~isempty(temp)
                iEvents(iFile) = temp;
            else
                iEvents(iFile) = 0; % This shows that that neuron didn't fire any spikes on that trial
            end
        end

        % Compute the averages of the appropriate indices
        STA_single_neuron = zeros(length(ChannelMat.Channel), length(time_segmentAroundSpikes));
        std_single_neuron = zeros(length(ChannelMat.Channel), length(time_segmentAroundSpikes));
        divideBy = 0;
        for iFile = 1:size(all_labels,2)
            if iEvents(iFile)~=0
                STA_single_neuron = STA_single_neuron + LFP_trials{iFile}(iEvents(iFile)).nSpikes * LFP_trials{iFile}(iEvents(iFile)).avgLFP; % The avgLFP are sum actually.
                divideBy = divideBy + LFP_trials{iFile}(iEvents(iFile)).nSpikes;

                % Here I have the assumption that the LFPs on all trials
                % have homogeneity in their variance (Cohen, 1988, p.67):
                % http://www.utstat.toronto.edu/~brunner/oldclass/378f16/readings/CohenPower.pdf
                % https://www.statisticshowto.datasciencecentral.com/pooled-standard-deviation/
                std_single_neuron = std_single_neuron + (LFP_trials{iFile}(iEvents(iFile)).nSpikes-1) * LFP_trials{iFile}(iEvents(iFile)).stdLFP.^2;
            end
        end
        % Divide by total number of averages
        STA_single_neuron = (STA_single_neuron./divideBy)';
        std_single_neuron = sqrt(std_single_neuron./(divideBy - size(all_labels,2)));

        % Get meaningful label from neuron name
        better_label = panel_spikes('GetChannelOfSpikeEvent', labelsNeurons{iNeuron});
        neuron = panel_spikes('GetNeuronOfSpikeEvent', labelsNeurons{iNeuron});
        if ~isempty(neuron)
            better_label = [better_label ' #' num2str(neuron)];
        end

        % ===== SAVE FILE =====
        % Prepare output file structure
        FileMat = db_template('datamat');
        FileMat.F           = STA_single_neuron';
        FileMat.Time        = time_segmentAroundSpikes;
        FileMat.Std         = 2 .* std_single_neuron; % MULTIPLY BY 2 TO GET 95% CONFIDENCE (ASSUMING NORMAL DISTRIBUTION)
        FileMat.Comment     = ['Spike Triggered Average: ' str_remove_parenth(DataMats{1}.Comment) ' (' better_label ')'];
        FileMat.DataType    = 'recordings';
        FileMat.ChannelFlag = ChannelFlag;
        FileMat.Device      = DataMats{1}.Device;
        FileMat.nAvg        = 1;
        FileMat.History     = DataMats{1}.History;

        % Add history field
        FileMat = bst_history('add', FileMat, 'compute', ['Spike Triggered Average: [' num2str(TimeWindow(1)) ', ' num2str(TimeWindow(2)) '] ms']);
        for iFile = 1:length(sInputs)
            FileMat = bst_history('add', FileMat, 'average', [' - ' sInputs(iFile).FileName]);
        end

        % Get output study
        [tmp, iTargetStudy] = bst_process('GetOutputStudy', sProcess, sInputs);
        sTargetStudy = bst_get('Study', iTargetStudy);
        % Output filename
        FileName = bst_process('GetNewFilename', bst_fileparts(sTargetStudy.FileName), 'data_STA');
        OutputFiles = {FileName};
        % Save output file and add to database
        bst_save(FileName, FileMat, 'v6');
        db_add_data(iTargetStudy, FileName, FileMat);
    end
end

%% ===== GET LFP =====
% Get the events that show neurons activity
function all = get_LFPs(trial, ChannelMat, TimeWindow, time_segmentAroundSpikes, sampling_rate)
    spikeEvents = []; % The spikeEvents variable holds the indices of the events that correspond to spikes.

    allChannelEvents = cellfun(@(x) panel_spikes('GetChannelOfSpikeEvent', x), {trial.Events.label}, 'UniformOutput', 0);
    for ielectrode = 1:length(ChannelMat.Channel)
        iEvents = find(strcmp(allChannelEvents, ChannelMat.Channel(ielectrode).Name)); % Find the index of the spike-events that correspond to that electrode (Exact string match)
        if ~isempty(iEvents)
            spikeEvents(end+1:end+length(iEvents)) = iEvents;
        end
    end

    % Get segments around each spike, FOR EACH NEURON
    all = struct();
    for iNeuron = 1:length(spikeEvents) % iNeuron is the iEvent
        % Check that the entire segment around the spikes [-150,150]ms is inside the trial segment and keep only those events
        iSel = trial.Events(spikeEvents(iNeuron)).times > trial.Time(1)   + abs(TimeWindow(1)) & ...
               trial.Events(spikeEvents(iNeuron)).times < trial.Time(end) - abs(TimeWindow(2));
        events_within_segment = round(trial.Events(spikeEvents(iNeuron)).times(iSel) .* sampling_rate);

        % Create a matrix that holds all the segments around the spike of that neuron, for all electrodes.
        allSpikeSegments_singleNeuron_singleTrial = zeros(length(events_within_segment),length(ChannelMat.Channel),length(time_segmentAroundSpikes));
        for ispike = 1:length(events_within_segment)
            allSpikeSegments_singleNeuron_singleTrial(ispike,:,:) = trial.F(:, ...
                round(abs(trial.Time(1))*sampling_rate) + events_within_segment(ispike) - round(abs(TimeWindow(1)) * sampling_rate) + 1 : ...
                round(abs(trial.Time(1))*sampling_rate) + events_within_segment(ispike) + round(abs(TimeWindow(2)) * sampling_rate) + 1);
        end

        all(iNeuron).label   = trial.Events(spikeEvents(iNeuron)).label;
        all(iNeuron).nSpikes = length(events_within_segment);
        all(iNeuron).avgLFP  = squeeze(sum(allSpikeSegments_singleNeuron_singleTrial,1));
        all(iNeuron).stdLFP  = squeeze(std(allSpikeSegments_singleNeuron_singleTrial,[],1));
        all(iNeuron).Used    = 0; % This indicates if this entry has already been used for computing the SFC (some spikes might not appear on every trial imported, so a new Neuron should be identified on a later trial).
    end

    % Check if any events had no spikes in the time-region of interest and remove them!
    %  Some spikes might be on the edges of the trial. Ultimately, the
    %  Spikes Channel i events would be considered in the STA (I mean the event group, not the events themselves),
    %  but there would be a zeroed avgLFP included. Get rid of those events
    iEventsToRemove = find([all.nSpikes]==0);
    all = all(~ismember(1:length(all),iEventsToRemove));
end

author: brainstorm-tools


text:
function ret = spm_ov_mesh(varargin)
% Mesh tool - plugin for spm_orthviews
%
% This routine is a plugin to spm_orthviews. For general help about
% spm_orthviews and plugins type
%             help spm_orthviews
% at the MATLAB prompt.
%__________________________________________________________________________
% Copyright (C) 2017-2019 Wellcome Trust Centre for Neuroimaging

% Torben Lund, Guillaume Flandin & Christian Gaser
% $Id: spm_ov_mesh.m 7678 2019-10-24 14:08:03Z guillaume $

switch lower(varargin{1})
    % Context menu and callbacks
    case 'context_menu'
        item0 = uimenu(varargin{3}, ...
            'Label', 'Mesh',...
            'UserData',...
                struct('style','','width',1,'id',varargin{2}));
        item1 = uimenu(item0,...
            'Label','Select mesh(es)...',...
            'Tag','menu',...
            'Callback',@(hObj,event) mesh_menu(hObj,event,varargin{2}));
        item2 = uimenu(item0,...
            'Label','Options');
        item2_1 = uimenu(item2,...
            'Label','Line style...',...
            'UserData',struct('str','Line style', 'field','style'),...
            'Callback',@(hObj,event) mesh_options(hObj,event,item0));
        item2_2 = uimenu(item2,...
            'Label','Line width...',...
            'UserData',struct('str','Line width', 'field','width'),...
            'Callback',@(hObj,event) mesh_options(hObj,event,item0));
        ret = item0;
    case 'display'
        if nargin < 2
            varargin{2} = spm_input('Select image', '!+1', 'e');
        end
        mesh_add(varargin{2:end});
        mesh_delete(varargin{2});
        mesh_display(varargin{2});
    case 'redraw'
        mesh_redraw(varargin{2:end});
    otherwise
end

%==========================================================================
function mesh_options(hObj,event,hM)
Finter = spm_figure('FindWin', 'Interactive');
spm_input('!DeleteInputObj',Finter);
UDm = get(hObj,'UserData');
UDc = get(hM,'UserData');
if isnumeric(UDc.(UDm.field))
    N = Inf;
    T = 'e';
else
    N = Inf;
    T = 's';
end
if T == 's' && iscell(UDc.(UDm.field))
    UDc.(UDm.field) = '';
end
UDc.(UDm.field) = spm_input(UDm.str, '!+1', T, UDc.(UDm.field), N);
if T == 's'
    if ~isempty(strfind(UDc.(UDm.field),'{'))
        UDc.(UDm.field) = eval(UDc.(UDm.field));
    end
end
set(hM,'UserData',UDc);
try, mesh_redraw(UDc.id); end

%==========================================================================
function mesh_menu(hObj,event,i,varargin)

global st

mesh_add(i);
mesh_redraw(i);

%-Display extra information in the presence of multiple meshes
if numel(st.vols{i}.mesh.meshes) > 1
    hM = findobj(st.vols{i}.ax{1}.cm,'Label','Mesh');
    UD = get(hM,'UserData');
    fprintf('\n');
    for j=1:numel(st.vols{i}.mesh.meshes)
        fprintf('%s: %s\n',UD.style{j},st.vols{i}.mesh.fname{j});
    end
end

%==========================================================================
function mesh_add(i,g)

global st

try
    m = st.vols{i}.mesh.meshes;
    f = st.vols{i}.mesh.fname;
catch
    m = [];
    f = {};
end

if nargin < 2
    [g,sts] = spm_select([1 Inf],'mesh','Select mesh(es)...');
    if ~sts, st.vols{i}.mesh.meshes = m; return; end
end
if ischar(g), fname = cellstr(g); else fname = repmat({''},1,numel(g)); end
g = export(gifti(g),'patch');
for j=1:numel(g)
    if ~isfield(g(j),'vertices') || ~isfield(g(j),'faces')
        error('Selected file must contain triangular meshes.');
    end
end

st.vols{i}.mesh.meshes = [m, g];
st.vols{i}.mesh.fname = [f; fname];

hM = findobj(st.vols{i}.ax{1}.cm,'Label','Mesh');
if numel(get(hM,'children')) < 3
    uimenu(hM,...
        'Label','Display 3D mesh(es)',...
        'Tag','tmp',...
        'Callback',@(hObj,event) mesh_render(i));
    uimenu(hM,...
        'Label','Remove mesh(es)',...
        'Tag','tmp',...
        'Callback',@(hObj,event) mesh_clear(i));
end

%==========================================================================
function mesh_clear(i)

global st

mesh_delete(i);
st.vols{i} = rmfield(st.vols{i},'mesh');
hM = findobj(st.vols{i}.ax{1}.cm,'Label','Mesh');
set(findobj(hM,'Tag','menu'),'Label','Select mesh(es)');
delete(findobj(get(hM,'Children'),'Tag','tmp'));

%==========================================================================
function mesh_display(i)

global st

if nargin < 1
    i = 1;
end
o = st.vols{i}.mesh.meshes;

hM = findobj(st.vols{i}.ax{1}.cm,'Label','Mesh');
set(findobj(hM,'Tag','menu'),'Label','Add mesh(es)');
UD = get(hM,'UserData');
linespec  = UD.style;
linewidth = UD.width;
if ~iscell(linespec)
    if ~isempty(linespec)
        disp('Wrong Linestyle format please use {''linestyle1'',''linestyle2'',...}.');
    end

    linespec = {'b-' 'g-' 'r-' 'c-' 'm-' 'y-' 'k-' 'w-' ...
        'b-.' 'g-.' 'r-.' 'c-.' 'm-.' 'y-.' 'k-.' 'w-.' ...
        'b--' 'g--' 'r--' 'c--' 'm--' 'y--' 'k--' 'w--' };
    if numel(o) > numel(linespec)
        % If there are more than 24 surfaces plot them all in red
        linespec = repmat({'r-'},1,numel(o));
    end
end
if numel(linespec) < numel(o)
    disp('Meshes with unspecified line style are displayed in black.')
    linespec(end+1:numel(o)) = {'k-'};
end
if numel(linewidth) < numel(o)
    linewidth(end+1:numel(o)) = linewidth(end);
end
UD.style = linespec;
UD.width  = linewidth;
set(hM,'UserData',UD);

% This requires more work: does not work in voxel space, ie handling of
% st.Space is incomplete
bb = st.bb;
is = inv(st.Space);
lh = {};
for n=1:numel(o)
    vert = (st.vols{i}.premul(1:3,:)*[o(n).vertices';ones(1,size(o(n).vertices,1))])';

    % 1
    set(st.vols{i}.ax{1}.ax,'NextPlot','add');
    s = spm_mesh_contour(...
        struct('faces',o(n).faces,'vertices',vert),...
        st.centre(3));
    for j=1:numel(s)
        c = is*[s(j).xdata; s(j).ydata;ones(2,size(s(j).xdata,2))];
        c(2,:) = c(2,:)-bb(1,2)+1;
        c(1,:) = c(1,:)-bb(1,1)+1;
        lh{end+1} = plot(st.vols{i}.ax{1}.ax,c(1,:),c(2,:),...
            linespec{n},'LineWidth',linewidth(n));
    end

    % 2
    set(st.vols{i}.ax{2}.ax,'NextPlot','add');
    s = spm_mesh_contour(...
        struct('faces',o(n).faces,'vertices',vert(:,[3 1 2])),...
        st.centre(2));
    for j=1:numel(s)
        c = is*[s(j).xdata;ones(1,size(s(j).xdata,2));s(j).ydata;ones(1,size(s(j).ydata,2))];
        c(3,:) = c(3,:)-bb(1,1)+1;
        c(1,:) = c(1,:)-bb(1,3)+1;
        lh{end+1} = plot(st.vols{i}.ax{2}.ax,c(3,:),c(1,:),...
            linespec{n},'LineWidth',linewidth(n));
    end

    % 3
    if st.mode == 0
        warning('Not handled.');
    else
        set(st.vols{i}.ax{3}.ax,'NextPlot','add');
        s = spm_mesh_contour(...
            struct('faces',o(n).faces,'vertices',vert(:,[2 3 1])),...
            st.centre(1));
        for j=1:numel(s)
            c = is*[ones(1,size(s(j).xdata,2));[s(j).xdata; s(j).ydata];ones(1,size(s(j).ydata,2))];
            c(2,:) = -c(2,:)+bb(2,2)+1;
            c(3,:) = c(3,:)-bb(1,3)+1;
            lh{end+1} = plot(st.vols{i}.ax{3}.ax,c(2,:),c(3,:),...
                linespec{n},'LineWidth',linewidth(n));
        end
    end

end

set(cat(1,lh{:}),'HitTest','off');

st.vols{i}.mesh.handles = lh;

%==========================================================================
function mesh_redraw(i,varargin) %i, TM0, TD, CM0, CD, SM0, SD

mesh_delete(i);
mesh_display(i);

%==========================================================================
function mesh_delete(i)

global st

try, delete(cat(1,st.vols{i}.mesh.handles{:})); end

%==========================================================================
function mesh_render(i)

global st

m = st.vols{i}.mesh.meshes;
for j=1:numel(m)
    H = spm_mesh_render(m(j));
    set(H.figure,'Name',spm_str_manip(st.vols{i}.mesh.fname{j},'a60'))
end
% data cursor should be linked to orthviews through spm_XYZreg
% see spm_ovhelper_3Dreg.m

author: spm


text:
function [structureInitS,tags] = dicomrt_d2c_voi(structureInitS,indexS,voi,xmesh,ymesh,zmesh,tags,scanUID)
% dicomrt_d2c_voi(structureInitS,indexS,voi,ct_zmesh)
%
% Convert DICOM VOI in CERR format
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check input data
[study,type,dummylabel]=dicomrt_checkinput(voi);

if strcmpi('VOI',type)~=1
    error('dicomrt_d2c_voi: input data does not have the right format. Exit now!');
else
    type='STRUCTURE';
end

% Accounting for different coordinate system between DICOM and RTOG
[study,xmesh,ymesh,zmesh,tags]=dicomrt_d2c_coordsystem(study,xmesh,ymesh,zmesh,tags);

% Get DICOM-RT toolbox dataset info
study_pointer=study{1,1};
study_array=study{2,1};

% Set parameters
n = indexS.scan;

% Scan info
for i=1:size(study_array,1); % loop through the number of VOIs
    % Fill structureInitS
    structureInitS(i).imageNumber            = tags.nimages + i;
    structureInitS(i).imageType              = type;
    structureInitS(i).caseNumber             = 1;                                         % LEAVE FOR NOW
    structureInitS(i).patientName            = [study_pointer.PatientName.FamilyName];
    try
        structureInitS(i).patientName        = [structureInitS(i).patientName,' ', ...
            study_pointer.PatientName.GivenName]; % may be not present in anonymized studies
    end
    structureInitS(i).structureName          = study_array{i,1};
    structureInitS(i).numberRepresentation   = 'CHARACTER';
    structureInitS(i).structureFormat        = 'SCAN-BASED';
    structureInitS(i).numberOfScans          = tags.nimages;
    structureInitS(i).transferProtocol       = 'DICOM';
    structureInitS(i).DICOMHeaders           = study_pointer;
    % Optional
    structureInitS(i).maximumNumberScans     = '';
    structureInitS(i).structureEdition       = '';
    structureInitS(i).unitNumber             = '';
    structureInitS(i).writer                 = '';
    structureInitS(i).dateWritten            = '';
    structureInitS(i).structureColor         = '';
    structureInitS(i).structureDescription   = '';
    structureInitS(i).studyNumberOfOrigin    = '';
    % build vector containing Z location of VOIs
    [voizdef,index] = dicomrt_getvoiz(study,i);
    for j=1:tags.nimages % loop through CT scans
        % locate position of VOI segment
        locate_voi=find(voizdef==zmesh(j));
        if isempty(locate_voi)==1
            structureInitS(i).contour(length(zmesh)-j+1).segments.points=[];
        else
            for jj=1:length(locate_voi)
                if ~isempty(study_array{i,2})% check for empty structures
                    structureInitS(i).contour(length(zmesh)-j+1).segments(jj).points=study_array{i,2}{index(locate_voi(jj))};
                else
                    structureInitS(i).contour(length(zmesh)-j+1).segments(jj).points=[];
                end
            end
        end
    end
    clear voizdef;
    structureInitS(i).rasterSegments         = '';
    structureInitS(i).DSHPoints              = '';
    structureInitS(i).orientationOfStructure = '';

    %DK
    structureInitS(i).strUID = createUID('STRUCTURE');
    if exist('scanUID')
        structureInitS(i).assocScanUID = scanUID;
    else
        structureInitS(i).assocScanUID = '';
    end
    %end DK
end

% Update nimages
tags.nimages = tags.nimages + size(study_array,1); % number of CT images

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef madgwick < handle
    methods (Static = true)

        function q = imu(q, Gyroscope, Accelerometer, SamplePeriod, Beta)

            % Normalise accelerometer measurement
            acc = Accelerometer;
            if(norm(acc) == 0), return; end   % handle NaN
            acc = acc / norm(acc);  % normalise magnitude

            % Gradient decent algorithm corrective step
            F = (quatrotate(q, [0 0 1]) - acc)';
%             F = [2*(q(2)*q(4) - q(1)*q(3)) - Accelerometer(1)
%                 2*(q(1)*q(2) + q(3)*q(4)) - Accelerometer(2)
%                 2*(0.5 - q(2)^2 - q(3)^2) - Accelerometer(3)];

            J = [-2*q(3),	2*q(4),    -2*q(1),	2*q(2)
                2*q(2),     2*q(1),     2*q(4),	2*q(3)
                0,         -4*q(2),    -4*q(3),	0    ];
            step = (J'*F);
            step = step / norm(step);	% normalise step magnitude

            % Compute rate of change of quaternion
            qd = 0.5 * quatmultiply(q, [0 Gyroscope]) - Beta * step';

            % Integrate to yield quaternion
            q = q + qd * SamplePeriod;
            q = q / norm(q); % normalise quaternion
        end

        function q = ahrs(q, Gyroscope, Accelerometer, Magnetometer, SamplePeriod, Beta)

            % Normalise accelerometer measurement
            if(norm(Accelerometer) == 0), return; end	% handle NaN
            acc = Accelerometer / norm(Accelerometer);	% normalise magnitude

            % Normalise magnetometer measurement
            if(norm(Magnetometer) == 0), return; end	% handle NaN
            mag = Magnetometer / norm(Magnetometer);	% normalise magnitude

            % Reference direction of Earth's magnetic feild
            h = quatrotate(quatconj(q), mag);
            h = [0 h];
            b = [0 norm([h(2) h(3)]) 0 h(4)];

            % Gradient decent algorithm corrective step
             F= (quatrotate(q, [0 0 1]) - acc)';
             F = [F;  (quatrotate(q, [b(2) 0 b(4)]) - mag)'];

%             F = [2*(q(2)*q(4) - q(1)*q(3)) - Accelerometer(1)
%                 2*(q(1)*q(2) + q(3)*q(4)) - Accelerometer(2)
%                 2*(0.5 - q(2)^2 - q(3)^2) - Accelerometer(3)
%                 2*b(2)*(0.5 - q(3)^2 - q(4)^2) + 2*b(4)*(q(2)*q(4) - q(1)*q(3)) - Magnetometer(1)
%                 2*b(2)*(q(2)*q(3) - q(1)*q(4)) + 2*b(4)*(q(1)*q(2) + q(3)*q(4)) - Magnetometer(2)
%                 2*b(2)*(q(1)*q(3) + q(2)*q(4)) + 2*b(4)*(0.5 - q(2)^2 - q(3)^2) - Magnetometer(3)];
            J = [-2*q(3),                 	2*q(4),                    -2*q(1),                         2*q(2)
                2*q(2),                 	2*q(1),                    	2*q(4),                         2*q(3)
                0,                         -4*q(2),                    -4*q(3),                         0
                -2*b(4)*q(3),               2*b(4)*q(4),               -4*b(2)*q(3)-2*b(4)*q(1),       -4*b(2)*q(4)+2*b(4)*q(2)
                -2*b(2)*q(4)+2*b(4)*q(2),	2*b(2)*q(3)+2*b(4)*q(1),	2*b(2)*q(2)+2*b(4)*q(4),       -2*b(2)*q(1)+2*b(4)*q(3)
                2*b(2)*q(3),                2*b(2)*q(4)-4*b(4)*q(2),	2*b(2)*q(1)-4*b(4)*q(3),        2*b(2)*q(2)];
            step = (J'*F);
            step = step / norm(step);	% normalise step magnitude

            % Compute rate of change of quaternion
            qd = 0.5 * quatmultiply(q, [0 Gyroscope]) - Beta * step';

            % Integrate to yield quaternion
            q = q + qd * SamplePeriod;
            q = q / norm(q); % normalise quaternion
        end
    end
end

author: yandld


text:
% ECEF2LLA - convert earth-centered earth-fixed (ECEF)
%            cartesian coordinates to latitude, longitude,
%            and altitude
%
% USAGE:
% [lat,lon,alt] = ecef2lla(x,y,z)
%
% lat = geodetic latitude (radians)
% lon = longitude (radians)
% alt = height above WGS84 ellipsoid (m)
% x = ECEF X-coordinate (m)
% y = ECEF Y-coordinate (m)
% z = ECEF Z-coordinate (m)
%
% Notes: (1) This function assumes the WGS84 model.
%        (2) Latitude is customary geodetic (not geocentric).
%        (3) Inputs may be scalars, vectors, or matrices of the same
%            size and shape. Outputs will have that same size and shape.
%        (4) Tested but no warranty; use at your own risk.
%        (5) Michael Kleder, April 2006

function [lat,lon,alt] = ecef2lla(x,y,z)

% WGS84 ellipsoid constants:
a = 6378137;
e = 8.1819190842622e-2;

% calculations:
b   = sqrt(a^2*(1-e^2));
ep  = sqrt((a^2-b^2)/b^2);
p   = sqrt(x.^2+y.^2);
th  = atan2(a*z,b*p);
lon = atan2(y,x);
lat = atan2((z+ep^2.*b.*sin(th).^3),(p-e^2.*a.*cos(th).^3));
N   = a./sqrt(1-e^2.*sin(lat).^2);
alt = p./cos(lat)-N;

% return lon in range [0,2*pi)
lon = mod(lon,2*pi);

% correct for numerical instability in altitude near exact poles:
% (after this correction, error is about 2 millimeters, which is about
% the same as the numerical precision of the overall function)

k=abs(x)<1 & abs(y)<1;
alt(k) = abs(z(k))-b;

return

author: Sable

repo: mcbench-benchmarks


text:
% test_itk_tri_rasterization.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Cube mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% cube mesh vertices, a bit away from voxel centres to avoid inconsistent
% results (see next section)
x = [
    0.5 -0.1 0
    0.5 -0.1 1.01
    0.5 1.1 0
    0.5 1.1 1.01
    1.505 -0.1 0
    1.505 -0.1 1.01
    1.505 1.1 0
    1.505 1.1 1.01
    ];

% apply alphashape
[~, as] = alphavol(x, 1.3);
tri = as.bnd;
clear as

% plot mesh
hold off
trimesh(tri, x(:,1), x(:,2), x(:,3))
axis([0.5 1.5 0 1 0 1])

% compute parameters of output image
res = [.2 .1 .1]; % r, c, s
origin = [0 -.5 -.5]; % x, y, z
final = [2 1.5 1.5]; % x, y, z
sz = 1 + round((final([2 1 3]) - origin([2 1 3]))./res); % r, c, s

% centers of voxels in the image
cx = linspace(origin(1), final(1), sz(2));
cy = linspace(origin(2), final(2), sz(1));
cz = linspace(origin(3), final(3), sz(3));

% rasterize the mesh
im = itk_tri_rasterization(tri, x, res, sz, origin);

% plot segmentation mask, and overlay mesh
subplot(1, 2, 1)
hold off
imagesc([origin(1) final(1)], [origin(2) final(2)], im(:, :, 6))
hold on
for I = 1:length(cx)
    plot(cx(I)*[1 1], [origin(2) final(2)])
end
for I = 1:length(cy)
    plot([origin(1) final(1)], cy(I)*[1 1])
end
xmin = min(x(:, 1));
xmax = max(x(:, 1));
ymin = min(x(:, 2));
ymax = max(x(:, 2));
plot([xmin xmax xmax xmin xmin], [ymin ymin ymax ymax ymin], 'w')
xlabel('x')
ylabel('y')
axis xy

% plot a vertical cut
subplot(1, 2, 2)
hold off
% imagesc([origin(1) final(1)], [origin(3) final(3)], squeeze(im(:, 11, :)))
imagesc([origin(2) final(2)], [origin(3) final(3)], squeeze(im(:, 11, :))')
hold on
for I = 1:length(cy)
    plot(cy(I)*[1 1], [origin(3) final(3)])
end
for I = 1:length(cz)
    plot([origin(2) final(2)], cz(I)*[1 1])
end
zmin = min(x(:, 3));
zmax = max(x(:, 3));
plot([ymin ymax ymax ymin ymin], [zmin zmin zmax zmax zmin], 'w')
xlabel('y')
ylabel('z')
axis xy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Note that voxelize offers inconsistent results if the mesh exactly goes
%% through a voxel centre
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% depending on whether the voxel is at the bottom, top, left or right of
% the segmentation, it will be or not set to 1. This can be seen using the
% following set of vertices.

% create a surface mesh that is a cube
x = [
    0.5 0 0
    0.5 0 1
    0.5 1 0
    0.5 1 1
    1.5 0 0
    1.5 0 1
    1.5 1 0
    1.5 1 1
    ];

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [x,fval,exitflag,info] = opti_lbfgsb(fun,grad,lb,ub,x0,opts)
%OPTI_LBFGSB Solve a bounded NLP using L-BFGS-B (Nocedal & Zhu)
%
%   min f(x)       subject to:   lb <= x <= ub
%    x
%
%   x = opti_lbfgsb(fun,grad,lb,ub,x0) solves a NLP where fun is the
%   objective function and grad is the gradient of the objective. lb and ub
%   are the decision variable bounds (required) and x0 is a starting guess.
%   Avoid Infinite bounds.
%
%   x = opti_lbfgsb(fun,...,x0,opts) uses opts to pass optiset options to
%   the solver.
%
%   [x,fval,exitflag,info] = opti_lbfgsb(...) returns the objective value
%   at the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR L-BFGS-B
%   See referenced BSD License

%   Copyright (C) 2012 Jonathan Currie (IPL)

if(nargin < 6), opts = optiset; end
if(nargin < 5), error('LBFGSB requires at least 5 arguments'); end

%Add in clpset settings
if(isfield(opts,'solverOpts') && ~isempty(opts.solverOpts))
    popts = lbfgsbset(opts.solverOpts);
else
    popts = lbfgsbset;
end
%Add in options from optiset
popts.maxiter = opts.maxiter;
popts.maxtime = opts.maxtime;
popts.maxfeval = opts.maxfeval;
popts.tolrfun = opts.tolrfun;
popts.iterfun = opts.iterfun;
popts.optiver = optiver;
%Setup display level
popts.display = dispLevel(opts.display);

%Ensure we have bounds
if(isempty(lb) || isempty(ub))
    error('L-BFGS-B only solves bounded NLPs - You must supply lb and ub to this function');
end
%Ensure we have a gradient
if(isempty(grad))
    error('L-BFGS-B requires a gradient function');
end

%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('L-BFGS-B requires an initial guess, x0!');
end

t = tic;
% Run L-BFGS-B
[x, fval, status, iter, feval] = lbfgsb(fun,grad,lb,ub,x0,popts);

%Collect Results
info.Iterations = iter;
info.FuncEvals = feval;
info.Time = toc(t);
info.Algorithm = 'L-BFGS-B: Limited Memory BFGS Bounded Optimization';

switch(status)
    case 0
        info.Status = 'Success';
        exitflag = 1;
    case 1
        info.Status = 'Abnormal Termination';
        exitflag = -2;
    case 2
        info.Status = 'Error On Input';
        exitflag = -3;
    case 3
        info.Status = 'Exceeded Iterations';
        exitflag = 0;
    case 4
        info.Status = 'Exceeded Maximum Fevals';
        exitflag = 0;
    case 5
        info.Status = 'Exceeded Maximum Time';
        exitflag = 0;
    case 6
        info.Status = 'User Exited';
        exitflag = -5;
    otherwise
        info.Status = 'L-BFGS-B Error';
        exitflag = -3;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function BW = im2bw_ent(IM)
% This fucntion convert an intensity image to a binary image
% by using Entropy-based method
%
% Usage: BW = im2bw_ent(IM)
% Input: IM -> Input Image
% Output: BW -> Segmented Image
%
% Example:
% IM = imread('house.jpg');
% BW = im2bw_ent(IM);
% imagesc(BW)
%
% Reference: E.R.Davies Machine Vision 3rd Edition

% size of input image
dim = size(IM);
% Change input image to gray scale
leng = length(dim);
if leng == 3
    IM = rgb2gray(IM);
end

% histgram of input image
ihis = imhist(IM);

leng = length(ihis);
para = zeros(1,leng);
for k = 2:leng-1
    % intensity of class A
    classa = ihis(1:k);
    ind = (classa==0);
    classa = classa+ind;
    clear ind
    % intensity of class B
    classb = ihis(k+1:end);
    ind = (classb==0);
    classb = classb+ind;
    clear ind
    % probability distribution of class A
    Pa = classa/(dim(1,1)*dim(1,2));
    % probability distribution of class B
    Pb = classb/(dim(1,1)*dim(1,2));
    % parameters to decide threshold
    para1 = log2(sum(Pa));
    para2 = log2(sum(Pb));
    logpa = log2(Pa);
    logpb = log2(Pb);
    para3 = -sum(Pa.*logpa)/sum(Pa);
    para4 = -sum(Pb.*logpb)/sum(Pb);
    % parameter which has to be maximized
    para(1,k) = abs(para1+para2+para3+para4);
    clear classa classb logpa logpb
end

% find threshold
[maxv,row] = max(para);
thresh = row-1;
% segment input image
BW = (IM>=thresh);


text:
function basename = bz_BasenameFromBasepath(basepath)
% gets basename from a basepath input to a function.  Uses name of the last
% directory on the path to give basename.  Takes care of the problem of the
% fact that if a fileseparator is on the end of the basepath, you have to
% take it off first.
%
% INPUTS
% basepath  is assumed to be a directory containing all data for a given
%           recording session.  It is assumed organization is
%           /whateverdirectory/basename/basename.xml and basepath is
%           equivalent to whateverdirectory/basename/.
%
% OUTPUTS
% basename  as specified above
%
% Brendon Watson, 2017

if strcmp(basepath(end),filesep);
    basepath = basepath(1:end-1);
end
[~,basename] = fileparts(basepath);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function PD = pwrDiagramPD(T, PC)
% function PD = pwrDiagramPD(T, PC)
%
% T: triangulation
% PC: power centers of triangles in T
%
% The output cell PD contains pieces of the power diagram, indexed by
% dimension. PD{mP} contains pieces of dimension zero (power centers that
% correspond to fully-dimensional pieces of the triangulation T) and PD{1}
% contains fully-dimensional regions of the power diagram (corresponding to
% vertices of the triangulation).

P = piecesPD(T);
mP = size(P,1);
PD = cell(mP,1);

for i=1:mP
    EA = edgeAttPD(T, P{i});
    for j=1:size(EA,1);
        EA{j} = PC(EA{j},:);
    end
    PD{i} = EA;
end

author: optimaltransport


text:
%
% [At,b,c,K]=readsdpa(fname)
%
% Reads in a problem in SDPA sparse format, and returns it in SeDuMi
% format.
%
% 7/20/07  Modified to handle comments and other cruft in the SDPA
% file.  In particular,
%
%    1. Initial comment lines beginning with " or * are ignored.
%    2. In the first three lines, any extraneous characters after
%       the numbers are ignored.
%    3. In the block description, any occurrences of "{", ",", "}", "(",
%       and ")" are converted into spaces.
%
function [At,b,c,K]=readsdpa(fname)
%
% Open the file for reading.
%
fid=fopen(fname,'r');
if (fid == -1),
  fprintf('file does not exist!\n');
  return;
end;
%
% Skip over the initial comments.
%
inputline=fgetl(fid);
while (inputline(1)=='"' || inputline(1)=='*')
  inputline=fgetl(fid);
end
%
% Read in m.
%
m=sscanf(inputline,'%d',1);
inputline=fgetl(fid);
%
% Read in nblocks.
%
nblocks=sscanf(inputline,'%d',1);
inputline=fgetl(fid);
%
% Read in the raw block sizes.
%
inputline=regexprep(inputline,'[\,(){}]',' ');
rawblocksizes=sscanf(inputline,'%d',nblocks);
inputline=fgetl(fid);
%
% Compute the actual block sizes, figure out block types, and figure out
% where in the vectors stuff will go.
%
blocksizes=abs(rawblocksizes);
blocktypes=zeros(nblocks,1);
for i=1:nblocks,
  if (rawblocksizes(i) < 0)
     blocktypes(i)=1;
  else
     blocktypes(i)=2;
  end;
end;
blockbases=zeros(nblocks,1);
lpbase=1;
for i=1:nblocks,
  if (blocktypes(i)==1)
    blockbases(i)=lpbase;
    lpbase=lpbase+blocksizes(i);
  end;
end;
K.l=lpbase-1;
sdpbase=lpbase;
K.s=[];
for i=1:nblocks,
  if (blocktypes(i)==2)
    blockbases(i)=sdpbase;
    sdpbase=sdpbase+blocksizes(i)^2;
    K.s=[K.s blocksizes(i)];
  end;
end;
n=sdpbase-1;
%
% Now, read in the right hand side.
%
b=zeros(m,1);
inputline=regexprep(inputline,'[\,(){}]',' ');
b=sscanf(inputline,'%le',m);
%
% Now, read in the entries in the constraints and c.
%
c=zeros(1,n);
At=sparse(n,m);
[entries,count]=fscanf(fid,'%d %d %d %d %le',[5,inf]);
count=count/5;
[entriesm,entriesn]=size(entries);
for i=1:entriesn,
  if (entries(1,i)==0),
    block=entries(2,i);
    if (blocktypes(block)==1)
%
% LP data.
%
      c(1,blockbases(block)+entries(3,i)-1)=entries(5,i);
    else
%
% SDP block entry
%
      c(1,blockbases(block)+(entries(3,i)-1)*blocksizes(block)+ ...
        entries(4,i)-1)=entries(5,i);
      c(1,blockbases(block)+(entries(4,i)-1)*blocksizes(block)+ ...
        entries(3,i)-1)=entries(5,i);
    end;
  else
%
% Constraint entry.
%
    block=entries(2,i);
    constraint=entries(1,i);
    if (blocktypes(block)==1)
%
% LP data.
%
      At(blockbases(block)+entries(3,i)-1,constraint)=entries(5,i);
    else
%
% SDP block entry
%
      At(blockbases(block)+(entries(3,i)-1)*blocksizes(block)+ ...
        entries(4,i)-1,constraint)=entries(5,i);
      At(blockbases(block)+(entries(4,i)-1)*blocksizes(block)+ ...
        entries(3,i)-1,constraint)=entries(5,i);
    end;
  end;
end;
%
% Fix up the sign of c to match SeDuMi's convention.  Also, make c
% a column vector to match SeDuMi's fromsdpa().
%
c=-c';


text:
function ex9bvp
%EX9BVP  Example 9 of the BVP tutorial.
%   This boundary value problem is the subject of Chapter 8 of
%   C.C. Lin and L.A. Segel, Mathematics Applied to Deterministic
%   Problems in the Natural Sciences, SIAM, Philadelphia, 1988.
%   The ODEs
%
%      v' = (C - 1)/n
%      C' = (vC - min(x,1))/eta
%
%   are solved on the interval [0, lambda].  The boundary conditions
%   are v(0) = 0, C(lambda) = 1, and continuity of v(x) and C(x) at
%   x = 1.  Accordingly, this is a three-point BVP that must be
%   reformulated for solution with the two-point BVP solver BVP4C.
%   This reformulation involves introducing unknowns y_1(x) for v
%   and y_2(x) for C on the interval 0 <= x <= 1 and unknowns y_3(x)
%   for v and y_4(x) for C on 1 <= x <= lambda. A new independent
%   variable is introduced for the second interval, tau =
%   (x - 1)/(lambda - 1), so that it also ranges from 0 to 1. The
%   differential equations for the four unknowns are then solved
%   on the interval [0, 1].  The continuity conditions on v and C
%   become boundary conditions on the new unknowns.  A plot of v(x)
%   and C(x) on [0, lambda] involves plotting the new unknowns over
%   the subintervals.
%
%   The quantity of most interest is the emergent osmolarity Os =
%   1/v(lambda).  The parameters are related to another parameter
%   kappa by eta = lambda^2/(n*kappa^2).  Lin and Segel develop an
%   approximate solution for Os valid for "small" n.  Here the BVP is
%   solved for a range of kappa when lambda = 2 and n = 0.005.  The
%   computed Os is compared to the approximation of Lin and Segel.

% Copyright 1999, The MathWorks, Inc.

n = 5e-2;
lambda = 2;

options = [];   % place holder

sol = bvpinit(linspace(0,1,5),[1 1 1 1]);

fprintf(' kappa    computed Os  approximate Os \n')
for kappa = 2:5
  eta = lambda^2/(n*kappa^2);

  sol = bvp4c(@ex9ode,@ex9bc,sol,options,n,lambda,eta);

  K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
  approx = 1/(1 - K2);
  computed = 1/sol.y(3,end);
  fprintf('  %2i    %10.3f    %10.3f \n',kappa,computed,approx);
end

% v and C are computed separately on 0 <= x <= 1 and 1 <= x <= lambda.
% A change of independent variable is used for the second interval,
% which must then be undone to obtain the corresponding mesh.
x = [sol.x sol.x*(lambda-1)+1];
y = [sol.y(1:2,:) sol.y(3:4,:)];

clf reset
plot(x,y(1,:),x,y(2,:),'--')
legend('v(x)','C(x)')
title('A three-point BVP.')
xlabel(['\lambda = ',num2str(lambda),', \kappa = ',num2str(kappa),'.'])
ylabel('v and C')
shg

% --------------------------------------------------------------------------

function dydx = ex9ode(x,y,n,lambda,eta)
%EX9ODE  ODE function for Example 9 of the BVP tutorial.
dydx = [ (y(2) - 1)/n
         (y(1)*y(2) - x)/eta
         (lambda - 1)*(y(4) - 1)/n
         (lambda - 1)*(y(3)*y(4) - 1)/eta ];

% --------------------------------------------------------------------------

function res = ex9bc(ya,yb,n,lambda,eta)
%EX9BC  Boundary conditions for Example 9 of the BVP tutorial.
res = [ ya(1)
        yb(4) - 1
        yb(1) - ya(3)
        yb(2) - ya(4)];


text:
function[yf,mf]=yearfrac(num)
%YEARFRAC  Converts a DATENUM into 'year.fraction' and 'month.fraction'.
%
%   YF=YEARFRAC(NUM) where NUM is an array of dates in Matlab's 'datenum'
%   format, returns the fraction of the year at each date.
%
%   FLOOR(YF) returns the year.  Note that the actual number of days in
%   each year is used, including leap years.
%
%   YF=YEARFRAC(NUM) where NUM is a cell array of numeric arrays, also
%   works.  YF is then a cell array of the same size as NUM.
%
%   [YF,MF]=YEARFRAC(NUM) also returns MF, the fraction of the current
%   month at each date.  FLOOR(MF) is the standard month number.
%
%   Note NaNs and Infs in NUM are passed through to the same values in YF.
%
%   See also DATENUM, DATEVEC.
%
%   Usage: yf=yearfrac(num);
%          [yf,mf]=yearfrac(num);
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 1998--2020 J.M. Lilly --- type 'help jlab_license' for details

if ~iscell(num)
    if strcmpi(num, '--t')
        yf2num('--t'),return
    end
end

if iscell(num)
    for i=1:length(num)
        if ~isempty(num{i})
            [yf{i,1},mf{i,1}]=yearfrac1(num{i});
        else
            yf{i}=[];
            mf{i}=[];
        end
    end
else
    [yf,mf]=yearfrac1(num);
    yf=reshape(yf,size(num));
    mf=reshape(mf,size(num));
end

function[yf,mf]=yearfrac1(num)
yf=[];
if ~isempty(num)
    index=find(isnan(num));
    if ~isempty(index)
        num(index)=0;
    end
    infindex=find(isinf(num));
    if ~isempty(infindex)
        num(infindex)=0;
    end

    [y,mo,d,h,mi,s] = datevec(num);

    %Number of days in this year?
    nd=datenum(y,12,31)-datenum(y-1,12,31);
    na=datenum(y,mo,d,h,mi,s)-datenum(y-1,12,31)-1;
    %The minus one is because for Jan 1, I add nothing to year.fraction

    yf=y+na./nd;

    %Number of days in this month?
    nd=datenum(y,mo+1,1)-datenum(y,mo,1);  %Yes, this also works if mo = 12, for January
    mf=mo+((d-1)+h/24)./nd;

    if ~isempty(index)
      yf(index)=nan;
      mf(index)=nan;
    end
    if ~isempty(infindex)
      yf(infindex)=inf;
      mf(infindex)=inf;
    end
end

function[num]=yf2num(yf)
%YF2NUM  Convert date in 'year.fraction' format to 'datenum' format.
%
%   YF2NUM(YF) where YF is an array of dates in 'year.fraction' format
%   returns the array in Matlab's 'datenum' format.
%
%   See also YEARFRAC, DATENUM, DATEVEC.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2009 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(yf, '--t')
  yf2num_test,return
end

y=floor(yf);

%Number of days in this year?
d0=datenum(y-1,12,31);
d1=datenum(y,12,31);
nd=d1-d0;
num=d0+nd.*(yf-y)+1;

function[]=yf2num_test

yearf=(1850:(1/360):2004)';
num=yf2num(yearf);
yearf2=yearfrac(num);
bool=maxmax(abs(yearf-yearf2))<1e-10;
reporttest('YF2NUM and YEARFRAC, daily resolution, 1e-10 cutoff',bool)

yearf=(1990:(1/360/24):2004)';
%yearf=(1850:(1/360/24):2004)';
num=yf2num(yearf);
yearf2=yearfrac(num);
bool=maxmax(abs(yearf-yearf2))<1e-10;
reporttest('YF2NUM and YEARFRAC, hourly resolution, 1e-10 cutoff',bool)

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [x_sil, y_sil] = removeSilentFrames(x, y, range, N, K)

x       = x(:);
y       = y(:);

frames  = 1:K:(length(x)-N);
w       = ml_hanning(N);
msk     = zeros(size(frames));

for j = 1:length(frames)
    jj      = frames(j):(frames(j)+N-1);
    msk(j) 	= 20*log10(norm(x(jj).*w)./sqrt(N));
end

msk     = (msk-max(msk)+range)>0;
count   = 1;

x_sil   = zeros(size(x));
y_sil   = zeros(size(y));

for j = 1:length(frames)
    if msk(j)
        jj_i            = frames(j):(frames(j)+N-1);
        jj_o            = frames(count):(frames(count)+N-1);
        x_sil(jj_o)     = x_sil(jj_o) + x(jj_i).*w;
        y_sil(jj_o)  	= y_sil(jj_o) + y(jj_i).*w;
        count           = count+1;
    end
end

x_sil = x_sil(1:jj_o(end));
y_sil = y_sil(1:jj_o(end));

author: mpariente

repo: pystoi

save_path: github-repos/MATLAB/mpariente-pystoi


text:
function [ pp, normal ] = plane_exp2normal_3d ( p1, p2, p3 )

%*****************************************************************************80
%
%% PLANE_EXP2NORMAL_3D converts an explicit plane to normal form in 3D.
%
%  Discussion:
%
%    The explicit form of a plane in 3D is
%
%      the plane through P1, P2 and P3.
%
%    The normal form of a plane in 3D is
%
%      PP, a point on the plane, and
%      N, the unit normal to the plane.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3), P2(3), P3(3), three points on the plane.
%
%    Output, real PP(3), a point on the plane.
%
%    Output, real NORMAL(3), a unit normal vector to the plane.
%
  dim_num = 3;

  pp(1:dim_num) = p1(1:dim_num);

  normal(1) = ( p2(2) - p1(2) ) * ( p3(3) - p1(3) ) ...
            - ( p2(3) - p1(3) ) * ( p3(2) - p1(2) );

  normal(2) = ( p2(3) - p1(3) ) * ( p3(1) - p1(1) ) ...
            - ( p2(1) - p1(1) ) * ( p3(3) - p1(3) );

  normal(3) = ( p2(1) - p1(1) ) * ( p3(2) - p1(2) ) ...
            - ( p2(2) - p1(2) ) * ( p3(1) - p1(1) );

  norm = sqrt ( sum ( normal(1:dim_num).^2 ) );

  if ( norm == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_EXP2NORMAL_3D - Fatal error!\n' );
    fprintf ( 1, '  The normal vector is null.\n' );
    fprintf ( 1, '  Two points coincide, or nearly so.\n' );
    error ( 'PLANE_EXP2NORMAL_3D - Fatal error!' );
  end

  normal(1:dim_num) = normal(1:dim_num) / norm;

  return
end


text:
function  test_sum_quadratic()

    clc;
    clear;
    close all;

    %% Set algorithms
    if 0
        algorithms = sgd_solver_list('ALL');
    else
        algorithms = {'SGD','SVRG','IQN'};
    end

    % prepare datasets
    N = 1000;   % number of random functions
    d = 100;    % dimension

    A = zeros(d,d,N);
    b = zeros(d,N);

    for i=1:N
        for j=1:(d/2)
           %A(j,j,i)=10^(0+floor(3*rand(1)));  % good Condition number
           A(j,j,i)=10^(-2+floor(3*rand(1))); %  bad Condition number
        end
        for j=((d/2)+1):d
            A(j,j,i)=10^(2+floor(3*rand(1)));
        end

    end

    for i=1:N
        for j=1:d
            b(j,i)=1000*rand(1);
        end
    end

    A_sum=zeros(d,d);
    for i=1:N
        A_sum=A_sum+A(:,:,i);
    end

    b_sum=zeros(d,1);
    for i=1:N
        b_sum=b_sum+b(:,i);
    end

    cn = max(eig(A_sum))/min(eig(A_sum));
    fprintf('condition number: %e\n', cn);

    % define problem definitions
    problem = sum_quadratic(A, b);

    % Calculate the solution
    A_inv=zeros(d,d);
    for i=1:d
        A_inv(i,i)=1/(A_sum(i,i));
    end

    w_opt=-A_inv*b_sum;
    f_opt = problem.cost(w_opt);
    fprintf('%f\n', f_opt);

    %% initialize
    w_init = zeros(d,1);
    batch_size = 1;
    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_optgap = 1e-12;
        options.max_iter = 100;
        options.verbose = true;
        options.f_opt = f_opt;
        options.store_w = false;
        options.permute_on = 1;
        options.f_opt = f_opt;
        options.batch_size = batch_size;

        switch algorithms{alg_idx}

           case {'SGD'}

                options.step_init = 0.0001 * options.batch_size;
                %options.step_alg = 'decay';
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = sgd(problem, options);

           case {'IQN'}

                options.step_init = 1;
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = iqn(problem, options);

            case {'SVRG'}

                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = svrg(problem, options);

            case {'SAG'}

                options.step_init = 0.000001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SAG';

                [w_list{alg_idx}, info_list{alg_idx}] = sag(problem, options);

            case {'SAGA'}

                options.step_init = 0.00000001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SAGA';

                [w_list{alg_idx}, info_list{alg_idx}] = sag(problem, options);

            case {'SVRG-SQN'}

                options.batch_hess_size = batch_size * 20;
                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SVRG-SQN';
                options.L = 20;
                options.r = 20;

                [w_list{alg_idx}, info_list{alg_idx}] = slbfgs(problem, options);

            case {'SVRG-LBFGS'}

                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SVRG-LBFGS';
                options.mem_size = 20;

                [w_list{alg_idx}, info_list{alg_idx}] = slbfgs(problem, options);

            case {'Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    %% plot all
    close all;

    % display iter vs cost/gnorm
    display_graph('iter','optimality_gap', algorithms, w_list, info_list);
    display_graph('time','optimality_gap', algorithms, w_list, info_list);
    display_graph('iter','gnorm', algorithms, w_list, info_list);

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
classdef prtUtilIntegerAssociativeArrayClassNames < prtUtilIntegerAssociativeArray
    % prtUtilIntegerAssociativeArray
    % xxx Need Help xxx

    methods
        function self = merge(self,other)
            temp = merge@prtUtilIntegerAssociativeArray(self,other);
            self = prtUtilIntegerAssociativeArrayClassNames(temp.integerKeys,temp.cellValues);
        end

        function obj = prtUtilIntegerAssociativeArrayClassNames(integerKeys,cellValues)

            if nargin == 0
                return;
            else
                obj.integerKeys = integerKeys;
                obj.cellValues = cellValues;
            end
        end

        function value = get(obj,key)
            value = cell(size(key));
            uKeys = unique(key);

            %this should be significantly faster for large vectors of keys,
            %which, btw, you should use instead of iterating over individual keys
            for uKeyInd = 1:length(uKeys)
                if isnan(uKeys(uKeyInd))
                    value(isnan(key)) = {'unlabeled'};
                else
                    keyInds = find(key == uKeys(uKeyInd));
                    refInds = find(uKeys(uKeyInd) == obj.integerKeys,1);

                    if isempty(refInds)
                        value(keyInds) = {sprintf('Class %d',uKeys(uKeyInd))};
                    else
                        value(keyInds) = {obj.cellValues{refInds}};
                    end
                end
            end

            %             for i = 1:length(key)
            %                 ind = find(key(i) == obj.integerKeys, 1);
            %                 if isempty(ind)
            %                     value{i} = sprintf('Class %d',key(i));
            %                 else
            %                     value{i} = obj.cellValues{ind};
            %                 end
            %             end
            %Don't return a cell in this case
%             if length(key) == 1
%                 value = value{1};
%             end

            % This has been changed 2013-06-12 because it is weird to do
            % this.
        end
    end
end


text:
% mff_getobj - import information from MFF object
%
% Usage:
%   values = mff_setobj(mffObj, structure);
%
% Inputs:
%  mffObj    - MFF object of any type
%  structure - description of variables to assign
%
% Outputs:
%  values    - values for the structure
%
% Note: this function may take any MFF object structure

% This file is part of mffmatlabio.
%
% mffmatlabio 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.
%
% mffmatlabio 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 mffmatlabio.  If not, see <https://www.gnu.org/licenses/>.

function matStruct = mff_getobj(mffObj, variables)

for iVar = 1:size(variables,1)
    varName  = variables{iVar,1};
    varType  = variables{iVar,2};
    varArray = variables{iVar,3};

    if ~isequal(varType, 'array')
        eval( [ 'matStruct.(lower(varName)) = ' varType '(mffObj.get' varName '());' ] );
    else
        if isempty(varName)
            for iArray = 1:mffObj.size
                matStruct.value(iArray) = mffObj.get(iArray-1);
            end
        else
            eval( [ 'tmparray = mffObj.get' varName '();' ] );
            if ~isempty(tmparray)
                for iArray = 1:tmparray.size
                    tmpval = tmparray.get(iArray-1);
                    if ~isempty(varArray)
                        matStruct.(lower(varName))(iArray) = mff_getobj(tmpval, varArray);
                    else
                        matStruct.(lower(varName))(iArray) = tmpval;
                    end
                end
            else
                matStruct.(lower(varName)) = [];
            end
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function rectangle_node = bezier_surface_rectangle_read ( rectangle_file_name, ...
  rectangle_num )

%*****************************************************************************80
%
%% BEZIER_SURFACE_RECTANGLE_READ reads rectangles from a rectangle file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string RECTANGLE_FILE_NAME, the name of the
%    rectangle file.
%
%    Input, integer RECTANGLE_NUM, the number of rectangles.
%
%    Output, real RECTANGLE_NODE(16,RECTANGLE_NUM),
%    the nodes that make up each rectangle.
%
  dim_num = 16;

  rectangle_node = i4mat_data_read ( rectangle_file_name, dim_num, ...
    rectangle_num );

  return
end


text:
% :Usage:
% ::
%
%    [c, alld] = getVertexColors(xyz, v, actcolor, [basecolor], [mind], 'vert', [xyz2], [actcolor2], 'vert', [xyz3], [actcolor3])
%
% given a point list of XYZ mm coordinates (3 columns)
% and a list of vertices in an isosurface,
% returns FaceVertexCData color values for brain near points and brain not near points.
% c is vertex color specification, 3 columns indicating RGB values
%
% Inputs:
% xyz   a 3-vol list of vertices to color
% v     can be a matrix of vertices
%       or a handle to a patch object containing vertices
%       if it's a handle, this function sets the color to interp
%       and the FaceVertexCData to the color matrix c
% actcolor
%       [r g b] activation color
% basecolor
%       [r g b] baseline color - optional.
% mind  optional - min distance to color vertex
%       Vertices within mind of an xyz coordinate will be colored
% colorscale
%       optional.  followed by vector of values by which to multiply input
%       color
%       these are scaled to be between .3 and one.
%       if entered, this will make the colors vary by, for example, Z score
%       so Z-scores are an acceptable input.
%       cscale should be in the same coordinate order as xyz
%       for ADDITIONAL clusters, repeat the 'colorscale', Z argument pair in the function call
%
%       YOU CAN ALSO pass true RGB values for each xyz coordinate in: 'colorscale', rgblist,
%       IF cscale is a 3-vector, it specifies the ACTUAL colors, and is not scaled to .3 - 1
%
%
% following basecolor and mind:
% additional xyz coordinate lists, with syntax:
% 'vert', xyz2 [your xyz input], [r g b] color for xyz plot
%
% also, you can enter 'ovlcolor' followed by [r g b] for overlaps between xyz sets
%   colors will ONLY appear in the overlap color if they share actual coordinates in common,
%   not necessarily if surface vertices are within the specified distance from both sets of coords.
%
% to get a good brain surface, try this:
%figure
%p = patch('Faces', faces, 'Vertices', vertices, 'FaceColor', [.5 .5 .5], ...
% 'EdgeColor', 'none', 'SpecularStrength', .2, 'FaceAlpha', 1, 'SpecularExponent', 200);
%lighting gouraud;camlight right
%axis image; myLight = camlight(0, 0);set(myLight, 'Tag', 'myLight');
%set(gcf, 'WindowButtonUpFcn', 'lightFollowView');lightfollowview
%drawnow
%
% ..
% by Tor Wager  August 25, 2002
% ..

function [c, alld] = getVertexColors_old_backup(xyz, v, actcolor, varargin)

    mind = 3;
    basecolor = [.5 .5 .5];
    alld = [];
    xyza = xyz;
    cscale = [];
    allda = [];
    vv = [];

    % -----------------------------------------------------------------------
    % * set up input arguments
    % -----------------------------------------------------------------------
    doalph = 0; cscale = [];
    if ~isempty(varargin), basecolor = varargin{1}; end
    if length(varargin) > 1, mind = varargin{2}; end
    ind = 1;
    for i = 3:length(varargin)
        if strcmp(varargin{i}, 'vert')
            vv{ind} = varargin{i+1};

            % intersections
            xyzb{ind, 1} = intersect(xyz, vv{ind}, 'rows');
            for j = ind-1:-1:1
                xyzb{ind, j} = intersect(vv{j}, vv{ind}, 'rows');
                xyza = intersect(xyza, xyzb{ind, j}, 'rows');
            end

            cc{ind} = varargin{i+2};
            ind = ind+1;
        elseif strcmp(varargin{i}, 'ovlcolor')
            ocol = varargin{i+1};
        elseif strcmp(varargin{i}, 'alphaone')
            doalph = 1;
        elseif strcmp(varargin{i}, 'allcolor')
            acol = varargin{i+1};
        elseif strcmp(varargin{i}, 'colorscale')
            cscale{end+1} = varargin{i+1};
            if min(size(cscale{end})) == 1
                % scale colors (may be necessary) - only if single vector, not RGB values
                cscale{end} = cscale{end} ./ max(cscale{end});
            end
            if any(cscale{end} < 0), error('Some color scale values are less than zero.'), end

        end
    end

    if isempty(cscale)
        cscale{1} = ones(size(xyz, 1), 1);
        for i = 1:length(vv)    % additional vertices
            cscale{end+1} = ones(size(vv{i}, 1), 1);
        end
    end

    %if ~isempty(cscale) don't need this.
    %    if length(cscale) < length(vv)
    %        cscale{length(vv)} = [];
    %    end
    %end

    if ishandle(v)
        p = v;
        v = get(p, 'Vertices');
        c = get(p, 'FaceVertexCData');   % get existing colors from surface
    else
        % uh-oh, not a handle!
        warning('Figure handle missing: Figure was closed?')
        p = findobj('Type', 'patch');
        v = get(p(1), 'Vertices');
        c = get(p, 'FaceVertexCData');   % get existing colors from surface
    end

    bad = any(size(c) - size(v));
    if bad
        c = repmat(basecolor, size(v, 1), 1);
    end

    % -----------------------------------------------------------------------
    % * main xyz color change
    % -----------------------------------------------------------------------

    t1 = clock;
    fprintf('Main color vertices: ')

    c = change_colors(c, xyz, v, mind, cscale, actcolor, p);
    drawnow();

    % -----------------------------------------------------------------------
    % * additional optional vertices
    % -----------------------------------------------------------------------
    if exist('vv', 'var')
        for j = 1:length(vv)
            fprintf('\nAdditional vertices: ')

            % cscale:
            % pass in vector of ones length coords for solid color
            % or scalar vals for color mapping
            % pass in cell array

            % cc{j} is color, cscale{j+1} is scaling vals for coords
            c = change_colors(c, vv{j}, v, mind, cscale(j+1), cc{j}, p);
        end
        drawnow();
    end

    % -----------------------------------------------------------------------
    % * figure out which vertices should be colored with ocol (overlap color)
    % -----------------------------------------------------------------------

    if exist('ocol', 'var') && exist('xyzb', 'var') && ~isempty(cat(1, xyzb{:}))

        alld = zeros(size(v, 1), 1);  % keeps track of overlap vertices
        xyzb = cat(1, xyzb{:});

        t1 = clock;
        fprintf('\nOverlap vertices: ')

        cscaletmp = {ones(size(xyzb, 1), 1)};
        c = change_colors(c, xyzb, v, mind, cscaletmp, ocol, p);
        drawnow();

        fprintf('%3.0f.done in %3.0f s\n', i, etime(clock, t1))
    end

    % -----------------------------------------------------------------------
    % * figure out which vertices should be colored with acol (all color)
    % -----------------------------------------------------------------------
    if exist('acol', 'var') && ~isempty(xyza) && length(vv)>1

        xyzall = xyza;
        for i = 2:length(vv)
            xyzall = intersect(xyzall, vv{i}, 'rows');
        end

        t1 = clock;
        fprintf('\nAll overlap vertices: ')

        cscaletmp = {ones(size(xyzall, 1), 1)};
        c = change_colors(c, xyzall, v, mind, cscaletmp, acol, p);
        drawnow();

        fprintf('%3.0f.done in %3.0f s\n', i, etime(clock, t1))

    end

    % -----------------------------------------------------------------------
    % * final color change
    % -----------------------------------------------------------------------

    if exist('p', 'var')
        set(p, 'FaceColor', 'interp')
        set(p, 'FaceVertexCData', c)
        drawnow()
    end

    %lightFollowView
end

% -----------------------------------------------------------------------
% * SUB-FUNCTIONS
% -----------------------------------------------------------------------

function c = change_colors(c, coords, v, mind, cscale, actcolor, p)

    if isempty(coords), disp('Coords is empty. Nothing to plot.'), return, end

    % select vertices that are even close
    cmax = max(coords, [], 1);
    cmin = min(coords, [], 1);
    fprintf('%3.0f vertices.  selecting: ', size(v, 1));
    wh = any(v - repmat(cmax, size(v, 1), 1) > mind, 2);
    wh2 = any(repmat(cmin, size(v, 1), 1) - v > mind, 2);

    % list vertices to test and possibly change color
    whverts = (1:size(v, 1))';
    whverts(wh | wh2) = [];     % vertex indices in big list
    smallv = v;
    smallv(wh | wh2,:) = [];     % vertices--restricted list

    fprintf('%3.0f\n', size(whverts, 1));

    if isempty(smallv), return, end

    % select coords that are even close
    % ----------------------------------
    cmax = max(smallv, [], 1);
    cmin = min(smallv, [], 1);
    fprintf('%3.0f coords.  selecting: ', size(coords, 1));
    wh = any(coords - repmat(cmax, size(coords, 1), 1) > mind, 2);
    wh2 = any(repmat(cmin, size(coords, 1), 1) - coords > mind, 2);

    % if cscale is matrix, must select these values of cscale as well!
    if size(cscale{1}, 1) == size(coords, 1)
        cscale{1}(wh | wh2,:) = [];
    end

    coords(wh | wh2,:) = [];

    if isempty(coords), return, end

    nc = size(coords, 1);
    fprintf('%3.0f\n', nc);

    % break up coords into list and run
    % ----------------------------------

    % break up coords into list
    xyz2 = {}; indx = 1;
    for kk = 1:1000:nc
        setwh{indx} = (kk:min(nc, kk + 1000 - 1))';
        xyz2{indx} = coords(setwh{indx},:);

        indx = indx + 1;
    end

    fprintf('Running %3.0f sets of coordinates: 000', length(xyz2));

    indxval = 1;
    wh_coords_near_surface = false(size(coords, 1), 1);

    for setno = 1:length(xyz2)
        fprintf('\b\b\b%03d', setno);

        for i = 1:size(xyz2{setno}, 1)
            % find vertices that are within range of point i in set  setno
            vertex_indices = find_in_radius(xyz2, setno, i, smallv, mind, whverts);

            % two modes: if cscale{1} is a matrix, treats as rgb values, and put in
            % color stored in cscale.  if cscale{1} is a vector, treat it as a scaling value for actcolor
            % In either case, indxval should index location of coordinate in FULL
            % list (corresponding to full list in cscale)

            c = color_change_vertices(c, mind, indxval, cscale, actcolor, vertex_indices);

            if ~isempty(vertex_indices), wh_coords_near_surface(indxval) = 1; end

            indxval = indxval + 1;
        end

        if exist('p', 'var')
            set(p, 'FaceColor', 'interp')
            set(p, 'FaceVertexCData', c)
        end
    end
end

function z = dist_tmp(w, p)

    %
    % if isstr(w)
    %   switch (w)
    %     case 'deriv',
    %       z = '';
    %     otherwise
    %       error('Unrecognized code.')
    %   end
    %   return
    % end

    % CALCULATION
    if nargin == 1
        p = w;
        w = w';
    end

    [S, R] = size(w);
    [R2, Q] = size(p);
    if (R ~= R2), error('Inner matrix dimensions do not match.'), end

    z = zeros(S, Q);
    if (Q<S)
        p = p';
        copies = zeros(1, S);
        for q=1:Q
            z(:,q) = sum((w-p(q+copies,:)).^2, 2);
        end
    else
        w = w';
        copies = zeros(1, Q);
        for i=1:S
            z(i,:) = sum((w(:,i+copies)-p).^2, 1);
        end
    end
    z = sqrt(z);

end

function [vertex_indices, d] = find_in_radius(xyz2, setno, i, smallv, mind, whverts)
    % output: indices of vertices in BIG list, and distances

    % get vertices v within box -- fast method
    wh = find(all(abs(bsxfun(@minus, xyz2{setno}(i,:), smallv)) <= mind, 2));
    d = dist_tmp(smallv(wh,:), xyz2{setno}(i,:)');

    % convert back to big list
    vertex_indices = whverts(wh(d < mind));
end

function c = color_change_vertices(c, mind, i, cscale, actcolor, vertex_indices)

    % two modes: if cscale{1} is a matrix, treats as rgb values, and put in
    % color stored in cscale
    % if cscale{1} is a vector, treat it as a scaling value for actcolor
    % In either case, the ith index point should correspond to the
    % coordinate being worked on.

    n = length(vertex_indices);

    if n
        if min(size(cscale{1})) > 1  % if we have rgb values rather than scaling values
            % this occurs if heatmap = yes and colorscale = no, we pass in rgb values
            c(vertex_indices,:) = repmat(cscale{1}(i,:), n, 1);
        else
            c(vertex_indices,:) = repmat(actcolor.*cscale{1}(i,:), n, 1);
        end
    end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function r = isempty(s)

%tstoolbox/@core/isempty
%   Syntax:
%     * r = isempty(s)
%
%   Input Arguments:
%     * s - core object
%
%   test if core contains no (valid) data
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

r = isempty(s.data);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa