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

%%
timer_proxA = tic;

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

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

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

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

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

% initialize cells/vectors
for j=1:nCha
    FTb{1,j} = iFFT3D(b{1,j});
    y{1,j} = real(FTb{1,j}); % needed for y_old
    y{1,j+nCha} = imag(FTb{1,j});
end;
z = FTb; % starting point

x_wave = cell(1,2*nCha);
x_helper = x_wave;
x_wave_helper = x_wave;
x_tv = x_wave;
x_nltv = x_wave;
for j=1:2*nCha
    x_nltv{1,j} = 0;
end;
x_g = x_wave;
x_g_helper = x_wave;
x_g_proxA = x_wave;
z_comp = x_wave;

%% MAD dependent lambdas:
flag_MAD = true;
if flag_MAD
    x_wavedec_2D = cell(1,nCha);
    x_wavedec_3D = cell(1,nCha);
    threshold_2D = zeros(nCha,1); % one threshold value is enough (only for testing purposes)
    threshold_3D = zeros(nCha,1); % one threshold value is enough (only for testing purposes)
    if flagRealAndImag
        for j=1:nCha
            for l = 1:nSlices
                x_wavedec_2D{1,j} = wavedec2(abs(z{1,j}),waveletStages,waveletFilterName_l12); % atm only based on l1-daubechie
                x_wave_fine_scale_2D = size(x_wavedec_2D{1,j},2) - (3*waveS_l12(waveletStages+1,1)*waveS_l12(waveletStages+1,2));
                threshold_2D(j) = mad(x_wavedec_2D{1,j}(x_wave_fine_scale_2D:end),1);
            end;
            x_wavedec_3D{1,j} = wavedec3(abs(z{1,j}),waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
            x_wave_fine_scale_3D = zeros(0,0);
            for k = (7*(x_wavedec_3D{1,j}.level-1)+2):maxWavecells
                x_wave_fine_scale_3D = [x_wave_fine_scale_3D  x_wavedec_3D{1,j}.dec{k}];
            end;
            [size1, size2, size3] = size(x_wave_fine_scale_3D);
            threshold_3D(j) = mad(reshape(x_wave_fine_scale_3D,size1*size2*size3,1));
        end;
    else
        for j=1:nCha
            for l = 1:nSlices
                x_wavedec_2D{1,j} = wavedec2(real(z{1,j}),waveletStages,waveletFilterName_l12); % atm only based on l1-daubechie
                x_wavedec_2D{1,j+nCha} = wavedec2(imag(z{1,j}),waveletStages,waveletFilterName_l12); % atm only based on l1-daubechie
                x_wave_fine_scale_2D_real = size(x_wavedec_2D{1,j},2) - (3*waveS_l12(waveletStages+1,1)*waveS_l12(waveletStages+1,2));
                threshold_2D(j) = mad(x_wavedec_2D{1,j}(x_wave_fine_scale_2D_real:end),1);
                x_wave_fine_scale_2D_imag = size(x_wavedec_2D{1,j},2) - (3*waveS_l12(waveletStages+1,1)*waveS_l12(waveletStages+1,2));
                threshold_2D(j+nCha) = mad(x_wavedec_2D{1,j}(x_wave_fine_scale_2D_imag:end),1);
            end;
            x_wavedec_3D{1,j} = wavedec3(real(z{1,j}),waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
            x_wavedec_3D{1,j} = wavedec3(imag(z{1,j}),waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
            x_wave_fine_scale_3D_real = zeros(0,0);
            x_wave_fine_scale_3D_imag = zeros(0,0);
            for k = (7*(x_wavedec_3D{1,j}.level-1)+2):maxWavecells
                x_wave_fine_scale_3D_real = [x_wave_fine_scale_3D_real  x_wavedec_3D{1,j}.dec{k+2}];
                x_wave_fine_scale_3D_imag = [x_wave_fine_scale_3D_imag  x_wavedec_3D{1,j}.dec{k+2}];
            end;
            [size1, size2, size3] = size(x_wave_fine_scale_3D_real);
            threshold_3D(j) = mad(reshape(x_wave_fine_scale_3D_real,size1*size2*size3,1));
            threshold_3D(j+nCha) = mad(reshape(x_wave_fine_scale_3D_imag,size1*size2*size3,1));
        end;
    end;
    clear x_wavedec
else
    threshold_2D(1:nCha) = 1;
    threshold_3D(1:nCha) = 1;
end;

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

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

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

    t_new = (1+sqrt(1+4*t_old^2))/2;
    t_old = t_new;

    y_old = y;   % y_old = y for complex case

    %% landweber step
    for j = 1:nCha
        x_helper{1,j} = iFFT3D(FFT3D_mask(z{1,j},mask)) -FTb{1,j} ;
        z{1,j} = z{1,j} - x_helper{1,j}/L;
        z_comp{1,j} = real(z{1,j});
        z_comp{1,j+nCha} = imag(z{1,j});
    end;

    %% l1-Wavelet
    if flagWave
        for j = 1:2*nCha
            x_wave_helper{1,j} = wavedec3(z_comp{1,j},waveletStages,waveletFilterName_l1);
        end;
        if flagRealAndImag
            for j = 1:nCha
                for k = 1:maxWavecells
                    dec_cell{j}{k} = x_wave_helper{1,j}.dec{k,1};
                    dec_cell{j+nCha}{k} = x_wave_helper{1,j+nCha}.dec{k,1};
                    [c1, c2, c3] = size(dec_cell{j}{k});
                    [dec_cell{j}{k}, dec_cell{j+nCha}{k}] = groupthresh_2vec(reshape(dec_cell{j}{k},c1*c2*c3,1,1),reshape(dec_cell{j+nCha}{k},c1*c2*c3,1,1),threshold_wave(j));
                    x_wave_helper{1,j}.dec{k,1} = reshape(dec_cell{j}{k},c1,c2,c3);
                    x_wave_helper{1,j+nCha}.dec{k,1} = reshape(dec_cell{j+nCha}{k},c1,c2,c3);
                    dec_cell{j}{k} = [];
                    dec_cell{j+nCha}{k} = [];
                end;
            end;
        else
            for j = 1:2*nCha
                x_wave_helper{1,j} = softthresh_real(x_wave_helper{1,j},threshold_wave(j));
                for k = 1:maxWavecells
                    dec_cell{j}{k} = x_wave_helper{1,j}.dec{k,1};
                    [c1, c2, c3] = size(dec_cell{j}{k});
                    dec_cell{j}{k} = softthresh_real(reshape(dec_cell{j}{k},c1*c2*c3,1,1),threshold_wave(j));
                    x_wave_helper{1,j}.dec{k,1} = reshape(dec_cell{j}{k},c1,c2,c3);
                    dec_cell{j}{k} = [];
                end;
            end;
        end;
        for j = 1:2*nCha
            x_wave{1,j} = waverec3(x_wave_helper{1,j});
        end;
    end;

    %% TV
    if flagTV
        if ~flagTV_iso

            for j = 1:2*nCha
                for i = 1:nSlices
                    x_tv{1,j}(:,:,i) = MTV_2D(z_comp{1,j}(:,:,i),threshold_TV(j),n1,n2);
                end;
            end;
        else

            for j = 1:2*nCha
                for i = 1:nSlices
                    if (itr==1)
                        [x_tv{1,j}(:,:,i), P]=denoise_TV_One((z_comp{1,j}(:,:,i)), threshold_TV(j),-inf,inf,[],parsin);
                    else
                        [x_tv{1,j}(:,:,i), P]=denoise_TV_One((z_comp{1,j}(:,:,i)), threshold_TV(j),-inf,inf,P,parsin);
                    end;
                end;
            end;
        end;
    end;

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

    %% l12-Wavelet
    if flagGroup
        for j = 1:2*nCha
            for i = 1:nSlices
                if flag_extendImage
                    z_proxA{1,j}(:,:,i) = extend_image(z_comp{1,j}(:,:,i), waveS_l12proxA, waveletStages, proxA_extend_y, proxA_extend_x);
                    x_g_helper{i}{1,j} = wavedec2(z_proxA{1,j}(:,:,i),waveletStages,waveletFilterName_l12);
                else
                    x_g_helper{i}{1,j} = wavedec2(z_comp{1,j}(:,:,i),waveletStages,waveletFilterName_l12);
                end;

                if flag_wavetree
                    x_g_helper{i}{2,j} = (G_prox*x_g_helper{i}{1,j}')';
                else
                    x_g_helper{i}{2,j} = zeros(1,size(x_g_helper{i}{1,j}(:,:,i),2));
                end;
            end;
        end;
        for i = 1:nSlices
            if flagRealAndImag
                x_g_helper{i}(1,:) = softthresh_proxA_cha(x_g_helper{i}(:,:),threshold_group,2*nCha,waveS_l12proxA,groupnorm_index);
            else
                x_g_helper{i}(1,1:nCha) = softthresh_proxA_cha(x_g_helper{i}(:,1:nCha),threshold_group_real,nCha,waveS_l12proxA,groupnorm_index);
                x_g_helper{i}(1,nCha+1:2*nCha) = softthresh_proxA_cha(x_g_helper{i}(:,1:nCha),threshold_group_imag,nCha,waveS_l12proxA,groupnorm_index);
            end;
        end;
        for j = 1:2*nCha
            for i = 1:nSlices
                x_g_proxA{1,j}(:,:,i) = waverec2(x_g_helper{i}{1,j},waveS_l12proxA,waveletFilterName_l12);
                x_g{1,j}(:,:,i) =  x_g_proxA{1,j}(1:end-proxA_extend_y,1:end-proxA_extend_x,i);
            end;
        end;
    end;

    %% add prox(.)
    for j = 1:2*nCha
        y{1,j} = zeros(n1,n2,nSlices);
        if flagWave y{1,j} = y{1,j} + x_wave{1,j}.*regularizerWeights(1); end;
        if flagTV y{1,j} = y{1,j} + x_tv{1,j}.*regularizerWeights(2); end;
        if flagGroup y{1,j} = y{1,j} + x_g{1,j}.*regularizerWeights(3); end;
        if flagNLTV y{1,j} = y{1,j} + x_nltv{1,j}.*regularizerWeights(4); end;

        if itr < itrNLTV
            y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagWave.*regularizerWeights(1) + flagGroup.*regularizerWeights(3));
        else
            y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagNLTV.*regularizerWeights(4) + flagWave.*regularizerWeights(1) + flagGroup.*regularizerWeights(3));
        end;

        if flag_fast
            y{1,j}=y{1,j}+((t_old-1)/t_new).*(y{1,j}-y_old{1,j});
        end;
    end;

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

    metrics.xtime(itr+1)= toc(timer_proxA);
    for j = 1:nCha
        z{1,j} = y{1,j} + 1i*y{1,j+nCha};
    end;
%     [metrics, ssim_map{1,2}] = get_metrics_itr( im_ref, im_ref_full, z, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma, nSlices );

end;
dispProgress('Proximal Average', 'Close');

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

end

author: thomaskuestner


text:
function geometry_test001 ( )

%*****************************************************************************80
%
%% TEST001 tests ANGLE_CONTAINS_POINT_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 6;
  n_angle = 12;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST001\n' );
  fprintf ( 1, '  ANGLE_CONTAINS_POINT_2D sees if a point\n' );
  fprintf ( 1, '  lies within an angle.\n' );
  fprintf ( 1, '\n' );
%
%  An acute angle (45 degrees)
%
  for j = 1 : ntest

    if ( j == 1 )

      p1(1:2,1) = [ 1.0; 0.0 ];
      p2(1:2,1) = [ 0.0; 0.0 ];
      p3(1:2,1) = [ 1.0; 1.0 ];

    elseif ( j == 2 )

      p1(1:2,1) = [ 1.0; 0.0 ];
      p2(1:2,1) = [ 0.0; 0.0 ];
      p3(1:2,1) = [ 0.0; 1.0 ];

    elseif ( j == 3 )

      p1(1:2,1) = [ 1.0; -1.0 ];
      p2(1:2,1) = [ 0.0;  0.0 ];
      p3(1:2,1) = [ 0.0;  1.0 ];

    elseif ( j == 4 )

      p1(1:2,1) = [  1.0; 0.0 ];
      p2(1:2,1) = [  0.0; 0.0 ];
      p3(1:2,1) = [ -1.0; 0.0 ];

    elseif ( j == 5 )

      p1(1:2,1) = [ 1.0;  0.0 ];
      p2(1:2,1) = [ 0.0;  0.0 ];
      p3(1:2,1) = [ 0.0; -1.0 ];

    elseif ( j == 6 )

      p1(1:2,1) = [ 1.0;  0.0 ];
      p2(1:2,1) = [ 0.0;  0.0 ];
      p3(1:2,1) = [ 1.0; -0.01 ];

    end

    r8vec_print ( 2, p1, '  Vertex P1' )
    r8vec_print ( 2, p2, '  Vertex P2' )
    r8vec_print ( 2, p3, '  Vertex P3' )

    fprintf ( 1, '\n' );
    fprintf ( 1, '       X            Y       Inside?\n' );
    fprintf ( 1, '\n' );

    for i = 0 : n_angle

      thetar = i * 2.0 * pi / n_angle;
%
%  For some bizarre MATLAB peculiarity, I can't say
%    p(1:2) = [ cos ( thetar ), sin ( thetar ) ];
%  MATLAB (Student 6.5) complains that COS has the wrong number
%  of arguments.  But write it out like this, and magically, everything
%  is fine.  People think I'm crabby for no reason...
%
      p(1,1) = cos ( thetar );
      p(2,1) = sin ( thetar );

      inside = angle_contains_point_2d ( p1, p2, p3, p );

      fprintf ( 1, '  %12f  %12f  %1d\n', p(1:2,1), inside );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% set values, load data
nx = 144;
ny = 256;
nc = 8;
niter = 250;
mu = 2^-3;
nlevels = 3;
beta = 2^8.5;

load braindat;
load kmask;
load brainxinf;
% load braindb4xinf;
mask = true(nx,ny);

% mask = abs(xinf) > 0.1*max(col(abs(xinf)));

% generate and sample data
A = (1/sqrt(nx*ny))*Gdft('mask', true(nx,ny), 'fftshift', 1, ...
    'ifftshift', 1);

dat = A*reshape(coilim, [nx*ny nc]); clear coilim;

kmask = kmask((256-nx)/2+1:256-(256-nx)/2,:);
dat = col(dat(col(kmask),:));

n = kmask; b = 1:nx*ny;
n = b(n); clear b;

% build system matrices
A = Apsm('knownfn', 'time', 'v', 1, 'n', n(:), 'smap', smap, 'immask', ...
    true(nx,ny), 'nk', nx*ny);
W = Godwt1(mask, 'level', nlevels, 'wname', 'haar');
P = Gdiag(1./(col(sum(abs(smap).^2,3)) + mu));

% don't regularize approx coeffs
beta = beta*ones(nx,ny);
beta(1:nx/2^nlevels,1:ny/2^nlevels) = 0;
beta = col(beta);

% initializations
x = A'*dat; x = x./ col(sum(abs(smap).^2,3));
x = x(mask);
z = W*x;
eta = zeros(size(z));

% system matrix
A = Apsm('knownfn', 'time', 'v', 1, 'n', n(:), 'smap', smap, 'immask', ...
    mask, 'nk', nx*ny);

% algorithm book-keeping
shrink = @(t, a) (t-a .* sign(t)) .* (abs(t) > a);
i = 0;
thetime(1) = 0;
xdist(1) = norm(col(x) - col(xinf(mask)));

% go
tic;
while i < niter
    printm('iteration %d of %d', i+1, niter);

    z = shrink(W*x + eta, beta./mu);
    x = qpwls_pcg1(x, [A; sqrt(mu)*W], 1, [dat(:); sqrt(mu)*(z - eta)], ...
        0, 'niter', 5, 'precon', P);
    eta = eta - (z - W*x);

    thetime(i+2) = toc;
    xdist(i+2) = norm(col(x) - col(xinf(mask)));
    i = i+1;
end
toc;
x = embed(x, mask);
im(x)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function gqn_sparse_test ( )

%*****************************************************************************80
%
%% GQN_SPARSE_TEST uses the GQN function to build a sparse grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2014
%
%  Author:
%
%    Original MATLAB version by Florian Heiss, Viktor Winschel.
%    This MATLAB version by John Burkardt.
%
%  Local parameters:
%
%    Local, integer D, the spatial dimension.
%
  d = 5;
  trueval = hermite_integral ( d );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GQN_SPARSE_TEST:\n' );
  fprintf ( 1, '  GQN sparse grid:\n' );
  fprintf ( 1, '  Gauss quadrature, Hermite weight over (-oo,+oo).\n' );
  fprintf ( 1, '  Exact integral is %g\n', trueval );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   D  Level   Nodes    SG estimate MC estimate\n' );
  fprintf ( 1, '                       SG error    MC error\n' );
  fprintf ( 1, '\n' );

  for k = 1 : 8
%
%  Compute sparse grid estimate.
%
    [ x, w ] = nwspgr ( 'gqn', d, k );
    [ n, ~ ] = size ( x );
    fx = hermite_integrand ( d, n, x );
    sgappr = w' * fx;
    sgerror = sqrt ( ( sgappr - trueval ) .^ 2 ) / trueval;
%
%  Average 1000 Monte Carlo estimates.
%
    sim = zeros(1000,1);
    for r = 1 : 1000
      x = randn ( n, d );
      fx = hermite_integrand ( d, n, x );
      sim(r) = mean ( fx );
    end
    simappr = mean ( sim );
    simerror = sqrt ( mean ( ( sim - trueval ) .^ 2 ) ) / trueval;

    fprintf ( '  %2d     %2d  %6d  %10.5g  %10.5g\n', ...
      d, k, n, sgappr, simappr );
    fprintf ( '                     %10.5g  %10.5g\n', ...
               sgerror, simerror );

  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment2B_K20_N30
close all;
global parameters

%STOMP PARAMETERS
%number of particles
K = 20;

%repeat experiment number of times
parameters.n_repeat = 30;
parameters.experiment_name = 'experiment2B_K20_N30.mat';
parameters.animate = 0;

%repeat the experiment E times
random_manips=[];
Gout = [];
for i=1:parameters.n_repeat
    close all
    [pk, final_manip] = path_planning_SCO_4DOF(K);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    save(parameters.experiment_name)
end

'ended'
parameters.experiment_name

repo: ARTE


text:
function pde = quadCurlDataSmooth1
% (curl)^4 u = f on \Omega = [0,pi]^3
% u x n = 0
% (curl u) x n = 0
% div u = g
%
pde = struct('exactu', @exactu,'g_D', @g_D, 'g', @divu, ...
             'curlu', @curlu, 'curlcurlu',@curlcurlu, ...
             'tricurlu',@tricurlu,'quadcurlu',@quadcurlu, 'pentacurlu',@pentacurlu);

Scaling = 10; % scaling

    function s = exactu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s(:,1) = 0*x;
        s(:,2) = 0*y;
        s(:,3) = (sin(x)).^2.*(sin(y)).^2.*sin(z);
        s = s/Scaling;
    end

    function s = curlu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s(:,1) =  sin(2*y).*(sin(x)).^2.*sin(z);
        s(:,2) = -sin(2*x).*(sin(y)).^2.*sin(z);
        s(:,3) = 0*z;
        s = s/Scaling;
    end

    function s = curlcurlu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s(:,1) = sin(2*x).*(sin(y)).^2.*cos(z);
        s(:,2) = (sin(x)).^2.*sin(2*y).*cos(z);
        s(:,3) = 4*(sin(x)).^2.*(sin(y)).^2.*sin(z) ...
               - 2*(cos(y)).^2.*(sin(x)).^2.*sin(z) ...
               - 2*(cos(x)).^2.*(sin(y)).^2.*sin(z);
        s = s/Scaling;
    end

    function s = tricurlu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s(:,1) = 7*(sin(x)).^2.*sin(2*y).*sin(z) - 2*(cos(x)).^2.*sin(2*y).*sin(z);
        s(:,2) = 2*(cos(y)).^2.*sin(2*x).*sin(z) - 7*sin(2*x).*(sin(y)).^2.*sin(z);
        s(:,3) = 0*z;
        s = s/Scaling;
    end

    function s = quadcurlu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s(:,1) = 7*cos(z).*sin(2*x).*(sin(y)).^2 - 2*(cos(y)).^2.*cos(z).*sin(2*x);
        s(:,2) = 7*cos(z).*(sin(x)).^2.*sin(2*y) - 2*(cos(x)).^2.*cos(z).*sin(2*y);
        s(:,3) = 8*(cos(x)).^2.*(cos(y)).^2.*sin(z) - 18*(cos(x)).^2.*(sin(y)).^2.*sin(z) ...
              - 18*(cos(y)).^2.*(sin(x)).^2.*sin(z) + 28*(sin(x)).^2.*(sin(y)).^2.*sin(z);
        s = s/Scaling;
    end

    function s = pentacurlu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s(:,1) = 53*(sin(x)).^2.*sin(2*y).*sin(z) - 28*(cos(x)).^2.*sin(2*y).*sin(z);
        s(:,2) = 28*(cos(y)).^2.*sin(2*x).*sin(z) - 53*sin(2*x).*(sin(y)).^2.*sin(z);
        s(:,3) = 0*z;
        s = s/Scaling;
    end

    function s = g_D(p) % Dirichlet boundary condition for u
        s = exactu(p);
    end

    function s = divu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s = (sin(x)).^2.*(sin(y)).^2.*cos(z);
        s = s/Scaling;
    end
end

repo: ifem


text:
% Load data "mu10.mat".
load mu10;

% Create a figure.
figure('DefaultAxesFontSize',16);
ax = gca;

% Plot the mu curves.
xbound = [1,1,1:500,500,500,fliplr(1:500)];
ybound = [mean(mu(:,1))+std(mu(:,1)),mean(mu(:,1))-std(mu(:,1)),mean(mu)-std(mu),...
    mean(mu(:,end))-std(mu(:,end)),mean(mu(:,end))+std(mu(:,end)),fliplr(mean(mu)+std(mu))];
hold on;
plot(mu','-','Color',[0.76,0.87,0.78]);
hold on;
s = fill(xbound,ybound,'y','EdgeColor','white');
alpha(s,0.20);
hold on;
h = plot(mean(mu),'Color','r','LineWidth',2);
hold on;
p = plot(39.01*ones(500,1),'-.','Color','magenta','LineWidth',2);

% Change settings and items.
box on;
grid on;
ax.XTick = [0,50,100,200,500];
ax.XTickLabel = [0,50,100,200,500];
ylim([38,40.5]);
ax.YTick = [38,38.5,38.75,39,39.25,39.50,39.75,40.5];
ax.YTickLabel = ['38.00';'38.50';'38.75';'39.00';'39.25';'39.50';'39.75';'40.50'];
ax.GridLineStyle = '-.';
xlabel('Iteration');
ylabel('Speed (km/h)');
legend([p,h],{'Mean of Observations', 'Model Parameter \mu'});

% Save the figure as "mu_curve10.pdf".
set(gcf, 'PaperSize', [6 4.5]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 6 4.5]);
saveas(gcf,'mu_curve10','pdf');

author: xinychen

repo: academic-drawing


text:
function res = killer(n,p)

persistent plist rd
plist = [plist p];

switch n
    case 1
        res = 'cooperate';
        plist = p;
        rd = 0;
    case 20
        res = 'defect';
    otherwise
        if any( (plist(n-1:n)==1) + (plist(n-1:n)==5) ) && rd<5
            res = 'defect';
            rd = rd+1;
        else
            res = 'cooperate';
            if ~rd<5, rd=0; end
        end
end

author: Sable

repo: mcbench-benchmarks


text:
% OP_F_V_TP: assemble the right-hand side vector r = [r(i)], with  r(i) = (f, v_i), exploiting the tensor product structure.
%
%   rhs = op_f_v_tp (spv, msh, coeff);
%
% INPUT:
%
%   spv:   object representing the function space (see sp_vector)
%   msh:   object defining the domain partition and the quadrature rule (see msh_cartesian)
%   coeff: function handle to compute the source function
%
% OUTPUT:
%
%   rhs: assembled right-hand side
%
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function rhs = op_f_v_tp (space, msh, coeff)

  for icomp = 1:space.ncomp_param
    for idim = 1:msh.ndim
      size1 = size (space.scalar_spaces{icomp}.sp_univ(idim).connectivity);
      if (size1(2) ~= msh.nel_dir(idim))
        error ('The discrete space is not associated to the mesh')
      end
    end
  end

  rhs = zeros (space.ndof, 1);

  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    sp_col  = sp_evaluate_col (space, msh_col);

    for idim = 1:msh.rdim
      x{idim} = reshape (msh_col.geo_map(idim,:,:), msh_col.nqn, msh_col.nel);
    end

    rhs = rhs + op_f_v (sp_col, msh_col, coeff (x{:}));
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function p02_header ( )

%*****************************************************************************80
%
%% P02_HEADER prints some information about problem 02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    None.
%
  m = 2;

  center = [ 0.0, 0.0 ];
  r1 = 1.0;
  r2 = 0.4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'P02:\n' );
  fprintf ( 1, '  Strang and Persson example #2\n' );
  fprintf ( 1, '  The unit circle, with a concentric hole.\n' );
  fprintf ( 1, '  Radius1 = %f\n', r1 );
  fprintf ( 1, '  Center = (%f, %f)\n', center(1:2) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A uniform mesh density is requested.\n' );
  fprintf ( 1, '  Element sizes tried were 0.4, 0.2, 0.1.\n' );

  fprintf ( 1, '\n' );
  boundary_segment_num = p02_boundary_segment_num ( );
  fprintf ( 1, '  Number of boundary segments = %d\n', boundary_segment_num );
  fixed_num = p02_fixed_num ( );
  fprintf ( 1, '  Number of fixed points = %d\n', fixed_num );
  hole_num = p02_hole_num ( );
  fprintf ( 1, '  Number of holes = %d\n', hole_num );

  return
end


text:
function computeAssociations_Textures(pathUnivariate,pathFeatures,cohorts,outcomes,featType,textType,textName,nPatient,scale_mat,algo_cell,Ng_mat)

startpath = pwd;
cd(pathUnivariate)
fid = fopen('Texture_SignificanceProportion.txt','w');

nameOutcomes = fieldnames(outcomes.(cohorts{1})); nOutcomes = numel(nameOutcomes); nCohort = numel(cohorts);
nFeatType = numel(featType); nTextType = numel(textType);
nText = 0;
for t = 1:nTextType
    nText = nText + numel(textName{t});
end

featureNames = {'Best PET texture - GTVp'; 'Best PET texture - GTVtot'; 'Best CT texture - GTVp'; 'Best CT texture - GTVtot'}; nFeature = numel(featureNames);
rsMat_Binary = zeros(nFeature,nOutcomes); pMat_Binary = zeros(nFeature,nOutcomes); pCell_Binary = cell(nFeature,nOutcomes); stringCell_Binary = cell(nFeature,nOutcomes); nameTextCell_Binary = cell(nFeature,nOutcomes);
scans = {'PT','CT'}; nScans = numel(scans); nParams = [numel(scale_mat)*numel(algo_cell)*numel(Ng_mat),numel(scale_mat)*numel(algo_cell)*numel(Ng_mat)]; paramSizes = {[numel(scale_mat),numel(algo_cell),numel(Ng_mat)],[numel(scale_mat),numel(algo_cell),numel(Ng_mat)]};
nPat = 0;
for c = 1:nCohort
    nPat = nPat + nPatient.(cohorts{c});
end
i = 0;
for scan = 1:nScans
    nParam = nParams(scan);
    paramSize = paramSizes{scan};
    textParamName = cell(1,nText*nParam); count = 0;
    for t = 1:numel(textType)
        for tStr = 1:numel(textName{t})
            for p = 1:nParam
                count = count + 1;
                [aa,bb,cc] = ind2sub(paramSize,p);
                textParamName{count} = [textType{t},'/',textName{t}{tStr},' -- ','Scale=',num2str(scale_mat(aa)),',Quant.algo=',algo_cell{bb},',Ng=',num2str(Ng_mat(cc))];
            end
        end
    end
    for type = 1:nFeatType
        i = i + 1; textMat = zeros(nPat,nText*nParam); countPat = 0;
        outcomeMat = zeros(nPat,nOutcomes);
        cd(pathFeatures)
        for c = 1:nCohort
            cohort = cohorts{c};
            text = load(['text_',cohort,'_',scans{scan},'_',featType{type}]); text = struct2cell(text); text = text{1};
            patientVect = (countPat+1):(countPat+nPatient.(cohort)); countPat = countPat + nPatient.(cohort);
            count = 0;
            for t = 1:numel(textType)
                for tStr = 1:numel(textName{t})
                    for p = 1:nParam
                        count = count + 1;
                        [aa,bb,cc] = ind2sub(paramSize,p);
                        textMat(patientVect,count) = text{aa,bb,cc}.(textType{t}).(textName{t}{tStr}).Data;
                    end
                end
            end
            for j = 1:nOutcomes
                outcomeMat(patientVect,j) = outcomes.(cohort).(nameOutcomes{j});
            end
        end
        cd(pathUnivariate)
        for j = 1:nOutcomes
            outcome = outcomeMat(:,j);
            [rsTemp,pTemp] = corr(textMat,outcome,'type','Spearman','rows','pairwise'); temp = abs(rsTemp);
            [significance] = benjamini_hochberg(pTemp,0.10); proportion = sum(significance)/numel(significance);
            if proportion > 0
                starB = '*';
            else
                starB = '';
            end
            fprintf(fid,['Proportion of significant textures: ',scans{scan},', ',featType{type},', ',nameOutcomes{j},', Binary: ',num2str(proportion),'\n']);
            [~,indMax] = max(temp);
            rsMat_Binary(i,j) = rsTemp(indMax); pMat_Binary(i,j) = pTemp(indMax); nameTextCell_Binary{i,j} = textParamName{indMax};
            pCell_Binary{i,j} = pTemp;
            if pMat_Binary(i,j) < 0.01
                stringCell_Binary{i,j} = [starB,'rs = ',num2str(rsMat_Binary(i,j),'%.2f'),', p = ',num2str(pMat_Binary(i,j),'%.2i'),starB];
            else
                stringCell_Binary{i,j} = [starB,'rs = ',num2str(rsMat_Binary(i,j),'%.2f'),', p = ',num2str(pMat_Binary(i,j),'%.2f'),starB];
            end
        end
    end
end
cd(pathUnivariate), fclose(fid);
Locoregional = stringCell_Binary(:,1); Distant = stringCell_Binary(:,2); Death = stringCell_Binary(:,3);
results.tableCorr = table(Locoregional,Distant,Death,'RowNames',featureNames);
Locoregional = nameTextCell_Binary(:,1); Distant = nameTextCell_Binary(:,2); Death = nameTextCell_Binary(:,3);
results.tableTextName = table(Locoregional,Distant,Death,'RowNames',featureNames);
results.rsMat = rsMat_Binary; results.pMat = pMat_Binary; results.pCell = pCell_Binary;
save('texturesBest_UniV_Binary','results'), clear results

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%% t_meshFibersOBJ
%
%  Shows how to create OBJ fiber files from AFQ fibers.
%
% BW/LMP

%% Download a small set of fibers in a pdb file

remote = 'http://scarlet.stanford.edu/validation/MRI/VISTADATA';
remoteF = 'diffusion/sampleData/fibers/leftArcuateSmall.pdb';
remoteF = fullfile(remote,remoteF);
tmp = [tempname,'.pdb'];
[fgFile, status] = urlwrite(remoteF,tmp);

% Read the fiber groups
fg = fgRead(fgFile);

% Render the Tract FA Profile for the left uncinate
% A small number of triangles (25 is the default).
nTriangles = 4;
[lgt , fSurf, fvc] = AFQ_RenderFibers(fg,'subdivs',nTriangles);
% mrvNewGraphWin; surf(fSurf.X{1},fSurf.Y{1},fSurf.Z{1},fSurf.C{1})
% mrvNewGraphWin; plot3(FV.vertices(:,1),FV.vertices(:,2),FV.vertices(:,3),'.');

%% Accumulate fascicles into a structure that we can write using objWrite

% The obj structure has one list of vertices and one list of normals.
% It can then have multiple groups of fascicles (e.g., f1, f2, ...)
% But we end up accumulating the vertices from all the fascicles, and when
% we descrive the faces for, say, f2, they have to refer to the vertices
% for f2, not the vertices for f1.  So, we need to add an offset value to
% the faces.
FV.vertices = [];
FV.faces    = [];
N           = [];
% select = 1:10;  % Small for debugging.  Select only some faces.
cnt = 0;
nFascicles = size(fvc,2);
c    = ones(nFascicles,3);
redF = round(nFascicles)/2;
c(1:redF,:) = repmat([1 0 0],redF,1);
% c(1:redF,:) = repmat([1 0 0],redF,1);

for ff = 1:2:size(fvc,2)

    % We expertimented with color, and this worked in meshLab but not in
    % brainbrowser
    %
    % When we add color, we do it this way by appending RGB to the
    % vertex, and dealing with the first case separately
    if isempty(FV.vertices)
        % FV.vertices = [fvc(ff).vertices repmat(c(ff,:),size(fvc(ff).vertices,1),1)];
        FV.vertices = [fvc(ff).vertices repmat(c(ff,:),size(fvc(ff).vertices,1),1)];
    else
        % Vertices of the triangles defining the fascicle mesh
        FV.vertices = [FV.vertices; [fvc(ff).vertices repmat(c(ff,:),size(fvc(ff).vertices,1),1)]];
    end

    % Cumulate the vertices
    % FV.vertices = [FV.vertices; fvc(ff).vertices];

    % Cumulate the normals for each vertex
    [Nx,Ny,Nz] = surfnorm(fSurf.X{ff},fSurf.Y{ff},fSurf.Z{ff});
    tmp =[Nx(:),Ny(:),Nz(:)];
    N = [N ; tmp];

    % Add an offset to the faces, to make them consistent with cumulating
    % vertices.
    FV.faces    = [FV.faces; fvc(ff).faces + cnt];

    % Update where we are
    cnt = size(FV.vertices,1);

end

%% Format the OBJ data and write them out

OBJ = objFVN(FV,N);

fname = '/Users/wandell/Desktop/testArcuate.obj';
% fname = '/home/wandell/Desktop/testArcuate.obj';
objWrite(OBJ,fname);

%% Copy the data onto SDM
pLink = 'https://sni-sdm.stanford.edu/api/acquisitions/558da2ba3113bb9e05daaf0f/file/1.3.12.2.1107.5.2.32.35381.2015012510504883990801657.0.0.0_nifti.nii.gz?user=';
uName = 'wandell@stanford.edu';

%
sdmPut(pLink,uName,fname);

%%  These commands create the data and put the method as a PDF into the SDM
%
%   pdfFile = publish('t_meshFibersOBJ.m','pdf');
%   sdmPut(pLink,uName,pdfFile); disp('Done');
%

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_suite = test_parcellfun
% tests for cosmo_cartprod
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_parcellfun_single_proc()
    nproc=1;
    helper_test_parcellfun(nproc);

function test_cosmo_parallel_get_nproc_available()
    warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('reset');
    cosmo_warning('off');

    nproc=cosmo_parallel_get_nproc_available();
    assert(nproc>=1);

    % should not have shown any warnings
    w=cosmo_warning();
    assert(isempty(w.shown_warnings),sprintf('warning shown: %s',...
                                      w.shown_warnings{:}));

function test_parcellfun_multi_proc()
    nproc=cosmo_parallel_get_nproc_available();
    if nproc==1
        cosmo_notify_test_skipped('No parallel process available');
        return
    end

    helper_test_parcellfun(nproc);

function helper_test_parcellfun(nproc)
    warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    % try various functions
    funcs={@func_identity,...
            @func_reverse,...
            @numel,...
            @(x)numel(x)>0,...
            };

    % try with both uniform output and without
    arg_cell={{},...
                {'UniformOutput',false}};

    % various number of dimensions
    ndim={0,1,2,3,4}; % 0=empty

    combis=cosmo_cartprod({funcs,arg_cell,ndim});
    n=size(combis,1);

    for c_i=1:n
        combi=combis(c_i,:);
        func_arg=combi{1};
        other_arg=combi{2};
        rand_str_ndim=combi{3};

        rand_cellstr=generate_rand_cellstr(rand_str_ndim);
        func=@()cosmo_parcellfun(nproc,func_arg,rand_cellstr,...
                                other_arg{:});

        ref_func=@() cellfun(func_arg,rand_cellstr,other_arg{:});
        assert_equal_result_or_both_exception_thrown(func,ref_func);
    end

function assert_equal_result_or_both_exception_thrown(f, g)
    try
        f_result=f();
        f_exception=false;
    catch
        f_exception=lasterror();
    end

    try
        g_result=g();
        g_exception=false;
    catch
        g_exception=lasterror();
    end

    f_threw_exception=isstruct(f_exception);
    g_threw_exception=isstruct(g_exception);

    if f_threw_exception
        if ~g_threw_exception
            f_exception.message=sprintf('only f threw exception: %s',...
                                        f_exception.message);
            rethrow(f_exception);
        end
    else
        if g_threw_exception
            g_exception.message=sprintf('only g threw exception: %s',...
                                        g_exception.message);
            rethrow(g_exception);
        end

        assertEqual(f_result,g_result);
    end

function test_illegal_arguments
    warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    aet=@(varargin)assertExceptionThrown(@()...
                            cosmo_parcellfun(varargin{:}),'');

    % first argument is not scalar integer
    aet(0,@numel,{1,2});
    aet(1.5,@numel,{1,2});
    aet([2 3],@numel,{1,2});

    % second argument is not a function handle
    aet(2,'a',{1,2});
    aet(2,cell(0),{1,2});

    % third argument is not a cell
    aet(2,@numel,[1,2]);
    aet(2,@numel,struct);

    % no uniform output
    % (note: Octave accepts output when using
    %           @(x)[x x]
    %  This may be a bug)
    aet(1,@(x) repmat(x,1,x),{1;2});
    aet(2,@(x) repmat(x,1,x),{1;2});

function rand_cellstr=generate_rand_cellstr(ndim)
    switch ndim
        case 0
            sz=[0,0];

        case 1
            sz=[0,10];

        case 2
            sz=[1,1];

        otherwise
            sz=floor(rand(1,1+ndim)*4);
    end

    randstr=@(unused)char(rand(1,floor(rand()*10))*24+65);
    rand_cellstr=arrayfun(randstr,zeros(sz),'UniformOutput',false);

function y=func_identity(x)
    y=x;

function y=func_reverse(x)
    y=x;
    y(end:-1:1)=x(:);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [ a, b ] = p23_lim ( dim_num )

%*****************************************************************************80
%
%% P23_LIM returns the integration limits for problem 23.
%
%  Discussion:
%
%    Because the integration region is the interior of the unit simplex,
%    the integration limits simply specify the limits of a box containing
%    the integration region.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end


text:
%  translate the inputs from RISE into inputs for csminwel
% --------------------------------------------------------
%

author: jmaih


text:
function id = getlmiid(F)

id = F.LMIid;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [Grad,t]=StdRect(tStart,tEnd,gAmp,step)

%tStart start time
%tEnd end time
%step increment steps
%gAmp amplitude

t=linspace(tStart,tEnd,step);
Grad=gAmp*ones(size(t)); % rectangle

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function S = oblateSurfaceArea(elli, varargin)
%OBLATESURFACEAREA  Approximated surface area of an oblate ellipsoid
%
%   S = oblateSurfaceArea(R1,R2)
%
%   Example
%   oblateSurfaceArea
%
%   See also
%   geom3d, ellipsoidSurfaceArea, prolateSurfaceArea
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2015-07-03,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2015 INRA - Cepia Software Platform.

%% Parse input argument

if size(elli, 2) == 7
    R1 = elli(:, 4);
    R2 = elli(:, 5);

elseif size(elli, 2) == 1 && ~isempty(varargin)
    R1 = elli(:, 1);
    R2 = varargin{1};
end

assert(R1 < R2, 'First radius must be smaller than second radius');

% surface theorique d'un ellipsoide oblate
% cf http://fr.wikipedia.org/wiki/Ellipso%C3%AFde_de_r%C3%A9volution
e = sqrt(R2.^2 - R1.^2) ./ R2;
S = 2 * pi * R2.^2 + pi * R1.^2 * log((1 + e) ./ (1 - e)) ./ e;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function writedata (x, y, yMax, yMin, p, filename, commentStr)
% WRITEDATA creates txt file with parameters and computed values
% writedata (x, y, yMax, yMin, p, filename, commentStr)
% x - values on x-axis
% y - values on y-axis
% p - parameters used for computation
% filename - name of the txt file
% commentStr - comment string
fid = fopen([filename '.txt'],'wt');

if nargin>6
    fprintf(fid,'%s \n', commentStr);
end

fprintf(fid,'Computation: %s \n\n', datestr(now));

fprintf(fid,'Parameters:\n');

fn =fieldnames(p);
for ii = 1: length(fn)
    name = fn{ii};
    value = getfield(p,fn{ii});
    if ischar (value)
        fprintf (fid, [name ' = ' value '\n']);
    else
        fprintf (fid, [name ' = %g \n'], value);
    end
end

fprintf (fid,'\nX \t\t\t Y \n');
fprintf (fid,'------------------\n');

for ii=1:length(y)
    fprintf (fid,'%-8.3f \t %-8.3f \n', x(ii), y(ii));
end

fclose(fid);

%envelopes

if ~isempty(yMax)
    fid = fopen([filename '_envelope_Max.txt'],'wt');
    if nargin>6
        fprintf(fid,'%s \n', commentStr);
    end
    fprintf(fid,'Computation: %s \n\n', datestr(now));
    fprintf(fid,'Maximum of %g simulations of Poisson process\n\n', p.Nsimul);
    fprintf (fid,'x \t\t\t max y \n');
    fprintf (fid,'------------------\n');
    for ii=1:length(y)
        fprintf (fid,'%-8.3f \t %-8.3f \n', x(ii), yMax(ii));
    end
    fclose(fid);
end

if ~isempty(yMin)
    fid = fopen([filename '_envelope_Min.txt'],'wt');
    if nargin>6
        fprintf(fid,'%s \n', commentStr);
    end
    fprintf(fid,'Computation: %s \n\n', datestr(now));
    fprintf(fid,'Minimum of %g simulations of Poisson process\n\n', p.Nsimul);
    fprintf (fid,'xK \t\t\t min K \n');
    fprintf (fid,'------------------\n');
    for ii=1:length(y)
        fprintf (fid,'%-8.3f \t %-8.3f \n', x(ii), yMin(ii));
    end
    fclose(fid);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function ihs_test03 ( )

%*****************************************************************************80
%
%% IHS_TEST03 tests the improved distributed hypercube sampling algorithm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2014
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  duplication = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'IHS_TEST03' );
  fprintf ( 1, '  IHS implements the IHS Algorithm\n' );
  fprintf ( 1, '  (Improved Distributed Hypercube Sampling)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Demonstrate the code for a fixed dimension\n' );
  fprintf ( 1, '  and duplication value, and increasing number of points.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension =        %d\n', dim_num );
  fprintf ( 1, '  Duplication factor =       %d\n', duplication );

  for i = 1 : 5

    point_num = 10 * 2^(i-1);

    opt = point_num / point_num^( 1.0 / dim_num );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of points =         %d\n', point_num );
    fprintf ( 1, '  Desired minimum distance = %f\n', opt );
%
%  Get the points.
%
    x = ihs ( dim_num, point_num, duplication );
%
%  Compute the covariance.
%
    [ average, sd, covc ] = covariance ( dim_num, point_num, x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Average minimum distance %f\n', average );
    fprintf ( 1, '  Standard deviation:      %f\n', sd );
    fprintf ( 1, '  Covariance:              %f\n', covc );

    fprintf ( 1, '\n' );

    for j = 1 : point_num

      if ( j <= 10 | point_num - 10 <= j )
        fprintf ( 1, '%4d    ', j );
        for i = 1 : dim_num
          fprintf ( 1, '%4d  ', x(i,j) );
        end
        fprintf ( 1, '\n');
      elseif ( j == 11 )
        fprintf ( 1, '....  ........\n' );
      end
    end

  end

  return
end


text:
function [ code, itree, more ] = tree_parent_next ( nnode, more )

%*****************************************************************************80
%
%% TREE_PARENT_NEXT generates, one at a time, all labeled trees.
%
%  Discussion:
%
%    The routine also returns the corresponding Pruefer codes.
%
%    There are N^(N-2) labeled trees on N nodes (Cayley's formula).
%
%    The number of trees in which node I has degree D(I) is the
%    multinomial coefficient: ( N-2; D(1)-1, D(2)-1, ..., D(N)-1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Parameters:
%
%    Input, integer NNODE, the number of nodes to be used in
%    the trees.
%
%    Output, integer CODE(NNODE).  The first NNODE-2 entries
%    of CODE contain the Pruefer code for the given labeled tree.
%
%    Output, integer ITREE(NNODE).  The first NNODE-1 entries
%    of ITREE describe the edges that go between the nodes.  Each pair
%    (I, ITREE(I)) represents an edge.  Thus if ITREE(5) = 3,
%    there is an edge from node 3 to node 5.
%
%    Input/output, logical MORE.  On the first call only, the
%    user is required to set MORE = .FALSE.  Then call the routine, and
%    it will return information about the first tree
%    as well as setting MORE to the value .TRUE.
%    Keep calling to get another tree until MORE is .FALSE.
%    on return, at which point there are no more trees.
%
  persistent iarray;

  if ( ~ more )
    iarray = [];
  end

  [ iarray, more ] = vec_next ( nnode-2, nnode, iarray, more );

  code(1:nnode-2) = iarray(1:nnode-2) + 1;

  itree = pruefer_to_tree_2 ( nnode, code );

  return
end


text:
function ppg_filt(up)
%PPG_FILT extracts respiratory signals using various filtering techniques
% from the PPG signal as specified in PC's literature review.
%	            ppg_filt(option, up)
%
%	Inputs:
%		option          the option which has led to this function being used
%       up              universal parameters structure
%
%	Outputs:
%       ...
%

fprintf('\n--- Extracting Respiratory Signals from PPG using Filtering Techniques ');
log_int_respSig = 0;             % Has value 1 unless this is a final respiratory signal

for subj = up.paramSet.subj_list

    %% Cycle through each ppg signal
    for sig_no = 1 : length(up.paramSet.ppg_sigs)

        %% Cycle through each method
        for filt_no = 1 : length(up.al.options.ppg_filt)

            %% Skip if this processing has been done previously
            eval(['save_name = ''' up.paramSet.ppg_sigs{sig_no}, up.paths.filenames.filt '_' up.al.options.ppg_filt{filt_no} ''';']);
            iden_resp_sig_file_ending
            savepath = [up.paths.data_save_folder, num2str(subj), ending];
            exist_log = check_exists(savepath, save_name);
            if exist_log
                continue
            end

            %% Load relevant data
            if ~exist('data', 'var')
                load([up.paths.data_load_folder, up.paths.data_load_filename]);
            end
            % Extract PPG data
            eval(['rel_data = data(subj).' up.paramSet.ppg_sigs{sig_no} ';']);

            %% Filter the raw signal using this method
            respWave = feval(up.al.options.ppg_filt{filt_no}, rel_data, up);

            %% Band-pass filter
            filtered_data = bpf_signal_to_remove_non_resp_freqs(respWave, respWave.fs, up);
            eval([save_name ' = filtered_data;']);

            %% Save processed data
            save_or_append_data

        end

    end

end

end

author: peterhcharlton

repo: RRest


text:
function [vertices, faces] = curveToMesh(curve, varargin)
%CURVETOMESH  Create a mesh surrounding a 3D curve.
%
%   [V, F] = curveToMesh(CURVE)
%   Computes the vertices and the faces of the mesh surrounding the
%   specified 3D curve.
%
%   [V, F] = curveToMesh(CURVE, THICKNESS)
%   Specifies the thickness of the mesh (distance between mesh vertices and
%   curve vertices). Default is 0.5.
%
%   [V, F] = curveToMesh(CURVE, THICKNESS, NCORNERS)
%   Also specifies the number of mesh vertices around each curve vertex.
%   Default is 8.
%
%
%   Example
%     % Creates a tubular mesh around a trefoil knot curve
%     t = linspace(0, 2*pi, 200)';
%     x = sin(t) + 2 * sin(2 * t);
%     y = cos(t) - 2 * cos(2 * t);
%     z = -sin(3 * t);
%     curve = [x, y, z];
%     [v2, f2] = curveToMesh(curve, .5, 16);
%     figure;
%     drawMesh(v2, f2);
%     axis equal; view(3);
%     axis([-4 4 -4 4 -2 2]);
%
%   See also
%     meshes3d, torusMesh, surfToMesh

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2015-01-07,    using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015 INRA - Cepia Software Platform.

radius = .1;
if nargin > 1
    radius = varargin{1};
end

nCorners = 8;
if nargin > 2
    nCorners = varargin{2};
end

nNodes = size(curve, 1);
nVerts = nNodes * nCorners;

vertices = zeros(nVerts, 3);

% create reference corners, that will be rotated and translated
t = linspace(0, 2*pi, nCorners + 1)';
t(end) = [];
baseCorners = radius * [cos(t) sin(t) zeros(size(t))];

for iNode = 1:nNodes
    % coordinate of current node
    node = curve(iNode, :);

    % compute local tangent vector
    iNext = mod(iNode, nNodes) + 1;
    tangentVector = normalizeVector3d(curve(iNext, :) - node);

    % convert to spherical coordinates
    [theta, phi, rho] = cart2sph2(tangentVector); %#ok<ASGLU>

    % apply transformation to place corners around current node
    rotY = createRotationOy(theta);
    rotZ = createRotationOz(phi);
    trans = createTranslation3d(node);
    transformMatrix = trans * rotZ * rotY;
    corners = transformPoint3d(baseCorners, transformMatrix);

    % concatenate with other corners
    vertices( (1:nCorners) + (iNode - 1) * nCorners, :) = corners;
end

% indices of vertices
inds = (1:nVerts)';
add1 = repmat([ones(nCorners-1, 1) ; 1-nCorners], nNodes, 1);

% generate faces
faces = [inds ...
    mod(inds + add1 - 1, nVerts) + 1 ...
    mod(inds + nCorners + add1 - 1, nVerts) + 1 ...
    mod(inds + nCorners - 1, nVerts) + 1];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_KUKA_KR16_2(robot, T)
%   Solves the inverse kinematic problem for the KUKA KR16 2 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_KUKA_KR16_2 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('kuka', 'KUKA_KR16_2');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function q = inversekinematic_kuka_kr16_2(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=abs(d(6));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));
% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist(robot, q(:,i), T, 1,'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));
A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

A2 =(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi - beta;
q3(2) = pi - phi + beta;


text:
function [maxVal] = maxIdx(a, idx)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

maxVal = max(a(idx, :), [], 1);

author: adikhosla

repo: feature-extraction


text:
function lebesgue_test04 ( )

%*****************************************************************************80
%
%% LEBESGUE_TEST04 looks at Chebyshev4 points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEBESGUE_TEST04:\n' );
  fprintf ( 1, '  Analyze Chebyshev4 points.\n' );

  n_max = 11;
  l = zeros ( n_max, 1 );
  xfun = linspace ( -1.0, +1.0, 501 );
  xfun = xfun';

  for n = 1 : n_max
    x = chebyshev4 ( n );
    l(n) = lebesgue_constant ( n, x, xfun );
  end

  r8vec_print ( n_max, l, '  Chebyshev4 Lebesgue constants for N = 1 to 11:' )
%
%  Examine one case more closely.
%
  n = 11;
  x = chebyshev4 ( n );
  r8vec_print ( n, x, '  Chebyshev4 points for N = 11' );

  label = 'Chebyshev4 points for N = 11';
  filename = 'chebyshev4.png';
  lebesgue_plot ( n, x, xfun, label, filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file saved as "s"\n', filename );

  return
end

repo: jburkardt-m


text:
% IMAGESCTC - DEPRECATED. never completed or documented.
% IMAGESCTC - imagesc in true color. Can help plot different
%               colormap on the same window.
%
% Usage: same as imagesctc
%
% Example:
%         figure;
%         colormap(jet);
%         subplot(1,2,1); imagesctc(rand(10,10));
%         colormap(gray);
%         subplot(1,2,2); imagesctc(rand(10,10));
%
% Author: Arnaud Delorme, CNL / Salk Institute, 31 July 2002
%
% See also: IMAGE, IMAGESC

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

if ~isempty(varargin)
	imagesc(a, varargin{:});
else
	imagesc(a);
end
c = caxis;

cm = colormap;
intervals = linspace(c(1), c(2), size(cm,1));
[tmp dest] = histc(a(:), intervals);

for wi = 1:size(a,2)
	for hi = 1:size(a,1)
		aa(hi,wi,1:3) = cm(dest(hi+(wi-1)*size(a,1)),1:3);
	end
end
if ~isempty(varargin)
	imagesc(aa, varargin{:});
else
	imagesc(aa);
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [cm] = mm2cm(mm)
% Convert length from millimeters to centimeters.
% Chad A. Greene 2012
cm = mm/10;


text:
%  INTERNAL FUNCTION: Computes kron(T1,T2)*vec(X)
%
%


text:
function f = getfigurepos(fig,units)

% function f = getfigurepos(fig,units)
%
% <fig> (optional) is the figure handle.  default: gcf.
% <units> (optional) is 'normalized' (default) | 'points'
%
% return the position of <fig> in units specified by <units>.
%
% example:
% figure; getfigurepos(gcf,'points')

% input
if ~exist('fig','var') || isempty(fig)
  fig = gcf;
end
if ~exist('units','var') || isempty(units)
  units = 'normalized';
end

% do it
  prev = get(fig,'Units');
  set(fig,'Units',units);
f = get(fig,'Position');
  set(fig,'Units',prev);

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function sigma_values_test ( )

%*****************************************************************************80
%
%% SIGMA_VALUES_TEST demonstrates the use of SIGMA_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIGMA_VALUES_TEST:\n' );
  fprintf ( 1, '  SIGMA_VALUES returns values of \n' );
  fprintf ( 1, '  the SIGMA function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N         SIGMA(N)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, fn ] = sigma_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %10d\n', n, fn );

  end

  return
end

author: johannesgerer


text:
classdef HBox < uix.HBox
    %uiextras.HBox  Arrange elements in a single horizontal row
    %
    %   obj = uiextras.HBox() creates a new horizontal box layout with
    %   all parameters set to defaults. The output is a new layout object
    %   that can be used as the parent for other user-interface components.
    %
    %   obj = uiextras.HBox(param,value,...) also sets one or more
    %   parameter values.
    %
    %   See the documentation for more detail and the list of properties.
    %
    %   Examples:
    %   >> f = figure();
    %   >> b = uiextras.HBox( 'Parent', f );
    %   >> uicontrol( 'Parent', b, 'Background', 'r' )
    %   >> uicontrol( 'Parent', b, 'Background', 'b' )
    %   >> uicontrol( 'Parent', b, 'Background', 'g' )
    %   >> set( b, 'Sizes', [-1 100 -2], 'Spacing', 5 );
    %
    %   >> f = figure();
    %   >> b1 = uiextras.VBox( 'Parent', f );
    %   >> b2 = uiextras.HBox( 'Parent', b1, 'Padding', 5, 'Spacing', 5 );
    %   >> uicontrol( 'Style', 'frame', 'Parent', b1, 'Background', 'r' )
    %   >> uicontrol( 'Parent', b2, 'String', 'Button1' )
    %   >> uicontrol( 'Parent', b2, 'String', 'Button2' )
    %   >> set( b1, 'Sizes', [30 -1] );
    %
    %   See also: uiextras.VBox
    %             uiextras.HBoxFlex
    %             uiextras.Grid

    %  Copyright 2009-2014 The MathWorks, Inc.
    %  $Revision: 1077 $ $Date: 2015-03-19 16:44:14 +0000 (Thu, 19 Mar 2015) $

    properties( Hidden, Access = public, Dependent )
        Enable % deprecated
        Sizes
        MinimumSizes
    end

    methods

        function obj = HBox( varargin )

            % Call uix constructor
            obj@uix.HBox( varargin{:} )

            % Auto-parent
            if ~ismember( 'Parent', varargin(1:2:end) )
                obj.Parent = gcf();
            end

        end % constructor

    end % structors

    methods

        function value = get.Enable( ~ )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

            % Return
            value = 'on';

        end % get.Enable

        function set.Enable( ~, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, {'on','off'} ) ), ...
                'uiextras:InvalidPropertyValue', ...
                'Property ''Enable'' must be ''on'' or ''off''.' )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

        end % set.Enable

        function value = get.Sizes( obj )

            % Get
            value = transpose( obj.Widths );

        end % get.Sizes

        function set.Sizes( obj, value )

            % Set
            obj.Widths = value;

        end % set.Sizes

        function value = get.MinimumSizes( obj )

            % Get
            value = transpose( obj.MinimumWidths );

        end % get.MinimumSizes

        function set.MinimumSizes( obj, value )

            % Get
            obj.MinimumWidths = value;

        end % set.MinimumSizes

    end % accessors

end % classdef

author: goGPS-Project


text:
function fx = p37_fun ( n, x )

%*****************************************************************************80
%
%% P37_FUN evaluates the integrand for problem 37.
%
%  Discussion:
%
%    The problem has a parameter ALPHA that can be set by calling
%    P37_PARAM_SET.  It had a default value of 5.0.
%
%    The integrand has a peak of height 4^ALPHA at X = PI/4.
%
%    Suggested values of ALPHA include 0 through 20.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    4^(-ALPHA) / ( ( x - PI/4 )^2 + 16^(-ALPHA) )
%
%  Exact Integral:
%
%    atan ( ( 4 - PI ) * 4^(ALPHA-1) ) + atan ( PI * 4^(ALPHA-1) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 83.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  alpha = p37_param_get ( );

  fx = 4.0^( -alpha ) ./ ( ( x - 0.25 * pi ).^2 + 16.0^(-alpha) );

  return
end


text:
function [gcjLat, gcjLng]  = bd2gcj(bdLat, bdLng)
    gcjLat = bdLat;
    gcjLng = bdLng;

    x_pi = pi * 3000.0 / 180.0;
    inChina = ~outOfChina(bdLat, bdLng);
    if ~any(inChina),return;end

    x = bdLng(inChina) - 0.0065;
    y = bdLat(inChina) - 0.006;

    z = hypot(x, y) - 0.00002 * sin(y * x_pi);
    theta = atan2(y, x) - 0.000003 * cos(x * x_pi);
    gcjLng(inChina) = z.*cos(theta);
    gcjLat(inChina) = z.*sin(theta);

end

author: googollee

repo: eviltransform


text:
function [S_sorted,s]=sort_ordinal(S)
% SORT_ORDINAL reorders nodes to emphasize persistent structure in an ordered multilayer partition
%
% Version: 2.2.0
% Date: Thu 11 Jul 2019 12:25:42 CEST
%
% Nodes are reordered using the optimal leave order for the
% average linkage hierarchical clustering tree based on Hamming distance
% between community assignments
%
% Call as:
%
%     [S_sorted, s] = sort_ordinal(S)
%
% Input:
%
%     S: multilayer partition (matrix of size NxT, where N is the number
%        of nodes and T is the number of layers)
%
% Output:
%
%     S_sorted: reordered multilayer partition
%
%     s: mapping of nodes to reordered nodes
%
% Note that S_sorted=S(s, :)

d1=pdist(S,'hamming');
Z1=linkage(d1,'average');
s=optimalleaforder(Z1,d1);
S_sorted=S(s,:);

author: GenLouvain

repo: GenLouvain


text:
% M1QN3  Solve a UNO using M1QN3
%
% THIS IS A LOW LEVEL FUNCTION - USE opti_m1qn3() INSTEAD!
%
%   [x,fval,exitflag,iter,feval] = m1qn3(fun,grad,x0,opts)
%
%   Input arguments:
%       fun - nonlinear function handle
%       grad - nonlinear gradient function handle (required)
%       x0 - initial solution guess
%       opts - solver options (see below)
%
%   Return arguments:
%       x - solution vector
%       fval - objective value at the solution
%       exitflag - exit status (see below)
%       iter - number of iterations taken by the solver
%       feval - number of function evaluations taken by the solver
%
%   Option Fields (all optional):
%       display - solver display level [0,1,2]
%       tolafun - gradient convergence tolerance
%       maxiter - maximum solver iterations {1000}
%       maxfeval - Maximum number of function evaluations {1500}
%       maxtime - Maximum solver execution time {1000s}
%       nupdate - Number of L-BFGS updates for Hessian (5-10) {5}
%       iterfun - Iteration callback function handle
%
%   Return Status:
%       0 - user exited
%       1 - successful termination
%       2 - input argument error
%       3 - line-search is blocked
%       4 - maximum iterations reached
%       5 - maximum function evaluations reached
%       6 - stop on dxmin during the line-search
%       7 - numerical error
%       8 - maximum time reached
%
%
%   Copyright (C) 2012 Jonathan Currie (I2C2)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function plotFunc()
  colors = ['r'; 'b';'m';'g'];
  lineWidth=2;
  markerSize=20;
  sigmas = [2.5, 5, 10];

  hold on;
  legendStrs = {};
  x=-20:0.2:20;
  for ii=1:length(sigmas)
    sigma = sigmas(ii);
    y = exp(-0.5*(x/sigma).^2);
    plot(x, y, colors(ii), 'MarkerSize', markerSize, 'linewidth', lineWidth);
    legendStrs{end+1} = ['\sigma=' num2str(sigma)];
  end

  % legend
  fontsize = 12;
  legendLocation = 'Best'; % 'NorthWest'; %
  legend(legendStrs,'FontSize', fontsize, 'Location', legendLocation);
end

author: lmthang


text:
function [resampSigM,TTHPv,SHPv,timeOutV] = halfPeak(ROIDataM,timeV,smoothFlag,resampFlag)
% =======================================================================================
% AI  09/28/16
% AI  10/20/16  Added smoothing/resampling options.
% ======================================================================================
% INPUTS
% ROIDataM   : (no. voxels x time pts) Matrix of temporal signals of
%              raster-ordered voxels along rows
% timeV      : Vector of time pts of DCE image acquisition.
% smoothFlag : Flag (1-Smooth    0-No smoothing )
% resampFlag : Flag (1-Resample  0-No resampling )
% =====================================================================================

%% Smooth/resample if selected
nVox = size(ROIDataM,1);
[resampSigM,timeOutV] = smoothResamp(ROIDataM,timeV,smoothFlag,resampFlag);
normalizedBaseline = 1;
relSigEnhM = resampSigM - normalizedBaseline;

%% Compute TTHP, SHP
%Identify first time point where changeSig > half-peak
gts50 = bsxfun(@ge,relSigEnhM,(max(resampSigM,[],2) - normalizedBaseline)/2);
[~,SHPcolIdx] = max(gts50,[],2);
SHProwIdx = (1:nVox).';
SHPidx = sub2ind([nVox,numel(timeOutV)],SHProwIdx,SHPcolIdx);
%Get signal at half-peak
SHPv = resampSigM(SHPidx).';
%Get time to half-peak
TTHPv = timeOutV(1,SHPcolIdx);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR

path: github-repos/MATLAB/cerr-CERR/CERR-d320754abad9dcb78508ab69f33ae9f644202114/CERR_core/PlanAnalysis/DCE-MR analysis/TTHP analysis for DCE-MR/halfPeak.m


text:
function sphere_fibonacci_grid_points_test ( )

%*****************************************************************************80
%
%% SPHERE_FIBONACCI_GRID_POINTS_TEST tests SPHERE_FIBONACCI_GRID_POINTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_FIBONACCI_GRID_POINTS_TEST\n' );
  fprintf ( 1, '  SPHERE_FIBONACCI_GRID_POINTS computes points on a sphere\n' );
  fprintf ( 1, '  that lie on a Fibonacci spiral.\n' );

  ng = 1000;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points NG = %d\n', ng );

  xg = sphere_fibonacci_grid_points ( ng );

  r8mat_print_some ( ng, 3, xg, 1, 1, 20, 3, '  Part of the grid array:' );
%
%  Display the nodes on a sphere.
%
  filename = sprintf ( 'sphere_fibonacci_grid_n%d.png', ng );

  sphere_fibonacci_grid_display ( ng, xg, filename );
%
%  Write the nodes to a file.
%
  filename = sprintf ( 'sphere_fibonacci_grid_n%d.xyz', ng );

  r8mat_write ( filename, ng, 3, xg )

  return
end


text:
function varargout = gsexpand(varargin)
%GSEXPAND Generalized scalar expansion.
%
%   [XE, YE, ZE, ...] = GSEXPAND(X, Y, Z, ...) performs a generalized scalar
%   expansion on the input arguments.  The output arguments will all have
%   the same size.
%
%   All arguments are expanded (replicated) along singleton dimensions to
%   match the size of the other arguments.
%
%   See also RESIZE, SEXPAND.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-06-06 16:13:35 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % check number of input and output arguments
   %

   nargsin  = nargin;
   nargsout = nargout;

   % check the number of input arguments
   if nargsin < 2
      error('Not enough input arguments.');
   end

   % check the number of output arguments
   if nargsout == 0
      % MATLAB convention: when a function is called with no output
      % arguments, return one output argument
      nargsout = 1;
   else
      if nargsout > nargsin
         error('Too many output arguments.');
      end
   end

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % find the size of the output arguments
   %

   % initialize common size vector and number of dimensions to that of a
   % scalar
   csize  = [1, 1];
   cdims  = 2;

   % iterate over the input arguments
   for i = 1:nargsin

      % size and number of dimensions for i'th input argument
      isize = size(varargin{i});
      idims = length(isize);

      if idims <= cdims
         % i'th argument has no more dimensions than the output will have;
         % only check the `idims' lowest dimensions
         dims = idims;
      else
         % i'th argument has more dimensions than the current common size
         % vector; only compare the `cdims' lowest dimensions
         dims = cdims;

         % update `cdims' for the next round
         cdims = idims;
      end

      % find non-singleton dimensions
      insdims = isize(1:dims) ~= 1;
      cnsdims = csize(1:dims) ~= 1;

      % check size compatibility
      if any( insdims & cnsdims & (isize(1:dims) ~= csize(1:dims)) )
         error('Lengths must match along non-singleton dimensions.');
      end

      csize(insdims) = isize(insdims);

   end

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % initialize output argument list and call other program
   %

   varargout = cell(1, nargsout);
   [varargout{:}] = resize(csize, varargin{:});

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function mtrFigureMTCorrelations(machineDir, saveDir)

subjVec = {'ss040804','mho040625','bg040719','md040714'};
threshVec = [1000];
paramNames = {'kLength','kSmooth','kMidSD'};
midP = [-2 18 0.175];
bSaveImages = 1;

for pp = 1:3
    figure;
    for ss = 1:length(subjVec)
        subjDir = [machineDir subjVec{ss}];

        fgDir = [subjDir '/conTrack/resamp_LMT'];
        disp(['cd ' fgDir]);
        cd(fgDir);
        % Correlation between CTs
        mtrPlotCorr('cc',threshVec,paramNames,midP,pp,'b'); hold on;
        % Correlation between CTs and STTs
        %mtrPlotCorr('cc_STT',threshVec,paramNames,midP,1,'g'); hold on;
        % Save out the figures to the image dir

        fgDir = [subjDir '/conTrack/resamp_RMT'];
        disp(['cd ' fgDir]);
        cd(fgDir);
        % Correlation between CTs
        mtrPlotCorr('cc',threshVec,paramNames,midP,pp,'b'); hold on;

        oldaxis = axis;
        if pp == 1
            axis([-4 1 0 1]);
        else
            axis([oldaxis(1) oldaxis(2) 0 1]);
        end
        % Correlation between CTs and STTs
        %mtrPlotCorr('cc_STT',threshVec,paramNames,midP,1,'g');
    end
    % Save out the figures to the image dir
    set(gcf,'Position',[395   447   289   240]);
    if bSaveImages
        figFilename = fullfile(saveDir,['param_' paramNames{pp} '_corrMT.png']);
        set(gcf,'PaperPositionMode','auto');
        print('-dpng', figFilename);
    end
end

author: vistalab

repo: vistasoft


text:
function es = ESupdate(mu,lambda,iter)

% Create vector of ES weights
es.mu = mu;
es.lambda = lambda;
es.iter = iter;
tot = es.mu + es.lambda;
w = ceil((1./sqrt(1:tot))/sum(1./sqrt(1:tot))*es.lambda);
nonzero = sum(w > 0);
while (sum(w) - es.lambda) > nonzero
    w = max(0, w - 1);
    nonzero = sum(w > 0);
end
delta = sum(w) - es.lambda;
lastnonzero = find(w > 0, 1, 'last');
w(max(1,lastnonzero-delta+1):lastnonzero) = w(max(1,lastnonzero-delta+1):lastnonzero) - 1;

% Create selection mask
cw = cumsum(w) - w + 1;
idx(cw) = 1;
es.selectmask = cumsum(idx(1:end-1));

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [voiContourHandles] = matRad_plotVoiContourSlice(axesHandle,cst,ct,ctIndex,selection,plane,slice,cMap,varargin)
% matRad function that plots the contours of the segmentations given in cst
%
% call
% voiContourHandles = matRad_plotVoiContourSlice(axesHandle,cst,ct,ctIndex,selection,plane,slice)
% voiContourHandles = matRad_plotVoiContourSlice(axesHandle,cst,ct,ctIndex,selection,plane,slice,cMap)
% voiContourHandles = matRad_plotVoiContourSlice(axesHandle,cst,ct,ctIndex,selection,plane,slice,cMap,...)
%
% input
%   axesHandle          handle to axes the slice should be displayed in
%   cst                 matRad cst cell array
%   ct                  matRad ct structure
%   ctIndex             index of the ct cube
%   selection           logicals defining the current selection of contours
%                       that should be plotted. Can be set to [] to plot
%                       all non-ignored contours.
%   plane               plane view (coronal=1,sagittal=2,axial=3)
%   slice               slice in the selected plane of the 3D cube
%   cMap                optional argument defining the colormap, default is
%                       colorcube
%   varargin            Additional Matlab Line-Property/value pairs
%
% output
%   voiContourHandles:  handles of the plotted contours
%   visibleOnSlice:     logicals defining if the contour is actually
%                       visible on the current slice
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[env, ~] = matRad_getEnvironment();

% overwrite colormap
if exist('cMap', 'var') && ~isempty(cMap)
    cMapScale = size(cMap,1)-1;
    %determine colors
    colors = cMap(round(linspace(1,cMapScale,size(cst,1))),:);
else
    for i = 1:size(cst,1)
        if isfield(cst{i,5},'visibleColor')
            colors(i,:) = cst{i,5}.visibleColor;
        else
            colors(i,:) = [0 0 0];
        end
    end
end

if isempty(selection) || numel(selection) ~= size(cst,1)
    selection = true(size(cst,1),1);
end

voiContourHandles = cell(0);
switch env
    case 'MATLAB'
        voiContourHandles = gobjects(0);
    case 'OCTAVE'
        voiContourHandles = [];
end

for s = 1:size(cst,1)

    if ~strcmp(cst{s,3},'IGNORED') && selection(s)
        %Check for precalculated contours
        C =[];
        if size(cst,2) >= 7 && ~isempty(cst{s,7})
            C = cst{s,7}{slice,plane};
        else
            %If we do not have precomputed contours available, then compute them
            mask = zeros(size(ct{ctIndex}));
            mask(cst{s,4}{ctIndex}) = 1;

            if plane == 1 && any(any(mask(slice,:,:) > 0))
                C = contourc(squeeze(mask(slice,:,:)),.5*[1 1]);
            elseif plane == 2 && any(any(mask(:,slice,:) > 0))
                C = contourc(squeeze(mask(:,slice,:)),.5*[1 1]);
            elseif plane == 3 && any(any(mask(:,:,slice) > 0))
                C = contourc(squeeze(mask(:,:,slice)),.5*[1 1]);
            end
        end

        % plot precalculated contourc data

         switch env
            case 'MATLAB'
                tmpLineHandle = gobjects(0);
            case 'OCTAVE'
                tmpLineHandle = [];
         end

         if any(C(:))
            lower = 1; % lower marks the beginning of a section
            while lower-1 ~= size(C,2)
                hold on
                steps = C(2,lower); % number of elements of current line section
                tmpLineHandle(end+1) = line(C(1,lower+1:lower+steps),...
                    C(2,lower+1:lower+steps),...
                    'Color',colors(s,:),'Parent',axesHandle,varargin{:});

                lower = lower+steps+1;
            end
            voiContourHandles{end+1} = tmpLineHandle;
         else
            % create empty line object
            voiContourHandles{end+1} = {};
         end

    end
end

end

repo: matRad


text:
function [x0] = ma_getX0ForEvent(event)
%ma_getX0ForEvent Summary of this function goes here
%   Detailed explanation goes here

    activeVars = event.vars(1,:) == 1;
    switch event.type
        case 'Coast'
            x0 = event.coastToValue(activeVars);
        case 'NBodyCoast'
            x0 = event.coastToValue(activeVars);
        case 'DV_Maneuver'
            dvOrig = event.maneuverValue(activeVars);
            x0 = dvOrig;
        case 'Set_State'
            if(strcmpi(event.subType,'setState'))
                epoch = event.epoch;
                rVect = event.rVect;
                vVect = event.vVect;
                gmu = event.centralBody.gm;

                [sma, ecc, inc, raan, arg, tru] = getKeplerFromState(rVect,vVect,gmu, true);

                x0 = [epoch, sma, ecc, inc, raan, arg, tru];
                x0 = x0(activeVars);
            elseif(strcmpi(event.subType,'estLaunch'))
                x0 = event.launch.launchValue(activeVars);
            end
        otherwise
            error('Unknown event type, could not determine x0 for optimization: %s', event.type);
    end
end