text:
function [out_vec] = tgn_model(in_vec)
% tgn_model.m - TGn channel model simulation
global numPkts PER_snr snc_idxs snc_dlys snr_idx H_iid H_dline P_tap Kf_tap Kv_tap R_tx R_rx;

% Channel Model Test Setup
%    0 - AWGN channel
%    1 - Channel D (no Doppler)
%    2 - Channel D (with Doppler)
test_setup = 1;

% PER graph settings (actual values are user-selected (GUI))
%%%%% numPkts = 200; %100
%%%%% PER_snr = 29:33; %15:1:21;

% Input parameters
Clk = in_vec(1);
Ntx = in_vec(2);
Nrx = in_vec(3);

% Channel inputs
in_len = length(in_vec(4:end))/4;
ch_in  = zeros(4, in_len);
ch_in(1,:) = in_vec(0*in_len+3+(1:in_len)).';
ch_in(2,:) = in_vec(1*in_len+3+(1:in_len)).';
ch_in(3,:) = in_vec(2*in_len+3+(1:in_len)).';
ch_in(4,:) = in_vec(3*in_len+3+(1:in_len)).';
testch_out = ch_in;
ch_in = ch_in(1:Ntx,:);  % keep only used Tx ant. inputs

% Channel outputs
ch_out = zeros(4, in_len);

if (test_setup==0) %%%%%%%%%%%%%%%  AWGN channel

    if (Clk==0)  snr_idx = 1;  end

    % Add AWGN to Rx streams (reduce by 5dB, due to upsampling)
    snr_val = PER_snr(snr_idx);
    for i=1:Nrx
        testch_out(i,:) = awgn(testch_out(i,:), snr_val-5, 'measured');
    end

    % Next SNR Value for PER Curve...
    if (mod(Clk,numPkts)==(numPkts-1) && snr_idx<length(PER_snr))
        snr_idx = snr_idx+1;
    end

    % Return output
    out_vec = reshape(testch_out.', 4*in_len, 1);

else  %%%%%%%%%%%%%%%  TGn Channel Model D

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% For Channel D (NLOS) %%%%%
chD_K = -inf;  %%% not used (since NLOS), K = -inf dB
chD_dly    = [    0     10     20     30     40     50     60     70     80     90    110    140    170    200    240    290    340    390];
chD_wdth   = [   10     10     10     10     10     10     10     10     10     20     30     30     30     40     50     50     50     50];
% Cluster 1 - tap powers, delays
chD_pow_c1 = [    0   -0.9   -1.7   -2.6   -3.5   -4.3   -5.2   -6.1   -6.9   -7.8   -9.0  -11.1  -13.7  -16.3  -19.3  -23.2   -inf   -inf];
chD_aoa_c1 = [158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9  158.9];
chD_asa_c1 = [ 27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7   27.7];
chD_aod_c1 = [332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1  332.1];
chD_asd_c1 = [ 27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4   27.4];
% Cluster 2 - tap powers, delays
chD_pow_c2 = [ -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -6.6   -9.5  -12.1  -14.7  -17.4  -21.9  -25.5   -inf];
chD_aoa_c2 = [320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2  320.2];
chD_asa_c2 = [ 31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4   31.4];
chD_aod_c2 = [ 49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3   49.3];
chD_asd_c2 = [ 32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1   32.1];
% Cluster 3 - tap powers, delays
chD_pow_c3 = [ -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf   -inf  -18.8  -23.2  -25.2  -26.7];
chD_aoa_c3 = [276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1  276.1];
chD_asa_c3 = [ 37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4   37.4];
chD_aod_c3 = [275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9  275.9];
chD_asd_c3 = [ 36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8   36.8];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num_taps = length(chD_dly);

% define imag. component
j = sqrt(-1);

% Initialize channel model
if (Clk==0)

    snr_idx = 1;

    % Initialize random generator
    %12345
    %randn('state',12345);

    % initialize variables
    snc_idxs = zeros(num_taps, 12, 16);
    snc_dlys = zeros(num_taps, 12, 16);
    P_tap  = zeros(num_taps, 1);
    Kf_tap = zeros(num_taps, 1);
    Kv_tap = zeros(num_taps, 1);
    R_tx   = zeros(Ntx, Ntx, num_taps);
    R_rx   = zeros(Nrx, Nrx, num_taps);

    % Initialize channel delay line
    tot_dly = sum(chD_wdth)/10;
    H_dline = zeros(Ntx, tot_dly);

    % Initialize calculations for each tap
    for m=1:num_taps

        % initialize Doppler-filtered, Rayleigh fading for tap
        for y=1:16 % max size = 4x4
            s_idxs = 0;
            s_dlys = 0;
            [out_val, s_dlys, s_idxs] = sinc_dop_rv(s_dlys, s_idxs, 256);
            snc_idxs(m,:,y) = s_idxs;
            snc_dlys(m,:,y) = s_dlys;
        end

        % Compute overall tap power
        P_tap(m) = 10^(chD_pow_c1(m)/10) + 10^(chD_pow_c2(m)/10) + 10^(chD_pow_c3(m)/10);
        P_tap(m) = sqrt(P_tap(m));
        % Compute K-factor for tap (convert dB to linear)
        if (m==1)  K = 10^(chD_K/10);
        else       K = 0;              end
        Kf_tap(m) = sqrt(K/(K+1));
        Kv_tap(m) = sqrt(1/(K+1));
        % Compute correlation matrices
        sq2 = sqrt(2);
        ds=pi;    % antenna spacing
        Dc = 4;   % matrix index for distance zero
        ci=1:3;   % cluster index vector
        mxx=1:20;  % 2nd sum index vector, Rxx
        mxy=0:20;  % 2nd sum index vector, Rxy
        dph=60*(pi/180); % delta phi = pi (-180 to +180 degrees)
        Q = [ 10^(chD_pow_c1(m)/10) 10^(chD_pow_c2(m)/10) 10^(chD_pow_c3(m)/10) ];
        otx  = (pi/180)*[ chD_asd_c1(m) chD_asd_c2(m) chD_asd_c3(m) ];
        orx  = (pi/180)*[ chD_asa_c1(m) chD_asa_c2(m) chD_asa_c3(m) ];
        phtx = (pi/180)*[ chD_aod_c1(m) chD_aod_c2(m) chD_aod_c3(m) ];
        phrx = (pi/180)*[ chD_aoa_c1(m) chD_aoa_c2(m) chD_aoa_c3(m) ];
        % Tx correlations...
        normQtx = sum(Q.*(1-exp(-sq2*dph./otx)));   % norm. constants
        Qtx = Q/normQtx;
        for D = -3:3
            Rxx_c1 = sum( besselj(2*mxx,D*ds)./((sq2/otx(1)).^2+(2*mxx).^2).*cos(2*mxx*phtx(1)).*...
                          (sq2/otx(1) + exp(-dph*sq2/otx(1)).*(2*mxx.*sin(2*mxx*dph) - (sq2/otx(1))*cos(2*mxx*dph))) );
            Rxx_c2 = sum( besselj(2*mxx,D*ds)./((sq2/otx(2)).^2+(2*mxx).^2).*cos(2*mxx*phtx(2)).*...
                          (sq2/otx(2) + exp(-dph*sq2/otx(2)).*(2*mxx.*sin(2*mxx*dph) - (sq2/otx(2))*cos(2*mxx*dph))) );
            Rxx_c3 = sum( besselj(2*mxx,D*ds)./((sq2/otx(3)).^2+(2*mxx).^2).*cos(2*mxx*phtx(3)).*...
                          (sq2/otx(3) + exp(-dph*sq2/otx(3)).*(2*mxx.*sin(2*mxx*dph) - (sq2/otx(3))*cos(2*mxx*dph))) );
            Rxx_c = [Rxx_c1 Rxx_c2 Rxx_c3];
            Rxx(Dc+D) = besselj(0,D*ds)+4*sum(Qtx./(otx.*sq2).*Rxx_c);
            % cross-correlation, real and imag
            Rxy_c1 = sum( besselj(2*mxy+1,D*ds)./((sq2/otx(1)).^2+(2*mxy+1).^2).*sin((2*mxy+1)*phtx(1)).*...
                          (sq2/otx(1) - exp(-dph*sq2/otx(1)).*((2*mxy+1).*sin((2*mxy+1)*dph) + (sq2/otx(1))*cos((2*mxy+1)*dph))) );
            Rxy_c2 = sum( besselj(2*mxy+1,D*ds)./((sq2/otx(2)).^2+(2*mxy+1).^2).*sin((2*mxy+1)*phtx(2)).*...
                          (sq2/otx(2) - exp(-dph*sq2/otx(2)).*((2*mxy+1).*sin((2*mxy+1)*dph) + (sq2/otx(2))*cos((2*mxy+1)*dph))) );
            Rxy_c3 = sum( besselj(2*mxy+1,D*ds)./((sq2/otx(3)).^2+(2*mxy+1).^2).*sin((2*mxy+1)*phtx(3)).*...
                          (sq2/otx(3) - exp(-dph*sq2/otx(3)).*((2*mxy+1).*sin((2*mxy+1)*dph) + (sq2/otx(3))*cos((2*mxy+1)*dph))) );
            Rxy_c = [Rxy_c1 Rxy_c2 Rxy_c3];
            Rxy(Dc+D) = 4*sum(Qtx./(otx.*sq2).*Rxy_c);
        end
        % Form Tx corr. matrix
        Rtx = [ Rxx(Dc)   Rxx(Dc+1) Rxx(Dc+2) Rxx(Dc+3); ...
                Rxx(Dc-1) Rxx(Dc)   Rxx(Dc+1) Rxx(Dc+2); ...
                Rxx(Dc-2) Rxx(Dc-1) Rxx(Dc)   Rxx(Dc+1); ...
                Rxx(Dc-3) Rxx(Dc-2) Rxx(Dc-1) Rxx(Dc)    ];
        Rtx = Rtx + j*...
              [ Rxy(Dc)   Rxy(Dc+1) Rxy(Dc+2) Rxy(Dc+3); ...
                Rxy(Dc-1) Rxy(Dc)   Rxy(Dc+1) Rxy(Dc+2); ...
                Rxy(Dc-2) Rxy(Dc-1) Rxy(Dc)   Rxy(Dc+1); ...
                Rxy(Dc-3) Rxy(Dc-2) Rxy(Dc-1) Rxy(Dc)    ];
        Rtx = Rtx(1:Ntx,1:Ntx); % keep only Ntx-by-Ntx elements
        R_tx(:,:,m) = Rtx^(1/2);

        % Rx correlations...
        normQrx = sum(Q.*(1-exp(-sq2*dph./orx)));   % norm. constants
        Qrx = Q/normQrx;
        for D = -3:3
            % cross-correlation, real parts
            Rxx_c1 = sum( besselj(2*mxx,D*ds)./((sq2/orx(1)).^2+(2*mxx).^2).*cos(2*mxx*phrx(1)).*...
                          (sq2/orx(1) + exp(-dph*sq2/orx(1)).*(2*mxx.*sin(2*mxx*dph) - (sq2/orx(1))*cos(2*mxx*dph))) );
            Rxx_c2 = sum( besselj(2*mxx,D*ds)./((sq2/orx(2)).^2+(2*mxx).^2).*cos(2*mxx*phrx(2)).*...
                          (sq2/orx(2) + exp(-dph*sq2/orx(2)).*(2*mxx.*sin(2*mxx*dph) - (sq2/orx(2))*cos(2*mxx*dph))) );
            Rxx_c3 = sum( besselj(2*mxx,D*ds)./((sq2/orx(3)).^2+(2*mxx).^2).*cos(2*mxx*phrx(3)).*...
                          (sq2/orx(3) + exp(-dph*sq2/orx(3)).*(2*mxx.*sin(2*mxx*dph) - (sq2/orx(3))*cos(2*mxx*dph))) );
            Rxx_c = [Rxx_c1 Rxx_c2 Rxx_c3];
            Rxx(Dc+D) = besselj(0,D*ds)+4*sum(Qrx./(orx.*sq2).*Rxx_c);
            % cross-correlation, real and imag
            Rxy_c1 = sum( besselj(2*mxy+1,D*ds)./((sq2/orx(1)).^2+(2*mxy+1).^2).*sin((2*mxy+1)*phrx(1)).*...
                          (sq2/orx(1) - exp(-dph*sq2/orx(1)).*((2*mxy+1).*sin((2*mxy+1)*dph) + (sq2/orx(1))*cos((2*mxy+1)*dph))) );
            Rxy_c2 = sum( besselj(2*mxy+1,D*ds)./((sq2/orx(2)).^2+(2*mxy+1).^2).*sin((2*mxy+1)*phrx(2)).*...
                          (sq2/orx(2) - exp(-dph*sq2/orx(2)).*((2*mxy+1).*sin((2*mxy+1)*dph) + (sq2/orx(2))*cos((2*mxy+1)*dph))) );
            Rxy_c3 = sum( besselj(2*mxy+1,D*ds)./((sq2/orx(3)).^2+(2*mxy+1).^2).*sin((2*mxy+1)*phrx(3)).*...
                          (sq2/orx(3) - exp(-dph*sq2/orx(3)).*((2*mxy+1).*sin((2*mxy+1)*dph) + (sq2/orx(3))*cos((2*mxy+1)*dph))) );
            Rxy_c = [Rxy_c1 Rxy_c2 Rxy_c3];
            Rxy(Dc+D) = 4*sum(Qrx./(orx.*sq2).*Rxy_c);
        end
        % Form Rx corr. matrix
        Rrx = [ Rxx(Dc)   Rxx(Dc+1) Rxx(Dc+2) Rxx(Dc+3); ...
                Rxx(Dc-1) Rxx(Dc)   Rxx(Dc+1) Rxx(Dc+2); ...
                Rxx(Dc-2) Rxx(Dc-1) Rxx(Dc)   Rxx(Dc+1); ...
                Rxx(Dc-3) Rxx(Dc-2) Rxx(Dc-1) Rxx(Dc)    ];
        Rrx = Rrx + j*...
              [ Rxy(Dc)   Rxy(Dc+1) Rxy(Dc+2) Rxy(Dc+3); ...
                Rxy(Dc-1) Rxy(Dc)   Rxy(Dc+1) Rxy(Dc+2); ...
                Rxy(Dc-2) Rxy(Dc-1) Rxy(Dc)   Rxy(Dc+1); ...
                Rxy(Dc-3) Rxy(Dc-2) Rxy(Dc-1) Rxy(Dc)    ];
        Rrx = Rrx(1:Nrx,1:Nrx); % keep only Nrx-by-Nrx elements
        R_rx(:,:,m) = (Rrx^(1/2)).';
    end
end

% Without Doppler, generate channel matrices
if (test_setup==1)
    H_tap = zeros(Nrx,Ntx,num_taps);
    for m=1:num_taps
        % Gaussian RV's for Rayleigh fading
        H_iid = randn(Nrx, Ntx) + j*randn(Nrx, Ntx);
        H_iid = H_iid/sqrt(2);

        % Fixed, LOS matrix
        H_f = ones(Nrx, Ntx);   % LOS (since AoA=AoD=45)

        % Variable, NLOS matrix (Rayleigh matrix)
        H_v = R_rx(:,:,m)*H_iid*R_tx(:,:,m);

        % Compute channel tap
        H_tap(:,:,m) = P_tap(m)*( Kf_tap(m)*H_f + Kv_tap(m)*H_v );
    end
end

% Process all input samples...
for i=1:in_len

    % Add new sample to delay line
    H_dline = [ch_in(:,i) H_dline(:,1:end-1)];

    % Apply channel filter
    H_dline_idx = 1;
    H_out = zeros(Nrx,1);
    for m=1:num_taps

        % With Doppler, update 'H_tap' every sample
        if (test_setup==2)
            % Generate Doppler-filtered, Rayleigh fading
            H_iid = zeros(1, Nrx*Ntx);
            for y=1:Nrx*Ntx
                s_idxs = snc_idxs(m,:,y);
                s_dlys = snc_dlys(m,:,y);
                [out_val, s_dlys, s_idxs] = sinc_dop_rv(s_dlys, s_idxs, 256);
                snc_idxs(m,:,y) = s_idxs;
                snc_dlys(m,:,y) = s_dlys;
                H_iid(y) = out_val;
            end
            H_iid = reshape(H_iid, Nrx, Ntx);

            % Fixed, LOS matrix
            H_f = ones(Nrx, Ntx);   % LOS (since AoA=AoD=45)

            % Variable, NLOS matrix (Rayleigh matrix)
            H_v = R_rx(:,:,m)*H_iid*R_tx(:,:,m);

            % Compute channel tap
            H_tap(:,:,m) = P_tap(m)*( Kf_tap(m)*H_f + Kv_tap(m)*H_v );
        end

        % Apply tap to inputs (considering tap width)
        for k = 1:(chD_wdth(m)/10)
if 1 % 10 ns (100 MHz) sampling
            H_out = H_out + H_tap(:,:,m)*H_dline(:,H_dline_idx);
else % 50 ns (20 MHz) sampling
            if mod(H_dline_idx,5)==1
                H_out = H_out + H_tap(:,:,m)*H_dline(:,1+floor(H_dline_idx/5));
            end
end
            H_dline_idx = H_dline_idx + 1;
        end
    end

    % Store channel output
    ch_out(1:Nrx,i) = H_out;
end

% Add AWGN to Rx streams (reduce by 5dB, due to upsampling)
snr_val = PER_snr(snr_idx);
for i=1:Nrx
    ch_out(i,:) = awgn(ch_out(i,:), snr_val-5, 'measured');
end

% Next SNR Value for PER Curve...
if (mod(Clk,numPkts)==(numPkts-1) && snr_idx<length(PER_snr))
    snr_idx = snr_idx+1;
end

% Return output
out_vec = reshape(ch_out.', 4*in_len, 1);

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function invert = spm_cfg_eeg_inv_invert
% Configuration file for configuring imaging source inversion
% reconstruction
%__________________________________________________________________________
% Copyright (C) 2010-2016 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_inv_invert.m 7076 2017-05-19 12:47:36Z vladimir $

D = cfg_files;
D.tag = 'D';
D.name = 'M/EEG datasets';
D.filter = 'mat';
D.num = [1 Inf];
D.help = {'Select the M/EEG mat files.'};

val = cfg_entry;
val.tag = 'val';
val.name = 'Inversion index';
val.strtype = 'n';
val.help = {'Index of the cell in D.inv where the forward model can be found and the results will be stored.'};
val.val = {1};

all = cfg_const;
all.tag = 'all';
all.name = 'All';
all.val  = {1};
all.help = {''};

condlabel = cfg_entry;
condlabel.tag = 'condlabel';
condlabel.name = 'Condition label';
condlabel.strtype = 's';
condlabel.val = {''};
condlabel.help = {''};

conditions = cfg_repeat;
conditions.tag = 'conditions';
conditions.name = 'Conditions';
conditions.help = {'Specify the labels of the conditions to be included in the inversion'};
conditions.num  = [1 Inf];
conditions.values  = {condlabel};
conditions.val = {condlabel};

whatconditions = cfg_choice;
whatconditions.tag = 'whatconditions';
whatconditions.name = 'What conditions to include?';
whatconditions.values = {all, conditions};
whatconditions.val = {all};
whatconditions.help = {'What conditions to include?'};

standard = cfg_const;
standard.tag = 'standard';
standard.name = 'Standard';
standard.help = {'Use default settings for the inversion'};
standard.val  = {1};

invtype = cfg_menu;
invtype.tag = 'invtype';
invtype.name = 'Inversion type';
invtype.help = {'Select the desired inversion type'};
invtype.labels = {'MSP (GS)',  'COH', 'IID', 'EBB'};
invtype.values = {'GS',        'LOR', 'IID', 'EBB'};
invtype.val = {'GS'};

woi = cfg_entry;
woi.tag = 'woi';
woi.name = 'Time window of interest';
woi.strtype = 'r';
woi.num = [1 2];
woi.val = {[-Inf Inf]};
woi.help = {'Time window to include in the inversion (ms)'};

foi = cfg_entry;
foi.tag = 'foi';
foi.name = 'Frequency window of interest';
foi.strtype = 'r';
foi.num = [1 2];
foi.val = {[0 256]};
foi.help = {'Frequency window (the same as high-pass and low-pass in the GUI)'};

hanning = cfg_menu;
hanning.tag = 'hanning';
hanning.name = 'PST Hanning window';
hanning.help = {'Multiply the time series by a Hanning taper to emphasize the central part of the response.'};
hanning.labels = {'yes', 'no'};
hanning.values = {1, 0};
hanning.val = {1};

priorsmask  = cfg_files;
priorsmask.tag = 'priorsmask';
priorsmask.name = 'Priors file';
priorsmask.filter = '(.*\.gii$)|(.*\.mat$)|(.*\.nii(,\d+)?$)|(.*\.img(,\d+)?$)';
priorsmask.num = [0 1];
priorsmask.help = {'Select a mask or a mat file with priors.'};
priorsmask.val = {{''}};

space = cfg_menu;
space.tag = 'space';
space.name = 'Prior image space';
space.help = {'Space of the mask image.'};
space.labels = {'MNI', 'Native'};
space.values = {1, 0};
space.val = {1};

priors = cfg_branch;
priors.tag = 'priors';
priors.name = 'Source priors';
priors.help = {'Restrict solutions to pre-specified VOIs'};
priors.val  = {priorsmask, space};

locs  = cfg_entry;
locs.tag = 'locs';
locs.name = 'Source locations';
locs.strtype = 'r';
locs.num = [Inf 3];
locs.help = {'Input source locations as n x 3 matrix'};
locs.val = {zeros(0, 3)};

radius = cfg_entry;
radius.tag = 'radius';
radius.name = 'Radius of VOI (mm)';
radius.strtype = 'r';
radius.num = [1 1];
radius.val = {32};
radius.help = {''};

mask  = cfg_files;
mask.tag = 'mask';
mask.name = 'Mask image';
mask.filter = '(.*\.nii(,\d+)?$)|(.*\.img(,\d+)?$)';
mask.num = [0 1];
mask.help = {'Select a mask image'};
mask.val = {{''}};

restrict = cfg_branch;
restrict.tag = 'restrict';
restrict.name = 'Restrict solutions';
restrict.help = {'Restrict solutions to pre-specified VOIs'};
restrict.val  = {locs, radius, mask};

custom = cfg_branch;
custom.tag = 'custom';
custom.name = 'Custom';
custom.help = {'Define custom settings for the inversion'};
custom.val  = {invtype, woi, foi, hanning, priors, restrict};

isstandard = cfg_choice;
isstandard.tag = 'isstandard';
isstandard.name = 'Inversion parameters';
isstandard.help = {'Choose whether to use standard or custom inversion parameters.'};
isstandard.values = {standard, custom};
isstandard.val = {standard};

modality = cfg_menu;
modality.tag = 'modality';
modality.name = 'Select modalities';
modality.help = {'Select modalities for the inversion (only relevant for multimodal datasets).'};
modality.labels = {'All', 'EEG', 'MEG', 'MEGPLANAR', 'EEG+MEG', 'MEG+MEGPLANAR', 'EEG+MEGPLANAR', 'EEG+MEG+MEGPLANAR'};
modality.values = {
    {'All'}
    {'EEG'}
    {'MEG'}
    {'MEGPLANAR'}
    {'EEG', 'MEG'}
    {'MEG', 'MEGPLANAR'}
    {'EEG', 'MEGPLANAR'}
    {'EEG', 'MEG', 'MEGPLANAR'}
    }';
modality.val = {{'All'}};

invert = cfg_exbranch;
invert.tag = 'invert';
invert.name = 'Source inversion';
invert.val = {D, val, whatconditions, isstandard, modality};
invert.help = {'Run imaging source reconstruction'};
invert.prog = @run_inversion;
invert.vout = @vout_inversion;
invert.modality = {'EEG'};

function  out = run_inversion(job)

D = spm_eeg_load(job.D{1});

inverse = [];
if isfield(job.whatconditions, 'condlabel')
    inverse.trials = job.whatconditions.condlabel;
end

if isfield(job.isstandard, 'custom')
    inverse.type = job.isstandard.custom.invtype;
    inverse.woi  = fix([max(min(job.isstandard.custom.woi), 1000*D.time(1)) min(max(job.isstandard.custom.woi), 1000*D.time(end))]);
    inverse.Han  = job.isstandard.custom.hanning;
    inverse.lpf  =  fix(min(job.isstandard.custom.foi));
    inverse.hpf  =  fix(max(job.isstandard.custom.foi));

    P = char(job.isstandard.custom.priors.priorsmask);
    if ~isempty(P)
        [p,f,e] = fileparts(P);
        switch lower(e)
            case '.gii'
                g = gifti(P);
                inverse.pQ = cell(1,size(g.cdata,2));
                for i=1:size(g.cdata,2)
                    inverse.pQ{i} = double(g.cdata(:,i));
                end
            case '.mat'
                load(P);
                inverse.pQ = pQ;
            case {'.img', '.nii'}
                S.D = D;
                S.fmri = P;
                S.space = job.isstandard.custom.priors.space;
                D = spm_eeg_inv_fmripriors(S);
                inverse.fmri = D.inv{D.val}.inverse.fmri;
                load(inverse.fmri.priors);
                inverse.pQ = pQ;
            otherwise
                error('Unknown file type.');
        end
    end

    if ~isempty(job.isstandard.custom.restrict.locs)
        inverse.xyz = job.isstandard.custom.restrict.locs;
        inverse.rad = job.isstandard.custom.restrict.radius;
    end

    P = char(job.isstandard.custom.restrict.mask);
    if ~isempty(P)
        inverse.mask = P;
    end
end

[mod, list] = modality(D, 1, 1);
if strcmp(job.modality{1}, 'All')
    inverse.modality  = list;
else
    inverse.modality  = intersect(list, job.modality);
end

if numel(inverse.modality) == 1
    inverse.modality = inverse.modality{1};
end

D = {};

for i = 1:numel(job.D)
    D{i} = spm_eeg_load(job.D{i});

    D{i}.val = job.val;

    D{i}.con = 1;

    if ~isfield(D{i}, 'inv')
        error('Forward model is missing for subject %d.', i);
    elseif  numel(D{i}.inv)<D{i}.val || ~isfield(D{i}.inv{D{i}.val}, 'forward')
        if D{i}.val>1 && isfield(D{i}.inv{D{i}.val-1}, 'forward')
            D{i}.inv{D{i}.val} = D{i}.inv{D{i}.val-1};
            warning('Duplicating the last forward model for subject %d.', i);
        else
            error('Forward model is missing for subject %d.', i);
        end
    end

    D{i}.inv{D{i}.val}.inverse = inverse;
end

D = spm_eeg_invert(D);

if ~iscell(D)
    D = {D};
end

for i = 1:numel(D)
    save(D{i});
end

out.D = job.D;

function dep = vout_inversion(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'M/EEG dataset(s) after imaging source reconstruction';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
function TestMappingDNN()

modeldir = my_dir('nnet');
modelfiles = findFiles(['nnet/' modeldir{1}], 'mat');
modelfiles = sort_nnet_by_itr(modelfiles);      % get the last iteration model

dnn = load(modelfiles{1});
para = dnn.para;
layer = dnn.layer;

[~, Data_cv, para] = LoadData_AU4_Mapping(para);

para.out_layer_idx = length(layer) + [-2 -3];   % you can specify which layers' activation will be outputed
output = FeatureTree2(Data_cv, para, layer);

for i=1:length(output)
    enhanced{i} = output{i}{1};
end

for i=1:length(enhanced)
    subplot(3,1,1); imagesc(Data_cv(1).data{i});    title('Noisy MFCC');
    subplot(3,1,2); imagesc(enhanced{i});          title('Enhanced MFCC');
    subplot(3,1,3); imagesc(Data_cv(2).data{i});       title('Clean MFCC');
    pause
end

end

author: singaxiong

repo: SignalGraph


text:
function [audio_obj, features] = get_files_and_stft_features(input_dir, l, fn_filter)
    current_dir = cd();
    cd(input_dir);
    files = get_files_by_filter(fn_filter);
    % label according to l-th character in file name
    audio_obj = miraudio(files, 'Label', l, 'Normal');
    a = mirgetdata(audio_obj);
    sr = get(audio_obj, 'Sampling');
    features = cell(length(a), 1);
    for i = 1:length(a)
        f = calc_stft(a{i}, sr{i});
%         f = lpc(a{i});
%         f = analyze_samples(a{i}, sr{i});
        features{i} = f;
%         features{i} = mfcc(a{i}, sr{i}, 512);
    end
    cd(current_dir);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function  [slcstack,  interfstack] = ImgRead_isce(InSAR_path,nline,reference_date,slcslist)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This file is part of TomoSAR.
%
%   TomoSAR is distributed in the hope that it will be useful,
%   but without warranty of any kind; without even the implied
%   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%   See the Apache License for more details.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author : Dinh Ho Tong Minh (INRAE) and Yen Nhi Ngo, Mar. 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N_int = size(slcslist,1);

for ii=1:N_int
    interfstack.datastack(:,:,ii)=single(freadbkj([InSAR_path,'/',num2str(slcslist(ii)),'/isce_minrefdem.int'],nline,'cpxfloat32'));
    interfstack.filename(ii,2)=slcslist(ii);

    tempstack(:,:,ii)=single(freadbkj([InSAR_path,'/',num2str(slcslist(ii)),'/secondary.slc'],nline,'cpxfloat32'));

    fprintf('Reading slc/interferogram %d / %d \n',ii,N_int);
end
interfstack.filename(:,1) = reference_date;

% insert reference_date in to slc list
slcslist = [reference_date; slcslist];
slcslist_datenum = sort(datenum(num2str(slcslist),'yyyymmdd'));

slcslist = datestr(slcslist_datenum,'yyyymmdd');
N_slc = size(slcslist,1);

[~,reference_ind]=ismember(datenum(num2str(reference_date),'yyyymmdd'),slcslist_datenum);
if reference_ind > 1
    slcstack.datastack(:,:,[1:reference_ind-1,reference_ind+1:N_slc]) = tempstack;
    slcstack.datastack(:,:,reference_ind) = single(freadbkj([InSAR_path,'/',num2str(slcslist(1,:)),'/reference.slc'],nline,'cpxfloat32'));
else
    slcstack.datastack(:,:,1) = single(freadbkj([InSAR_path,'/',num2str(slcslist(2,:)),'/reference.slc'],nline,'cpxfloat32'));
    slcstack.datastack(:,:,[2:N_slc]) = tempstack;
end

for ii = 1:N_slc
    slcstack.filename(ii,1) = str2double(slcslist(ii,:));
end

author: DinhHoTongMinh

repo: TomoSAR


text:
syms x
I=int(1/(1+sqrt(1-x^2)))
pretty(I)

author: Eurus-Holmes


text:
function [assocScansV, relStructNumV] = getStructureAssociatedScan(structsV, planC)
%"getStructureAssociatedScan"
%   Returns a vector with the corresponding associated scan for each
%   structure number passed in structsV.
%
%   Structures without an associatedScan field, structures with an empty
%   associatedScan field are assumed to be associated with scan 1.
%
%   Also returns relStructNum, a list of "relative" structure indices for
%   the passed structsV.  This is the number of the structure in the list
%   of all structures with the same associatedScan.
%
%   If planC is not specified, the global planC is used.
%
%JRA 2/15/05
%
%Usage:
%   [assocScansV, relStructNum] = getStructureAssociatedScan(structsV, planC)

% if ~exist('planC')
%     global planC
% end
indexS = planC{end};

%Return if no structures were passed in.
if isempty(structsV)
    assocScansV = [];
    relStructNumV = [];
    return;
end

%If no associated scan field, assume scan 1.
if ~isfield(planC{indexS.structures}, 'associatedScan')
    allAssocScansV = ones(1, length(planC{indexS.structures}));
else
    %If field is blank, assume scan 1, else take field contents.
    for i=1:length(planC{indexS.structures})
        if isempty(planC{indexS.structures}(i).associatedScan)
            allAssocScansV(i) = 1;
        else
            allAssocScansV(i) = planC{indexS.structures}(i).associatedScan;
        end
    end
end

%Return only requested structures' assocScans.
assocScansV = allAssocScansV(structsV);

%Calculate relativeStructNum for all structures.
allScans = unique(allAssocScansV);
for i=1:length(allScans)
    scanNum = allScans(i);
    matchingIndV = find(allAssocScansV == scanNum);
    allRelStructNumV(matchingIndV) = 1:length(matchingIndV);
end

%Only return relStructNum for the requested structs.
relStructNumV = allRelStructNumV(structsV);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function g = compute_gratio(D)

% Ikeda M, Oka Y. Brain Behav, 2012. "The relationship between nerve conduction velocity and fiber morphology during peripheral nerve regeneration."
g = 0.220 .* log10(2*D) + 0.508;
% g = 0.76; % if you want a constant g-ratio

% figure
% plot(2.*R, g)

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function bookplot(fig, plotfile)
%bookplot   saves IFISS figure as postscript file
%   bookplot(fig, plotfile);
%   input
%          fig        figure number of figure to be saved
%          plotfile   character string naming the file
%
%   uses the MATLAB script exportfig
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
fprintf('\nCreating plotfiles ...')
gohome, cd plotfiles
figure(fig)
exportfig(gcf,[plotfile,'ef.eps'],'width',4.5,'fontmode','fixed',...
          'fontsize',8,'color','bw')
fprintf('\n%s',[plotfile,'ef.eps'])
exportfig(gcf,[plotfile,'ef.epsc'],'width',4.5,'fontmode','fixed',...
          'fontsize',8,'color','cmyk')
fprintf('\n%s\n',[plotfile,'ef.epsc'])
return

author: johannesgerer

repo: jburkardt-m


text:
function im = loadImage(filename)
%LOADIMAGE      Load an image file.
%
% DESCRIPTION:
%       loadImage loads an external image file using imread and returns a
%       two-dimensional image matrix scaled between 0 and 1. If a colour
%       image is loaded, the colour channels are summed before scaling. If
%       no input is given for filename, a selection dialog is invoked.
%
% USAGE:
%       im = loadImage()
%       im = loadImage(filename)
%
% OPTIONAL INPUTS:
%       filename    - filename (and pathname if not in the same directory)
%                     of the image to load
%
% OUTPUTS:
%       im          - scaled image matrix
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 30th June 2009
%       last update - 25th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also imread, resize

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

% load the file if not selected
if nargin == 0
    [filename, pathname] = getFilename();
    filename = strcat(pathname, filename);
end

% extract the filetype
split_filename = regexp(filename, '\.', 'split');
filetype = lower(split_filename{2});

% check the filetype
if ~strcmp(filetype, 'jpg') && ~strcmp(filetype, 'bmp') && ~strcmp(filetype, 'gif')  && ~strcmp(filetype, 'png')
    disp('WARNING: image filetype may not be supported - see loadImage.m');
end

% load the image
im = imread(filename);

% if there are multiple colour channels, sum them up
if numDim(im) > 2
    im = squeeze(double(im(:, :, 1)) + double(im(:, :, 2)) + double(im(:, :, 3)));
else
    im = double(im);
end

% scale pixel values from 0 -> 1
im = max(im(:)) - im;
im = im .* (1/max(im(:)));

repo: Aero-matlab


text:
function l = readascii(filename)

% l = readascii(filename)
% read each line of an ASCII file into list l

l = list;		% create empty list
fid = fopen(filename, 'r');

if fid ~= -1
	while 1
		line = fgetl(fid);
		if ~isstr(line), break, end
			l = TSTL_append(l, line);
		end
	fclose(fid);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function res = badchannels(this, varargin)
% Method for getting/setting bad channels
% FORMAT res = badchannels(this)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: badchannels.m 4329 2011-05-19 20:54:36Z christophe $

if length(varargin) == 2 && ~isempty(varargin{1})
    % make sure that the two inputs for set are the same length
    if ~(length(varargin{2}) == 1 | (length(varargin{1}) == length(varargin{2})))
        error('Use either same vector length or scalar for value');
    end
end

if numel(varargin) >= 1  && ~isempty(varargin{1})
    if ~(all(varargin{1} >= 1) && all(varargin{1} <= nchannels(this)))
        error('Channel number out of range.');
    end
end

if numel(varargin) >= 2
    ubad = unique(varargin{2});
    if isempty(ubad) | ~all(ismember(ubad, [0 1]))
        error('Illegal bad flags (should be 0 or 1)');
    end
end

res = getset(this, 'channels', 'bad', varargin{:});

if isempty(varargin)
    if iscell(res)
        res = [res{:}];
    end

    res = find(res);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function save_hdf5(ndx,outfilename)
% Saves the Ndx object in an hdf5 file.
% Inputs:
%   ndx: The Ndx object to save.
%   outfilename: The name of the output hdf5 file.

assert(nargin==2)
assert(isa(ndx,'Ndx'))
assert(isa(outfilename,'char'))
assert(ndx.validate())

nummods = length(ndx.modelset);
numsegs = length(ndx.segset);

assert(size(ndx.trialmask,1)==nummods)
assert(size(ndx.trialmask,2)==numsegs)

trialmask = uint8(ndx.trialmask);

hdf5write(outfilename,'/trial_mask',trialmask,'V71Dimensions',true);
hdf5write(outfilename,'/ID/row_ids',ndx.modelset,'V71Dimensions',true,'WriteMode','append');
hdf5write(outfilename,'/ID/column_ids',ndx.segset,'V71Dimensions',true,'WriteMode','append');

author: nesl


text:
% Digital circuit sizing for an inverter chain (GP)
% Boyd, Kim, Patil, and Horowitz, "Digital circuit optimization
% via geometric programming"
% Written for CVX by Almir Mutapcic 02/08/06
%
% We consider a chain of N inverters driving a load capacitance CL.
% The problem is to find optimal scale factors for the inverter
% that minimize the sum of them (area), while obeying constraints
% on the maximum delay through the circuit, and minimum and maximum
% limits on scale factors. There are no limits on the total power.
% (For more details about the inverter chain see sec. 2.1.11 in the paper.)
%
%   minimize   sum(x)
%       s.t.   T_j <= Dmax          for j an output gate
%              T_j + d_i <= T_i     for j in FI(i)
%              x_min <= x <= x_max
%
% where variables are x and T.
% Here we use data structures and digital circuit models from the
% referenced paper.

%********************************************************************
% problem data
%********************************************************************
N  = 8;      % number of inverters
CL = 20;     % capacitance load
Dmax = 20;   % maximum delay through the circuit
x_min = 1;   % minimum scale factor
x_max = 20;  % maximum scale factor

% circuit labeling convention:
% label primary input (input to the first inverter in the chain) with N+1
% label primary output (output of the last inverter in the chain) with N+2
% label inverters in the chain with 1,2,...,N based on their location

% primary input and primary output labels (start with N+1)
primary_inputs  = [N+1];
primary_outputs = [N+2];
M = N + length( primary_inputs ) + length( primary_outputs );

% fan-in cell array for a straight chain of inverters
FI{1} = [N+1];   % fan-in of the first inverter is the primary input
for k = 2:N
  FI{k} = [k-1]; % fan-in of other inverters is the inverter feeding into them
end
FI{N+2} = [N];   % fan-in of the primary output is the last inverter in the chain

% fan-out cell array
% (will be computed from the fan-in cell array, no need to modify)
FO = cell(M,1);
for gate = [1:N primary_outputs]
  preds = FI{gate};
  for k = 1:length(preds)
    FO{preds(k)}(end+1) = gate;
  end
end

% input and internal capacitance of gates and the driving resistance
Cin_norm  = ones(N,1);
Cint_norm = ones(N,1);
Rdrv_norm = ones(N,1);

% place extra capacitance before the input of the 5th inverter
Cin_norm(5) = 80;

% primary output has Cin capacitance (but has no Cload)
Cin_po = sparse(M,1);
Cin_po(primary_outputs) = CL;

% primary input has Cload capacitance (but has no Cin)
Cload_pi = sparse(M,1);
Cload_pi(primary_inputs) = 1;

%********************************************************************
% optimization
%********************************************************************
cvx_begin gp
  % optimization variables
  variable x(N)                 % sizes
  variable T(N)                 % arrival times

  % minimize the sum of scale factors subject to above constraints
  minimize( sum(x) )
  subject to

    % input capacitance is an affine function of sizes
    Cin  = Cin_norm.*x;
    Cint = Cint_norm.*x;

    % driving resistance is inversily proportional to sizes
    R = Rdrv_norm./x;

    % gate delay is the product of its driving resistance and load cap.
    Cload = cvx( zeros(N,1) );
    for gate = 1:N
      if ~ismember( FO{gate}, primary_outputs )
        Cload(gate) = sum( Cin(FO{gate}) );
      else
        Cload(gate) = Cin_po( FO{gate} );
      end
    end

    % delay
    D = 0.69*ones(N,1).*R.*( Cint + Cload );

    % create timing constraints
    for gate = 1:N
      if ~ismember( FI{gate}, primary_inputs )
        for j = FI{gate}
          % enforce T_j + D_j <= T_i over all gates j that drive i
          D(gate) + T(j) <= T(gate);
        end
      else
        % enforce D_i <= T_i for gates i connected to primary inputs
        D(gate) <= T(gate);
      end
    end

    % circuit delay is the max of arrival times for output gates
    output_gates = [FI{primary_outputs}];
    circuit_delay = max( T(output_gates) );

    % collect all the constraints
    circuit_delay <= Dmax;
    x_min <= x <= x_max;
cvx_end

% message about extra capacitance and result display
disp(' ')
disp(['Note: there is an extra capacitance between the 4th and 5th inverter'...
     ' in the chain.'])
fprintf(1,'\nOptimal scale factors are: \n'), x

% plot scale factors and maximum delay for inverter i
close all;
subplot(2,1,1); plot([1:N],T,'g--',[1:N],T,'bo');
ylabel('maximum delay T')
subplot(2,1,2); stem([1:N],x);
ylabel('scale factor x')
xlabel('inverter stage')

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function c=ref_dstiii(f)
%REF_DSTIII  Reference Discrete Sine Transform type III
%   Usage:  c=ref_dstiii(f);
%
%   This is the inverse of DSTII

L=size(f,1);
W=size(f,2);

% Create weights.
w=ones(L,1);
w(L)=1/sqrt(2);
w=w*sqrt(2/L);

% Create transform matrix.
F=zeros(L);

for m=0:L-1
  for n=0:L-1
    F(m+1,n+1)=w(n+1)*sin(pi*(n+1)*(m+.5)/L);
  end;
end;

% Compute coefficients.
c=F*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% Downsampling procedure.
%
% Arguments:
%   'I': image
%   downsampling filter 'filter', should be a 2D separable filter.
%   'border_mode' should be 'circular', 'symmetric', or 'replicate'. See 'imfilter'.
%   subwindow indices 'subwindow', given as [r1 r2 c1 c2] (optional)
%
% tom.mertens@gmail.com, August 2007
% sam.hasinoff@gmail.com, March 2011  [handle subwindows, reweighted boundaries]
%

function [R,subwindow_child] = downsample(I, filter, subwindow)

r = size(I,1);
c = size(I,2);
if ~exist('subwindow','var')
    subwindow = [1 r 1 c];
end
subwindow_child = child_window(subwindow);

border_mode = 'reweighted';
%border_mode = 'symmetric';

switch border_mode
    case 'reweighted'
        % low pass, convolve with 2D separable filter
        R = imfilter(I,filter);

        % reweight, brute force weights from 1's in valid image positions
        Z = imfilter(ones(size(I)),filter);
        R = R./Z;

    otherwise
        % low pass, convolve with 2D separable filter
        R = imfilter(I,filter,border_mode);
end

% decimate
reven = mod(subwindow(1),2)==0;
ceven = mod(subwindow(3),2)==0;
R = R(1+reven:2:r, 1+ceven:2:c, :);

end

author: drakeguan


text:
function [ newp ] = tsgRecycleGrid( lGrid, depth, order, type, anisotropy )
%
% [ newp ] = tsgRecycleGrid( lGrid, depth, order, type, anisotropy )
%
% overwrites the grid in lGrid with a new grid having the same parameters,
% but with updated depth, order, type and anisotropy. Recycling also reuses
% all of the old points values.
%
% INPUT:
%
% lGrid: a grid list created by tsgMakeGrid(...)
%
% depth: (all grids) integer bigger than or equal to 0 (0 is accepted only by type = basis)
%        the deinsity of the grid, i.e. levels
%
% order: (local polynomial or wavelet grids only) integer 1, 2 or 3 (2 is for polynomial grids only)
%        the order of the grid, i.e. linear, quadratic, cubic
%
% type: (global grids only) string from <level,basis>
%       level used the classical Smolyak tensor selection (depth == level)
%       basis gives the minimum subset of the tensors requires to capture a
%             polynomial of certain degree,
%             clenshaw-curtis and chebychev consider interolation
%             gauss-legendre focuses on interation
%
% anisotropy: (global grids only) vector of size dim with min( anisotropy ) == 1
%             gives linear weighting on the selection of the tensors
%
% OUTPUT:
%
% newp: (optional) the new points associated with the new grid in an array
%       of dimension [ num_new_poits, dim ]
%
% NOTE: lGrid is overwritten
%

if ( nargin < 8 )
    anisotropy = [];
end

[ sFiles, sTasGrid ] = tsgGetPaths();
[ gFile, qFile, vFile, wFile, pFile, rFile ] = tsgMakeFilenames( lGrid.name );

sCommand = [sTasGrid,' -recycle'];

sCommand = [ sCommand, ' -gridfile ', gFile];

sCommand = [ sCommand, ' -depth ',num2str(depth)];

sCommand = [ sCommand, ' -order ',num2str(order)];

sCommand = [ sCommand, ' -type ',type];

if ( size( anisotropy ) ~= [ 0 0 ] )
    tsgWriteMatrix( wFile, anisotropy );
    sCommand = [ sCommand, ' -af ',wFile];
end

if ( nargout > 1 )
    sCommand = [ sCommand, ' -of ',qFile];
end

[status, cmdout] = system(sCommand);

if ( size( findstr( 'ERROR', cmdout ) ) ~= [ 0, 0 ] )
    disp(cmdout);
    fprintf(2,['Erors were encountered\n']);
    return;
else
    if ( nargout > 1 )
        newp = tsgReadMatrix( qFile );
    end
end

end


text:
function [videoMeta] = readInVideoCsv(filename)
% Read the csv and convert it into a table
csvData = readtable(filename, 'Delimiter', ',');
% Initialize the tracks struct
videoMeta = {};

% Declare and initialize the extracted_meta_dictionary
videoMeta.id = int32(csvData.id);
videoMeta.frameRate = int32(csvData.frameRate);
videoMeta.locationId = int32(csvData.locationId);
videoMeta.speedLimit = double(csvData.speedLimit);
videoMeta.month = char(csvData.month);
videoMeta.weekDay = char(csvData.weekDay);
videoMeta.startTime = char(csvData.startTime);
videoMeta.duration = char(csvData.duration);
videoMeta.totalDrivenDistance = double(csvData.totalDrivenDistance);
videoMeta.totalDrivenTime = double(csvData.totalDrivenTime);
videoMeta.numVehicles = int32(csvData.numVehicles);
videoMeta.numCars = int32(csvData.numCars);
videoMeta.numTrucks = int32(csvData.numTrucks);
try
    % Extract upper Lanes
    upperCell = regexp(csvData.upperLaneMarkings, ';', 'split');
    upperCell = upperCell{1};
    cellSize = size(upperCell, 2);
    upperLanes = double.empty(cellSize, 0);
    for i = 1:size(upperCell, 2)
        cellValue = upperCell(1, i);
        upperLanes(i) = str2double(cellValue);
    end
    videoMeta.upperLanes = upperLanes;

    % Extract lower Lanes
    lowerCell = regexp(csvData.lowerLaneMarkings, ';', 'split');
    lowerCell = lowerCell{1};
    cellSize = size(lowerCell, 2);
    lowerLanes = double.empty(cellSize, 0);
    for i = 1:size(lowerCell, 2)
        cellValue = lowerCell(1, i);
        lowerLanes(i) = str2double(cellValue);
    end
    videoMeta.lowerLanes = lowerLanes;
catch
    printf("Something failed when reading the lane markings.");
end
end

author: RobertKrajewski

repo: highD-dataset


text:
function poisson_simulation_test02 ( )

%*****************************************************************************80
%
%% POISSON_SIMULATION_TEST02 simulates waiting for a given length of time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 September 2012
%
%  Author:
%
%    John Burkardt
%
  lambda = 0.5;
  t = 1000.0;
  test_num = 20000;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POISSON_SIMULATION_TEST02:\n' );
  fprintf ( 1, '  POISSON_FIXED_EVENTS simulates a Poisson process\n' );
  fprintf ( 1, '  counting the number of events that occur during\n' );
  fprintf ( 1, '  a given time.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Simulate a Poisson process, for which, on average,\n' );
  fprintf ( 1, '  LAMBDA events occur per unit time.\n' );
  fprintf ( 1, '  Run for a total of %g time units.\n', t );
  fprintf ( 1, '  LAMBDA = %g\n', lambda );

  n = zeros ( test_num, 1 );
  for test = 1 : test_num
    n(test) = poisson_fixed_time ( lambda, t );
  end

  n_mean = mean ( n );
  n_var = var ( n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Mean number of events = %g\n', n_mean );
  fprintf ( 1, '  Variance = %g,  STD = %g\n', n_var, sqrt ( n_var ) );

  bins = 30;
  hist ( n, bins );
  xlabel ( '<--Number of events-->' )
  ylabel ( '<--Frequency-->' )
  grid on
  title ( ' Number of Poisson events observed over 1000 time units' );

  filename = 'poisson_events.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file saved as "%s"\n', filename );

  return
end


text:
function varargout = ellipseToPolygon(ellipse, N)
%ELLIPSETOPOLYGON Convert an ellipse into a series of points
%
%   P = ellipseToPolygon(ELL, N);
%   converts ELL given as [x0 y0 a b] or [x0 y0 a b theta] into a polygon
%   with N edges. The result P is a N-by-2 array containing the coordinates
%   of the N vertices of the polygon.
%
%   P = ellipseToPolygon(ELL);
%   Use a default number of edges equal to 72. This results in one point
%   for each 5 degrees.
%
%   [X, Y] = ellipseToPolygon(...);
%   Return the coordinates of vertices in two separate arrays.
%
%   Example
%     poly = ellipseToPolygon([50 50 40 30 20], 60);
%     figure; hold on;
%     axis equal; axis([0 100 10 90]);
%     drawPolygon(poly, 'b');
%     drawPoint(poly, 'bo');
%
%   See also:
%   ellipses2d, drawEllipse, circleToPolygon, rectToPolygon
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/04/2005.
%

% HISTORY
% 2011-03-30 use angles in degrees, add default value for N
% 2011-12-09 rename to ellipseToPolygon
% 2017-08-31 returns N vertices instead of N+1

% default value for N
if nargin < 2
    N = 72;
end

% angle of ellipse
theta = 0;
if size(ellipse, 2) > 4
    theta = ellipse(:,5);
end

% get ellipse parameters
xc = ellipse(:,1);
yc = ellipse(:,2);
a  = ellipse(:,3);
b  = ellipse(:,4);

% create time basis
t = linspace(0, 2*pi, N+1)';
t(end) = [];

% pre-compute trig functions (angles is in degrees)
cot = cosd(theta);
sit = sind(theta);

% position of points
x = xc + a * cos(t) * cot - b * sin(t) * sit;
y = yc + a * cos(t) * sit + b * sin(t) * cot;

% format output depending on number of a param.
if nargout == 1
    varargout = {[x y]};
elseif nargout == 2
    varargout = {x, y};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function x = bin2hex(BinStr)

% BIN2HEX (BinStr,...output)
%
% Converts binary strings of any length to hexadecimal pairs. Adds leading
% zeros if there is not an even number of hex bits.
%
% Now also supports cell array inputs.
%
%
% Author: Richard Medlock 2001

% Check to see what the input is:

if iscell(BinStr)

    NCells = length(BinStr)

    for C = 1:NCells

        Str = BinStr{C};

        x{C} = doconvert(Str);

    end

elseif ~ischar(BinStr)

    error('Input must be a string or character or cell array.')

else

    x = doconvert(BinStr);

end

%-----------------------------------------------
function out = doconvert(BinStr)

Bits = length(BinStr);
Words = Bits/4;
WholeWords = floor(Words);
PartWords = Words-WholeWords;

while PartWords > 0

    BinStr = ['0' BinStr];
    Bits = length(BinStr);
    Words = Bits/4;
    WholeWords = floor(Words);
    PartWords = Words-WholeWords;

end

% For each WORD (4-bits), convert it to HEX based on the following rules:

Words = length(BinStr)/4;

HEXi = {'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'};

HexStr = [];

for W = 1:Words

    %For each word, convert it first to decimal:

    P = (W*4)-3;

    Word = BinStr(P:P+3);
    DEC = bin2dec(Word);
    HEX = HEXi{DEC+1};

    HexStr = [HexStr HEX];

end

LHex = length(HexStr); %If the Hex String has an odd number of bits, add 0 to the front.

if isodd(LHex)

    HexStr = ['0' HexStr];

end

out = HexStr;

author: Sable

repo: mcbench-benchmarks


text:
% example for convolvepoints.m
%
% First you need to generate psf (just a small image 11x11 (should be odd numbers...) pixels in enough):
% eg:
psf = psfgen('lambda', 600, 'na', 1.2, 'pixelsize', 110, 'sizevec', [11 11], 'method', 'airy');

% your x and y are column vectors of coordiantes. If units of x and y are
% [nm], then set pixelsize = 1. If units of x and y are in pixels then set
% pixelsize = [size of the pixel in nm].
% Here I generate radnom coordinates of N=1000 points with values in nm (in the range 0:10000 nm)
N = 500;
pixelsize = 100; % in nanometers
x = 1e4*rand(N,1); % N random coordiantes (in nm) in the range 0:10000
y = 1e4*rand(N,1);
np = 100 + 20*randn(N,1); % N random intesnities with mean 100 pixels and variance 50

% generate the image:
imout1 = convolvepoints(x,y,psf,np,pixelsize);
figure;
imagesc(imout1)
title('imout1')
axis off
set (gca, 'DataAspectRatio',[1 1 1]);
colormap('gray')

% if all intensities are the same than you can pass np as a single number
% eg:
np = 1000;
% this will crate an image where all psf have the same intesnity 1000;
imout2 = convolvepoints(x,y,psf,np,pixelsize);
figure;
imagesc(imout2)
title('imout2')
axis off
set (gca, 'DataAspectRatio',[1 1 1]);
colormap('gray')

% if you don't pass np and pixelsize argumetn then defaulet values (np =1
% and pixelsize = 1 are used. This is for the case where all psf have same
% intensity and x and y are in pixels:
% eg: to create image with psf randomly scattered in an image 100x100
% pixels:
x = 100*rand(N,1); % in pixels
y = 100*rand(N,1);
imout3 = convolvepoints(x,y, psf);
% this is equivalent to:
% convolvepoints(x,y, psf, 1,1);
figure;
imagesc(imout3)
title('imout3')
axis off
set (gca, 'DataAspectRatio',[1 1 1]);
colormap('gray')

% This is for showing the location of the points...
[imout4, xout, yout] = convolvepoints(x,y, psf);
figure;
imagesc(imout4)
hold on
scatter(xout, yout,'r')
title('imout4')
axis off
set (gca, 'DataAspectRatio',[1 1 1]);
colormap('gray')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef prtRvMultinomial < prtRv & prtRvMemebershipModel
    % prtRvMultinomial  Multinomial random variable
    %
    %   RV = prtRvMultinomial creates a prtRvMultinomial object
    %   with an unknown number of categories with unspecified probabilities
    %   These properties can be set manually or by using the MLE method.
    %
    %   prtRvMultinomial operates on count matrices. Therfore, the DRAW()
    %   method outputs a matrix that is N x nCategories and has a single 1
    %   in each row. To draw integer categories you can use the
    %   DRAWINTEGER() method.  Similarly, the MLE function takes a count
    %   matrix as an input. Type help prtRvMultinomial.mle for more
    %   information.
    %
    %   RV = prtRvMultinomial(PROPERTY1, VALUE1,...) creates a
    %   prtRvMultinomial object RV with properties as specified by
    %   PROPERTY/VALUE pairs.
    %
    %   A prtRvMultinomial object inherits all properties from the
    %   prtRv class. In addition, it has the following properties:
    %
    %   nCategories   - number of integers modeled by the RV
    %   probabilities - A 1 x nCategories vector of doubles less than 1
    %                   that sum to 1, representing the probability of
    %                   each of the integers
    %
    %  A prtRvMultinomial object inherits all methods from the prtRv
    %  class. The MLE  method can be used to set the parameters from data.
    %  In addition, it has the the following methods:
    %
    %   x = R.drawIntegers(N) - Draws N integers with the corresponding
    %                           probabilities
    %
    %  Example:
    %
    %  data = rand(100,5);                  % Uniformly random data
    %  X = bsxfun(@eq,data,max(data,[],2)); % Generate data that has a
    %                                       % single 1 in each row
    %
    %  RV = prtRvMultinomial;               % Generate a prtRvMultinomial
    %  RV = mle(RV,X);                      % Estimate the parameters
    %
    %  RV.plotPdf()                         % Plot the pdf (pmf)
    %
    %   See also: prtRv, prtRvMvn, prtRvGmm, prtRvVq, prtRvKde

    properties (SetAccess = private)
        name = 'Multinomial Random Variable';
        nameAbbreviation = 'RVMulti';
    end

    properties (SetAccess = protected)
        isSupervised = false;
        isCrossValidateValid = true;
    end

    properties (Dependent)
        probabilities
    end

    properties (Dependent = true)
        nCategories  % The number of categories
    end

    properties (Hidden = true, SetAccess='private', GetAccess='private')
        probabilitiesDepHelper
    end

    properties (Hidden = true, Dependent = true)
        nDimensions
    end

    properties (Hidden = true)
        approximatelyEqualThreshold = 1e-4;
        nDrawsPerObservationDraw = [];
    end

    methods
        % The Constructor
        function R = prtRvMultinomial(varargin)
            R = constructorInputParse(R,varargin{:});
        end

        function R = mle(R,X)
            % RV = RV.mle(X) computes the maximum likelihood estimate based
            % the data X. X should be nObservations x nDimensions. X must
            % be a count matrix, consisting of only zeros and ones. A one
            % in the ith column indicates that the sample in the jth row is
            % of class i.
            %

            %             if ~isequal(unique(X(:)),[0;1])
            %                 error('prtRvMultinomial:invalidData','Input matrix must contain only ones and zeros');
            %             end

            % Do we really want to or need to enforce this?
            % Can't we allow learning of the probability from multinomial
            % draws that had N greater than 1?
            % I don't see why not. - KDM 2013-03-01
            %if ~prtUtilApproxEqual(sum(X,2),1,1e-6);
            %    error('prtRvMultinomial:invalidData','Rows of input data must contain no more than one "1"');
            %end

            X = R.dataInputParse(X); % Basic error checking etc

            N_bar = sum(X,1);
            R.probabilities = N_bar./sum(N_bar(:));
        end

        function vals = pdf(R,X)
            % PDF Output the pdf of the random variable evaluated at the points specified
            %
            % pdf = RV.pdf(X) returns  the pdf of the prtRv
            % object evaluated at X. X must be an N x nDims matrix, where
            % N is the number of locations to evaluate the pdf, and nDims
            % is the same as the number of dimensions, nDimensions, of the
            % prtRv object RV.

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'PDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            X = R.dataInputParse(X); % Basic error checking etc
            assert(isnumeric(X) && ndims(X)==2,'X must be a 2D numeric array.');
            assert(size(X,2) == R.nCategories,'Incorrect number of categories for this RV object. This RV object is defined to have %d categories, but the input data has only %d columns. Remember that prtRvMultinomial operates on count matrices.', R.nCategories, size(X,2))

            vals = sum(bsxfun(@times,X,R.probabilities),2);
        end

        function vals = logPdf(R,X)
            % LOGPDF Output the log pdf of the random variable evaluated at the points specified
            %
            % logpdf = RV.logpdf(X) returns the logarithm of value of the
            % pdf of the prtRv object evaluated at X. X must be an N x
            % nDims matrix, where N is the number of locations to evaluate
            % the pdf, and nDims is the same as the number of dimensions,
            % nDimensions, of the prtRv object RV.

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'LOGPDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            vals = log(pdf(R,X));
        end

        function vals = draw(R,N)
            % DRAW  Draw random samples from the distribution described by the prtRv object
            %
            % VAL = RV.draw(N) generates N random samples drawn from the
            % distribution described by the prtRv object RV. VAL will be a
            % N x nDimensions vector, where nDimensions is the number of
            % dimensions of RV.

            assert(numel(N)==1 && all(N==floor(N)) && all(N > 0),'N must be a positive integer scalar.')

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'DRAW cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            if ~isempty(R.nDrawsPerObservationDraw)
                nObsToDraw = R.nDrawsPerObservationDraw;
                vals = zeros(N,R.nCategories);
                for iBag = 1:N
                    cVals = zeros(nObsToDraw,R.nCategories);
                    cVals(sub2ind([nObsToDraw, R.nCategories], (1:nObsToDraw)',drawIntegers(R,nObsToDraw))) = 1;
                    vals(iBag,:) = sum(cVals,1);
                end
            else
                vals = zeros(N,R.nCategories);
                vals(sub2ind([N, R.nCategories], (1:N)',drawIntegers(R,N))) = 1;
            end
        end

        function vals = drawIntegers(R,N)
            % DRAW  Draw random integer samples from the distribution described by the prtRv object
            %
            % VAL = RV.draw(N) generates N random integer samples drawn from the
            % distribution described by the prtRv object RV. VAL will be a
            % N x nDimensions vector, where nDimensions is the number of
            % dimensions of RV.

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'DRAWINTEGERS cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            assert(numel(N)==1 && N==floor(N) && N > 0,'N must be a positive integer scalar.')

            vals = prtRvUtilRandomSample(R.probabilities, N);
        end

        function varargout = plotPdf(R,varargin)
            %plotPdf Plot the pdf of the RV
            %
            % rv.plotPdf() plots the pdf of rv
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'plotPdf cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            h = bar(1:R.nCategories,R.probabilities,'k');
            ylim([0 1])
            xlim(R.plotLimits());

            varargout = {};
            if nargout
                varargout = {h};
            end
        end

        function plotCdf(R,varargin) %#ok<MANU>
            % plotCDF Not implemented for this prtRv
            error('prt:prtRvMultinomial','plotCdf is not implimented for this prtRv');
        end

    end

    % Set methods
    methods
        function R = set.probabilities(R,probs)
            assert(abs(sum(probs)-1) < R.approximatelyEqualThreshold,'Probability vector must must sum to 1.')

            R.probabilitiesDepHelper = probs(:)';
        end
        function R = set.nCategories(R,val) %#ok<MANU,INUSD>
            error('prt:prtRvMultinomial','nCategories is a dependent property that cannot be set by the user. To set the number of categories, set "probabilities" to be a vector of the desired length.');
        end
    end

    % Get methods
    methods
        function val = get.probabilities(R)
            val = R.probabilitiesDepHelper;
        end
        function val = get.nCategories(R)
            if ~isempty(R.probabilities)
                val = length(R.probabilities);
            else
                val = [];
            end
        end
        function val = get.nDimensions(R)
            val = R.nCategories;
        end
    end

    methods (Hidden = true)
        function [val, reasonStr] = isValid(R)
            if numel(R) > 1
                val = false(size(R));
                for iR = 1:numel(R)
                    [val(iR), reasonStr] = isValid(R(iR));
                end
                return
            end

            val = ~isempty(R.probabilities);

            if val
                reasonStr = '';
            else
                reasonStr = 'because probabilities has not been set';
            end
        end

        function val = plotLimits(R)
            [isValid, reasonStr] = R.isValid;
            if isValid
                val = [0.5 0.5+R.nCategories];
            else

                error('multinomial:plotLimits','Plotting limits can not yet be determined. This RV is not yet valid %s.',reasonStr)
            end
        end

        function val = isPlottable(R) %#ok
            val = true; % Always plottable
        end

        function R = weightedMle(R,X,weights)
            assert(numel(weights)==size(X,1),'The number of weights must mach the number of observations.');

            weights = weights(:);

            N_bar = sum(bsxfun(@times,X,weights),1);

            R.probabilities = N_bar./sum(N_bar(:));
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function title = p18_title ( title )

%*****************************************************************************80
%
%% P18_TITLE returns the title of problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = '10^14 * (x-1)^7, written term by term.';

  return
end


text:
function [x,fval,exitflag,info,Opt] = opti_lsqcurvefit(fun,x0,xdata,ydata,lb,ub,opts)
%OPTI_LSQCURVEFIT Solve a NLS using an OPTI NLS Solver (Matlab Overload)
%
%   [x,fval,exitflag,info] = opti_lsqcurvefit(fun,x0,xdata,ydata,lb,ub)
%   solves the nonlinear least squares problem sum[(f(x,xdata)-ydata)^2]
%   where fun is the nonlinear function to be fitted [fun(x,xdata)],
%   starting at x0, to the sample data ydata. Optional bounds lb and ub can
%   be placed on the decision variables x.
%
%   [x,fval,exitflag,info] = opti_lsqcurvefit(fun,...,ub,opts) allows the
%   user to specify optiset options. This includes specifying a solver via
%   the 'solver' field of optiset.
%
%   [x,...,info,Opt] = opti_lsqcurvefit(fun,...) returns the internally
%   built OPTI object.

%   Copyright (C) 2011 Jonathan Currie (IPL)

% Handle missing arguments
if nargin < 7, opts = optiset; end
if nargin < 6, ub = []; end
if nargin < 5, lb = []; end
if nargin < 4, error('You must supply at least 4 arguments to opti_lsqcurvefit'); end

%Build OPTI Object
Opt = opti('fun',fun,'data',xdata,ydata,'bounds',lb,ub,'x0',x0,'options',opts);

%Solve
[x,fval,exitflag,info] = solve(Opt);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [offsetsRows, offsetsCols, distances] = templateMatchingIntegralImage(row,...
    col,patchSize, searchWindowSize, image)
% This function should for each possible offset in the search window
% centred at the current row and col, save a value for the offsets and
% patch distances, e.g. for the offset (-1,-1)
% offsetsX(1) = -1;
% offsetsY(1) = -1;
% distances(1) = 0.125;

% The distance is simply the SSD over patches of size patchSize between the
% 'template' patch centred at row and col and a patch shifted by the
% current offset

% This time, use the integral image method!
% NOTE: Use the 'computeIntegralImage' function developed earlier to
% calculate your integral images
% NOTE: Use the 'evaluateIntegralImage' function to calculate patch sums

image = double(image);
[rows columns d] = size(image);

% The intialization is the same as the TemplateMatching in the naive case
delta = floor(patchSize/2);
delta_window = floor(searchWindowSize/2);

distances = zeros(1,searchWindowSize*searchWindowSize);
offsetsRows = zeros(1,searchWindowSize*searchWindowSize);
offsetsCols = zeros(1,searchWindowSize*searchWindowSize);
distances_index = 1;

% I store all my integral images in a Cell Array
integralImagesCell = cell(searchWindowSize, searchWindowSize);

% all the possible offsets are all the possible combinations of indices in
% thesearch window
    for row_searchWindow = -delta_window:delta_window
        for column_searchWindow = -delta_window:delta_window

            xOffset = row_searchWindow;
            yOffset = column_searchWindow;

            % HOW WE CAN SHIFT THE IMAGE
            % CASE ONE -- AFFINE TRANFSORMATION
            %T = maketform('affine', [1 0 0; 0 1 0; yOffset xOffset  1]);
            %shifted_image = imtransform(image, T,'XData',[1 size(image,2)],'YData',[1 size(image,1)]);

            % CASE TWO -- CIRCULAR SHIFT
            %shifted_image = circshift(image , [xOffset,yOffset]);

            %TRANSLATE THE IMAGE
            shifted_image = imtranslate(image,[yOffset, xOffset]);

            % Let's compute the integral image for the difference squared of the
            % two images (how I compute the integral image is explained in
            % the method "computeIntegralImage").
            integral_image = computeIntegralImage((double(shifted_image-image)).^2, false);

            % Store the result, note the cell has no negative indices
            integralImagesCell{xOffset+delta_window+1, yOffset+delta_window+1} = integral_image;

            %NOTE : we can also store the integral images ina matrix or in a
            %       dictionary
            %c([num2str(xOffset) ' ' num2str(yOffset)]) = integral_image;
            %integralImages(:,:, index:index+2) = integral_image;

        end
    end

    start_rows = max(row-delta_window, 1+delta);
    end_rows = min(row+delta_window, rows-delta);
    start_columns = max(col-delta_window, 1+delta);
    end_columns = min(col+delta_window, columns-delta);

            for row_searchWindow = start_rows : end_rows
                for column_searchWindow = start_columns : end_columns

                    % SOME DEBUG CODE
                    %disp('###################');
                    %disp(row_searchWindow-delta);
                    %disp(column_searchWindow);
                    %disp(delta);
                    %disp(column_searchWindow-delta);
                    %disp(column_searchWindow+delta);

                    xOffset = row_searchWindow - row;
                    yOffset = column_searchWindow - col;

                    % Get the corresponding integral image at the two
                    % offsets
                    integral_image = integralImagesCell{xOffset+delta_window+1, yOffset+delta_window+1};

                    % compute the distance between the two patches
                    distance = evaluateIntegralImage(integral_image, row_searchWindow, column_searchWindow, delta);

                    % Store the results
                    distances(distances_index) = distance;
                    offsetsRows(distances_index) = xOffset;
                    offsetsCols(distances_index) = yOffset;
                    distances_index = distances_index + 1;
                end
            end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function  A1 = UpdataArchive(A1,New,V,mu,NI)
% Update archive

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

% This function is written by Cheng He

    %% Delete duplicated solutions
    All       = [A1.decs;New.decs];
    [~,index] = unique(All,'rows');
    ALL       = [A1,New];
    Total     = ALL(index);

    %% Select NI solutions for updating the models
	if length(Total)>NI
        [~,active] = NoActive(New.objs,V);
        Vi         = V(setdiff(1:size(V,1),active),:);
        % Select the undeplicated solutions without re-evaluated solutions
        index = ismember(Total.decs,New.decs,'rows');
        Total = Total(~index);
        % Since the number of inactive reference vectors is smaller than
        % NI-mu, we cluster the solutions instead of reference vectors
        PopObj = Total.objs;
        PopObj = PopObj - repmat(min(PopObj,[],1),length(Total),1);
        Angle  = acos(1-pdist2(PopObj,Vi,'cosine'));
        [~,associate] = min(Angle,[],2);
        Via    = Vi(unique(associate)',:);
        Next   = zeros(1,NI-mu);
        if size(Via,1) > NI-mu
            [IDX,~] = kmeans(Via,NI-mu);
            for i = unique(IDX)'
                current = find(IDX==i);
                if length(current)>1
                    best = randi(length(current),1);
                else
                    best = 1;
                end
                Next(i)  = current(best);
            end
        else
            % Cluster solutions based on objective vectors when the number
            % of active reference vectors is smaller than NI-mu
            [IDX,~] = kmeans(Total.objs,NI-mu);
            for i   = unique(IDX)'
                current = find(IDX==i);
                if length(current)>1
                    best = randi(length(current),1);
                else
                    best = 1;
                end
                Next(i)  = current(best);
            end
        end
        A1 = [Total(Next),New];
    else
        A1 = Total;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%
% Description:
%
% The default behavior for an arbitrary function
% (ECGtask_arbitrary_function) of the strat function is to do nothing.
%
%
% See also default_finish_function, ECGtask_arbitrary_function
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 17/4/2015
% Last update: 17/4/2015
% Copyright 2008-2017
%
function bStarted = default_start_function(obj, ECG_header)

% Flag to indicate correct initialization.
bStarted = true;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% Author: James A. Douthwaite

classdef agent_2D_formation_HRVO < agent_2D_formation_RVO & agent_2D_HRVO
%% INITIALISE THE AGENT SPECIFIC PARAMETERS
    properties
        % PROPERTITIES UNIQUE TO AGENT CLASS
    end
    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function obj = agent_2D_formation_HRVO(varargin)
            % CALL THE SUPERCLASS CONSTRUCTOR
            obj = obj@agent_2D_formation_RVO(varargin);

            % //////////////// Check for user overrides ///////////////////
            [obj] = obj.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % NOTES:
        % This class simply overrides the 'process_timeCycle' function of
        % the superclass 'agent_2D_formation_RVO' such that it considers the
        % 'avoidanceCorrection' of the second super class 'agent_2D_HRVO'.
    end
end

repo: OpenMAS


text:
% @author: Maziar Raissi

function params_list = Navier_Stokes()
% quantile(params_list,[0.025 0.25 0.50 0.75 0.975])
clc; close all;

plt = 1;
save_plt = 0;

addpath ..
addpath ../Utilities
addpath ../Kernels/Navier_Stokes
addpath ../Utilities/export_fig

function CleanupFun()
    rmpath ../Utilities
    rmpath ../Kernels/Navier_Stokes
    rmpath ../Utilities/export_fig
end

finishup = onCleanup(@() CleanupFun());

rng('default')

set(0,'defaulttextinterpreter','latex')

%% Load Data
load('../Data/cylinder_fine.mat', 't_star', 'X_star', 'U_star', 'w_star')
N_star = size(X_star,1);
nsteps = size(t_star,1) - 1;

%% Setup
noise = 0.00;
U_data = U_star;
u_star = U_star(:,1,:);
U_data(:,1,:) = U_star(:,1,:) + noise*std(u_star(:))*randn(size(U_star(:,1,:)));
v_star = U_star(:,2,:);
U_data(:,2,:) = U_star(:,2,:) + noise*std(v_star(:))*randn(size(U_star(:,2,:)));

N0 = 251;
N1 = 249;
%% Optimize model
params_list = zeros(nsteps,2);
hyp = [log([1.0 1.0 1.0]) 0.0 -1.0 log([1.0 1.0 1.0]) -4.0];
idx1 = randsample(N_star, N0);
step = 1;
for i = 1:step:nsteps
    dt = t_star(i+step) - t_star(i);

    idx0 = idx1;
    X0 = X_star(idx0,:);
    U0 = U_data(idx0,:,i);

    idx1 = randsample(N_star,N1);
    X1 = X_star(idx1,:);
    U1 = U_data(idx1,:,i+step);

    model = HPM(X1, U1, X0, U0, dt, hyp);
    model = model.train(20);

    hyp = model.hyp;
    params_list(i,:) = [hyp(4) exp(hyp(5))];

    [pred_n_star, var_n_star] = model.predict(X_star);
    pred_U_star = reshape(pred_n_star,N_star,2);
    var_n_star = abs(diag(var_n_star));

    error = norm(pred_U_star(:,1) - U_star(:,1,i+step))/norm(U_star(:,1,i+step));

    fprintf(1,'=========================\n');
    fprintf(1,'Step: %d, Time = %.2f\n\nNLML = %.2f, Error = %.2e\n\n', i, ...
        t_star(i+step), model.NLML, error);

    str = sprintf('%.4f  ', params_list(i,:));
    fprintf('Parameters: %s\n\n', str)

    str = sprintf('%.4f  ', median(params_list(1:step:i,:),1));
    fprintf('Median: %s\n', str)
    fprintf(1,'=========================\n\n');

    if plt == 1
        if ~exist('fig','var')
            fig = figure(2);
        end
        set(fig,'units','normalized','outerposition',[0 0 1 1])
        clf

        subplot(2,2,1)
        plot_surface_griddata(X_star, U_star(:,1,i+step),'','','')

        subplot(2,2,2);
        plot_surface_griddata(X_star, U_star(:,2,i+step),'','','')

        subplot(2,2,3)
        plot_surface_griddata(X_star, pred_U_star(:,1),'','','')

        subplot(2,2,4);
        plot_surface_griddata(X_star, pred_U_star(:,2),'','','')

        drawnow()
    end

end

if save_plt == 1
    export_fig ./Figures/Navier_Stokes.png -r300
end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function res = vis_eigmodes(varargin)
%
% Visualize the results of an eigenmode decomposition of a VAR[p] model as
% computed using est_eigmode().
% This plots the characteristic frequencies and damping times of the
% eigenmodes, optionally sorted by dynamical importance [1].
%
% TODO: this function is incomplete. Plotting functionality not yet
% implemented in this function.
%
% Inputs:
%
%   EEG:        EEGLAB dataset containing CAT.EIGMODE substructure
%
% Optional:
%
%     TopPercentEigenmodes: Top percent eigenmodes (0-100)
%                           The percent of eigmodes with highest dynamical importance to plot
%                           Input Range  : [0  100]
%                           Default value: 100
%                           Input Data Type: real number (double)
%
%     ConvertPeriodToHz:    Convert res.period to Hertz
%                           Input Data Type: boolean
%
%     sortEigmodes:         Sort eigmodes by dynamical importance
%                           If true, sort for each window independently.If false, and topPercent < 100, sort each window,
%                           determine the topPercent most frequently occuring eigenmodes (mode across windows),and return these
%                           eigenmodes, unsorted across windows
%                           Input Data Type: boolean
%
%     plotResults:          Plot results
%                           Input Data Type: boolean
%
% Output:
%
%   res:        results structure
%
% See Also: est_eigenmode()
%
% References:
%
% [1] Neumaier A, Schneider T (2001) Estimation of parameters and
%   eigenmodes of multivariate autoregressive models. ACM Transactions on
%   Mathematical Software (TOMS) 27:57
%
% Author: Tim Mullen, 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    g = arg_define([0 1],varargin, ...
            arg_norep('EEG',mandatory), ...
            arg({'topPercent','TopPercentEigenmodes'},100,[0 100],'Top percent eigenmodes (0-100). The percent of eigmodes with highest dynamical importance to plot'), ...
            arg({'convertPeriodToHz','ConvertPeriodToHz'},true,[],'Convert res.period to Hertz.'), ...
            arg({'sortEigmodes'},true, [], ['Sort eigmodes by dynamical importance. ' ...
                                            'If true, sort for each window independently.' ...
                                            'If false, and topPercent < 100, sort each window, determine the topPercent most frequently occuring eigenmodes (mode across windows),' ...
                                            'and return these eigenmodes, unsorted across windows']), ...
            arg({'plotResults'},false,[],'Plot results') ...
            );

    % commit some variables to the workspace
    [data g] = hlp_splitstruct(g,{'EEG'});
    arg_toworkspace(data);
    clear data;

    [nch nmodes] = size(EEG.CAT.EIGMODE.modes{1});
    nwins        = length(EEG.CAT.EIGMODE.modes);

    % convert cell data to matrix format
    res.eigmodes    = cell2mat(EEG.CAT.EIGMODE.modes);
    res.eigmodes    = reshape(res.eigmodes,[nch nmodes nwins]);

    if ~isempty(EEG.CAT.EIGMODE.modeconf{1})
        res.eigmodesConf = cell2mat(EEG.CAT.EIGMODE.modeconf);
        res.eigmodesConf = reshape(res.eigmodesConf,[nmodes nwins]);
    else
        res.eigmodesConf = zeros(nmodes, nwins);
    end

    res.period = cell2mat(EEG.CAT.EIGMODE.period);
    res.period = reshape(res.period,[2 nmodes nwins]);
    res.periodConf = squeeze(res.period(2,:,:));
    res.period = squeeze(res.period(1,:,:));

    res.dampingTime = cell2mat(EEG.CAT.EIGMODE.dampingTime);
    res.dampingTime = reshape(res.dampingTime,[2 nmodes nwins]);
    res.dampingTimeConf = squeeze(res.dampingTime(2,:,:));
    res.dampingTime = squeeze(res.dampingTime(1,:,:));

    res.eigenvalues      = cell2mat(EEG.CAT.EIGMODE.lambda);

    % sort the eigmodes by dynamical importance
    res.exctn = cell2mat(EEG.CAT.EIGMODE.exctn);
    res.exctn = reshape(res.exctn,[nmodes nwins]);

    % sort eigmodes by dynamical importance
    if g.sortEigmodes
        [dummy res.topEigmodesIndex] = sort(res.exctn,1,'descend');
    elseif g.topPercent < 100
        % determine which modes occur most frequently across windows
        % and sort rows -- for all windows jointly -- in this order
        [dummy topmodes] = sort(res.exctn,1,'descend');
        topmodes = mode(topmodes,2);

        % remove duplicate topmodes
        for i=1:length(topmodes)
            val = topmodes(i);
            if ~isnan(val)
                topmodes(topmodes==topmodes(i))=NaN; % mark all duplicates
                topmodes(i)=val;
            end
        end
        topmodes(isnan(topmodes))=[];

%         [b m n] = unique_bc(topmodes);
%         mostFreqModes = topmodes(sort(m));
        res.topEigmodesIndex = repmat(topmodes,1,nwins);
    else
        res.topEigmodesIndex = repmat((1:nmodes)',1,nwins);
    end

    % extract the top percent of most important eigenmodes
    numTopEigmodes = round((g.topPercent/100)*size(res.topEigmodesIndex,1));
    res.topEigmodesIndex = res.topEigmodesIndex(1:numTopEigmodes,:);

    % extract the data corresponding to the top eigenmodes
    colshift = nmodes*(0:nwins-1);
    eigmodeLinearIndex = res.topEigmodesIndex + repmat(colshift,numTopEigmodes,1);
    res.exctn       = res.exctn(eigmodeLinearIndex);
    res.eigmodesConf = res.eigmodesConf(eigmodeLinearIndex);
    res.period       = res.period(eigmodeLinearIndex);
    res.periodConf   = res.periodConf(eigmodeLinearIndex);
    res.dampingTime  = res.dampingTime(eigmodeLinearIndex);
    res.dampingTimeConf = res.dampingTimeConf(eigmodeLinearIndex);
    res.eigenvalues     = res.eigenvalues(eigmodeLinearIndex);
    tmp = zeros(nch,numTopEigmodes,nwins);
    for ch=1:nch
        tmpeigmode  = squeeze(res.eigmodes(ch,:,:));
        tmp(ch,:,:) = tmpeigmode(eigmodeLinearIndex);
    end
    res.eigmodes = tmp;

    res.convertPeriodToHz = g.convertPeriodToHz;
    if g.convertPeriodToHz
        res.period = 1./res.period;
        res.periodConf = 1./res.periodConf;
    end

repo: REST


text:
function store = getStore(problem, x, storedb) %#ok<STOUT,INUSD>

    error('This file was removed from Manopt. Please use the StoreDB class.');

end

author: NicolasBoumal

repo: manopt


text:
function [ t_start, y_start ] = p17_start ( neqn )

%*****************************************************************************80
%
%% P17_START returns the starting point for problem p17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  delta = p17_param ( 'GET', 'DELTA', [] );
  y_start(1) = 1.0 - delta;
  y_start(2) = 0.0;
  y_start(3) = 0.0;
  y_start(4) = sqrt ( ( 1.0 + delta ) / ( 1.0 - delta ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function res = randomSol(model)
    res = randperm(model.city+model.veh-1);
end

author: lzane

repo: VRP-using-SA-with-Matlab


text:
function x = uplus( x )
    %UPLUS Unary plus.
    %   X = UPLUS(X) does nothing.
    %
    %   See also UMINUS.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

end


text:
% StackExchange Signal Processing Q76644
% https://dsp.stackexchange.com/questions/76644
% Simple and Effective Method to Estimate the Frequency of a Sine Signal in
% White Noise.
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     09/08/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

% Signal Parameters
numSamples      = 200;
samplingFreq    = 10; %<! The CRLB is for Normalized Frequency

% Sine Signal Parameters
expFreq    = 1; %<! [Hz]
expAmp     = 10; %<! High value to allow high SNR
expPhase   = pi * rand();

% Analysis Parameters
numRealizations = 500;
% SNR of the Analysis (dB)
vSnrdB = linspace(-10, 40, 100).';

%% Generate Data

angFreq = 2 * pi * (expFreq / samplingFreq);

vS = expAmp * exp(1i * ((angFreq * (0:(numSamples - 1))) + expPhase));
vS = vS(:);

numNoiseStd = length(vSnrdB);
vNoiseStd   = zeros(numNoiseStd, 1);

for ii = 1:numNoiseStd
    vNoiseStd(ii) = sqrt((expAmp * expAmp) / (10 ^ (vSnrdB(ii) / 10)));
end

tFreqErr = zeros(numRealizations, numNoiseStd, 2);

%% Analysis

for estType = 1:2
    for jj = 1:numNoiseStd
        noiseStd = vNoiseStd(jj);
        for ii = 1:numRealizations
            vW = (noiseStd / sqrt(2)) * (randn(numSamples, 1) + 1i * randn(numSamples, 1));
            vX = vS + vW;
            tFreqErr(ii, jj, estType) = expFreq - EstimateHarmonicFreqKay(vX, samplingFreq, estType);
        end
    end
end

mFreqErr = reshape(mean(tFreqErr .^ 2, 1), numNoiseStd, 2);

expMse     = expAmp * expAmp;
vNoiseVar  = vNoiseStd .^ 2;
vSnr       = expMse ./ vNoiseVar;

% See Steven Kay Estimation Theory (Pg. 57)
vFreqMseCrlb = (6 * samplingFreq * samplingFreq) ./ (((2 * pi) ^ 2) * vSnr * ((numSamples ^ 3) - numSamples));
vFreqMseCrlb = vFreqMseCrlb(:);
% In order to have any sampling rate multiply it by Fs ^ 2.

%% Display Results

figureIdx = figureIdx + 1;

hFigure = figure('Position', figPosLarge);
hAxes   = axes(hFigure);
hLineObj = plot(vSnrdB, 10 * log10([vFreqMseCrlb, mFreqErr]));
set(hLineObj, 'LineWidth', lineWidthNormal);
% set(hLineObj(1), 'LineStyle', 'none', 'Marker', '*');
% set(hLineObj(2), 'LineStyle', 'none', 'Marker', 'x');
set(get(hAxes, 'Title'), 'String', {['MSE of Harmonic Exponential Frequency Estimation']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['SNR [dB]']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['MSE']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['CRLB'], ['Kay Estimator Type 1'], ['Kay Estimator Type 2']});

if(generateFigures == ON)
    % saveas(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [azimuth, dist]=indirectRhumbProblem(latLonStart,latLonEnd,height,useHeightApprox,a,f,numSteps4Circ)
%INDIRECTRHUMBPROBLEM Given a starting and ending latitude and longitude on
%                     a reference ellipsoid, determine the heading (in
%                     radians East of North) and the distance one must
%                     travel on the shortest constant-heading course to go
%                     from the starting point to the stopping point. A
%                     constant heading course follows a rhumb line
%                     (loxodrome) and is usually not the shortest path
%                     between two points.
%
%INPUTS: latLonStart A 2X1 vector of the starting ellipsoidal latitude
%                    (North) and longitude (East) in radians. This cannot
%                    be a pole.
%          latLonEnd A 2X1 vector of the ending ellipsoidal latitude
%                    (North) and longitude (East) in radians.
%             height The height above the reference ellipsoid at which the
%                    trajectory should be determined. This changes the
%                    distance traveled, but not the azimuthal angle of
%                    departure. If this parameter is omitted, then the
%                    default value of 0 is used.
%    useHeightApprox If true, and the height is not zero, then an
%                    approximation is made for how dist scales with
%                    altitude. Specifically, an equiatorial trajectory will
%                    scale as (a+height)/a. Thus, this scaling factor is
%                    applied to any trajectory to scale dist with altitude.
%                    If height is false, a significantly slower iterative
%                    optimization technique is used. The default value is
%                    true. The difference made when useHeightApprox=false is
%                    can generally be assumed to be less than 80m. This
%                    parameter is ignored if height=0.
%                  a The semi-major axis of the reference ellipsoid. If
%                    this argument is omitted, the value in
%                    Constants.WGS84SemiMajorAxis is used.
%                  f The flattening factor of the reference ellipsoid. If
%                    this argument is omitted, the value in
%                    Constants.WGS84Flattening is used.
%      numSteps4Circ If height!=0 then an algorithm propagating a state
%                    in ECEF coordinates around the curved Earth is used to
%                    solve the direct geodetic problem as a step in solving
%                    the indirect geodetic problem. This parameter
%                    determines the number of steps that would be needed in
%                    the direct geodetic problem for a target that
%                    circumnavigates the globe around the equator. The
%                    default value if this parameter is not provided is
%                    2000. A value of 6000 appears to be about the best
%                    number for overall precision. Reducing the number of
%                    steps will speed up the function. This parameter is not
%                    used if height=0.
%
%OUTPUTS: azimuth The constant heading in radians East of North that one
%                 must travel to go on a constant-heading course from
%                 latLonStart to latLonEnd.
%            dist The distance that one must travel on a constant-heading
%                 course to go from latLonStart to latLonEnd.
%
%If height=0, the algorithm is mostly taken from [1]. However, a formula
%using isometric latitudes, which are described in Chapter 3 of [2] to get
%the azimuth angle was used, because it is simpler. The formula is also
%explicitly mentioned in Equation 3 of [3]. However, the expression for
%computing the distance from that paper is only for a sphere, not for an
%ellipsoid, which is why the Carlton-Wippern distance computation using an
%incomplete elliptic integral of the second kind is preferred.
%
%When the azimuth found by the technique is very close to +/-pi/2 (when one
%is traveling at nearly a constant latitude), the distance computation
%switches to assume that the latitude is indeed constant even if
%latLonStart(1)!=latLonStart(2). This avoid precision problems that arise
%as a very small number is multiplied by a very large number. However, this
%reduces the accuracy of the method.
%
%Generally, calling directRhumbProblem or directRhumbProbGen with the
%azimuth and dist returned by this function will return latLonStart.
%However, if the stopping point is at a pole, then directRhumbProblem will
%correctly return a polar location, but the longitude will generally be
%wrong.
%
%When height!=0, the algorithm can be significantly slower if no
%approximation is used. Around the equator, the distance scales as
%(a+height)/a*dist as one changes the ellipsoidal height. This scaling
%applied to dist in  non-equatorial trajectories is the approximation if
%useHeightApprox=true. When useHeightApprox=false, the approximate value
%is used to determine the bounds around which the fminbnd function
%searches. The maximum error in the approximation is expected to be less
%than 80m. The search region for the value of dist used in the fminbnd
%function was set to 0.9*dist to 1.1*dist, where dist is the distance
%obtained after scaling the distance from the zero-altitude solution.
%
%EXAMPLE:
%A trajectory that crosses the international date line and and goes from
%the Northern hemisphere to the southern hemisphere. We also compute the
%reverse path and show that the azimuth angles in each direction are
%consistent with each other. We then plot the trajectory on an image of the
%spherical Earth.
% N=100;
% latStart=degMinSec2Rad(37,47.5);
% lonStart=degMinSec2Rad(-122,-27.8);
% latEnd=degMinSec2Rad(-33,-51.7);
% lonEnd=degMinSec2Rad(151,12.7);
%
% latLonStart=[latStart;lonStart];
% latLonEnd=[latEnd;lonEnd];
% [azimuth,dist]=indirectRhumbProblem(latLonStart,latLonEnd);
% [azEnd,distRev]=indirectRhumbProblem(latLonEnd,latLonStart);
%
% %If the forward and reverse estimates agree, then these values will
% %ideally be zero.
% azimuth-(-azEnd)
% dist-distRev
%
% distVals=linspace(0,dist,N);
% latLonWayPoints=directRhumbProblem(latLonStart,azimuth,distVals);
%
% %Show that the approximate direct algorithm reaches nearly the same
% %endpoint as the indirect algorithm.
% xEndWay=ellips2Cart([latLonWayPoints(:,end);0]);
% xEnd=ellips2Cart([latLonEnd;0]);
% max(abs(xEndWay-xEnd))
% max(abs(wrapRange(latLonWayPoints(:,end)-latLonEnd,-pi,pi)))
%
% xStartCart=ellips2Cart([latLonStart;0]);
% xEndCart=ellips2Cart([latLonEnd;0]);
% %The path is displayed slightly above the Earth's surface to make it
% %easier to see.
% pathPoints=ellips2Cart([latLonWayPoints;0.02*ones(1,N)]);
%
% figure(1)
% clf
% hold on
% plotMapOnEllipsoid([]);
% scatter3(xStartCart(1),xStartCart(2),xStartCart(3),100,'filled')
% scatter3(xEndCart(1),xEndCart(2),xEndCart(3),100,'filled')
% plot3(pathPoints(1,:),pathPoints(2,:),pathPoints(3,:),'-r','linewidth',4)
%
%REFERENCES:
%[1] K. C. Carlton-Wippern, "On loxodromic navigation," Journal of
%    Navigation, vol. 45, no. 2, pp. 292-297, May 1992.
%[2] J. P. Snyder, "Map projections- a working manual," U.S. Geological
%    Survey, Tech. Rep. 1395, 1987.
%[3] J. Alexander, "Loxodromes: A rhumb way to go," Mathematics Magazine,
%    vol. 77, no. 5, pp. 349-356, Dec. 2004.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<7||isempty(numSteps4Circ))
   numSteps4Circ=2000;
end

if(nargin<6||isempty(f))
    f=Constants.WGS84Flattening;
end

if(nargin<5||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

if(nargin<4||isempty(useHeightApprox))
   useHeightApprox=true;
end

if(nargin<3||isempty(height))
   height=0;
end

%Extract the components
latStart=latLonStart(1);
lonStart=latLonStart(2);
latEnd=latLonEnd(1);
lonEnd=latLonEnd(2);

%The first numerical eccentricity of the ellipsoid.
e=sqrt(2*f-f^2);

%Convert the ellipsoidal latitudes to reduced co-latitudes. A co-latitude
%is pi/2 minus the latitude. Also,
nu1=pi/2-ellipsLat2ReducedLat(latStart,f);
nu2=pi/2-ellipsLat2ReducedLat(latEnd,f);

%Though not mentioned in the above papers, the difference in the longitudes
%must be wrapped to the range of -pi/pi or else one will get useless
%results when crossing the -pi/pi boundary.
num=wrapRange(lonEnd-lonStart,-pi,pi,false);
%Equation 11 in the paper provides an expression to get the azimuth.
%however, it is simpler if one just uses isometric latitudes.
val1=ellipsLat2IsoLat(latStart,f);
val2=ellipsLat2IsoLat(latEnd,f);

if(~isfinite(val1))
    warning('The starting point is located at a geographic pole. Azimuth values will be inaccurate.')
end

if(isfinite(val1)||isfinite(val2))
    %This will properly return 0 and pi for infinite values of val1 when
    %val2 is finite and vice versa, which corresponds to headings to or
    %from a pole.
    denom=val2-val1;
    azimuth=atan2(num,denom);
else%If neither is finite, then that is the case where one is going from
    %pole to pole. In such an instance, just set the heading to 0, if going
    %North, and pi, if going South.

    if(latStart>latEnd)
        azimuth=pi;
    else
        azimuth=0;
    end
end

%The distance
if(abs(abs(azimuth)-pi/2)>2e-8)
    %Equation 12 in the paper.
    dist=a*abs(sec(azimuth))*abs(ellipIntInc2Kind(nu2,e^2)-ellipIntInc2Kind(nu1,e^2));
else
    %Equation 14b in the paper.
    dist=a*abs(sin(nu1))*abs(wrapRange(lonEnd-lonStart,-pi,pi));
end

%If a non-zero height is given, then iterate over the direct rhumb
%problem at height to determine the solution.
if(height~=0)
    endCart=ellips2Cart([latLonEnd;height],a,f);

    %The approximate scaling for the height.
    dist=((a+height)/a)*dist;

    %If a computationally-intensive but more precise algorithm to search
    %for the true distance at altitude should be used instead of a
    %simple approximation of scaling the distance.
    if(useHeightApprox==false)
        %Assume that the correct height-adjusted distance is within 10% of the
        %scaled distance value.
        distFun=@(distCur)distCostFunc(distCur,endCart,latLonStart,azimuth,height,a,f,numSteps4Circ);
        dist=fminbnd(distFun,0.9*dist,1.1*dist);
    end
end
end

function cost=distCostFunc(distCur,endCart,latLonStart,azStart,height,a,f,numSteps4Circ)
    latLonCalc=directRhumbProbGen(latLonStart,azStart,distCur,height,false,a,f,numSteps4Circ);
    cost=norm(ellips2Cart([latLonCalc;height],a,f)-endCart);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function value = angle_rad_2d ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_RAD_2D returns the angle swept out between two rays in 2D.
%
%  Discussion:
%
%    Except for the zero angle case, it should be true that
%
%      ANGLE_RAD_2D(P1,P2,P3) + ANGLE_RAD_2D(P3,P2,P1) = 2 * PI
%
%        P1
%        /
%       /
%      /
%     /
%    P2--------->P3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2), P2(2), P3(2), define the rays
%    P1 - P2 and P3 - P2 which in turn define the angle.
%
%    Output, real VALUE, the angle swept out by the rays, measured
%    in radians.  0 <= VALUE < 2*PI.  If either ray has zero length,
%    then VALUE is set to 0.
%
  p(1) = ( p3(1) - p2(1) ) * ( p1(1) - p2(1) ) ...
       + ( p3(2) - p2(2) ) * ( p1(2) - p2(2) );

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

  if ( p(1) == 0.0 & p(2) == 0.0 )
    value = 0.0;
    return
  end

  value = atan2 ( p(2), p(1) );

  if ( value < 0.0 )
    value = value + 2.0 * pi;
  end

  return
end