text:
%% patchConnectivity
% Below is a demonstration of the features of the |patchConnectivity| function

%% Syntax
% |[C]=patchConnectivity(F,V,conType);|

%% Description
% This functions creates connectivity matrices for the input patch data
% defined by the faces F and the vertices V. The output is a structure
% containing the connectivity matrices:
%
% C.vertex.vertex
% C.vertex.face
% C.vertex.edge
%
% C.edge.face
% C.edge.vertex
% C.edge.edge
%
% C.face.vertex
% C.face.face
% C.face.edge
%
% If the 3rd optional input conType is not provided its default value is
% 'all' and all connectivity matrices are output in the structure. If not
% all types are desired the user may request only particular types by
% setting conType. The following conTypes can be specified:
% conType (v stands for vertex, f for face, e for edge):
% 'all','vv','vf','ve','ev','ef','ee','fv','ff','fe'

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
markerSize=50;

%% Example 1: Demonstrating |patchConnectivity| for different patch types

for testCase=1:3

    switch testCase
        case 1
            [F,V]=geoSphere(1,1);
        case 2
            [F,V]=quadSphere(2,1);
        case 3
            r=1; %Sphere radius
            rc=1.5; %Central radius
            nr=12;
            nc=18;
            patchType='honey';
            [F,V]=patchTorus(r,nr,rc,nc,patchType);
    end

    %%
    % Using |patchConnectivity| to compute connectivity arrays

    [C]=patchConnectivity(F,V);

    %%
    E=C.edge.vertex;
    Z=V(:,3);
    [~,indPlotVertex]=max(Z);
    [~,indPlotFace]=max(mean(Z(F),2));
    [~,indPlotEdge]=max(mean(Z(E),2));

    %%
    % Visualize vertex connectivity

    indVertexVertex=C.vertex.vertex(indPlotVertex,:);
    indVertexVertex=indVertexVertex(indVertexVertex>0);
    indVertexFace=C.vertex.face(indPlotVertex,:);
    indVertexEdge=C.vertex.edge(indPlotVertex,:);

    cFigure;
    subplot(1,3,1); hold on;
    title('Vertex-vertex connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=plotV(V(indPlotVertex,:),'g.','MarkerSize',markerSize);
    hl(2)=plotV(V(indVertexVertex,:),'r.','MarkerSize',markerSize);

    legend(hl,{'Example point','Connected points'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    subplot(1,3,2); hold on;
    title('Vertex-face connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=plotV(V(indPlotVertex,:),'g.','MarkerSize',markerSize);
    hl(2)=gpatch(F(indVertexFace,:),V,'r','r',1);

    legend(hl,{'Example point','Connected faces'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    subplot(1,3,3); hold on;
    title('Vertex-edge connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=plotV(V(indPlotVertex,:),'g.','MarkerSize',markerSize);
    hl(2)=gpatch(E(indVertexEdge,:),V,'none','r',1);
    hl(2).LineWidth=3;
    legend(hl,{'Example point','Connected Edges'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    drawnow;

    %%
    % Visualize face connectivity

    indFaceVertex=C.face.vertex(indPlotFace,:);
    indFaceVertex=indFaceVertex(indFaceVertex>0);
    indFaceFace=C.face.face(indPlotFace,:);
    indFaceFace=indFaceFace(indFaceFace>0);
    indFaceEdge=C.face.edge(indPlotFace,:);

    cFigure;
    subplot(1,3,1); hold on;
    title('Face-vertex connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=gpatch(F(indPlotFace,:),V,'g','g',1);
    hl(2)=plotV(V(indFaceVertex,:),'r.','MarkerSize',markerSize);

    legend(hl,{'Example face','Connected points'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    subplot(1,3,2); hold on;
    title('Face-face connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=gpatch(F(indPlotFace,:),V,'g','g',1);
    hl(2)=gpatch(F(indFaceFace,:),V,'r','r',1);

    legend(hl,{'Example face','Connected faces'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    subplot(1,3,3); hold on;
    title('Face-edge connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=gpatch(F(indPlotFace,:),V,'g','g',1);
    hl(2)=gpatch(E(indFaceEdge,:),V,'none','r',1);
    hl(2).LineWidth=3;
    legend(hl,{'Example face','Connected Edges'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    drawnow;

    %%
    % Visualize edge connectivity

    indEdgeVertex=C.edge.vertex(indPlotEdge,:);
    indEdgeVertex=indEdgeVertex(indEdgeVertex>0);
    indEdgeFace=C.edge.face(indPlotEdge,:);
    indEdgeFace=indEdgeFace(indEdgeFace>0);
    indEdgeEdge=C.edge.edge(indPlotEdge,:);
    indEdgeEdge=indEdgeEdge(indEdgeEdge>0);

    cFigure;
    subplot(1,3,1); hold on;
    title('Edge-vertex connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=gpatch(E(indPlotEdge,:),V,'g','g',1); hl(1).LineWidth=3;
    hl(2)=plotV(V(indEdgeVertex,:),'r.','MarkerSize',markerSize);

    legend(hl,{'Example edge','Connected points'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    subplot(1,3,2); hold on;
    title('Edge-face connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=gpatch(E(indPlotEdge,:),V,'g','g',1); hl(1).LineWidth=3;
    hl(2)=gpatch(F(indEdgeFace,:),V,'r','r',1);

    legend(hl,{'Example edge','Connected faces'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    subplot(1,3,3); hold on;
    title('Edge-edge connectivity');
    gpatch(F,V,'kw','k',0.5)

    hl(1)=gpatch(E(indPlotEdge,:),V,'g','g',1); hl(1).LineWidth=3;
    hl(2)=gpatch(E(indEdgeEdge,:),V,'none','r',1); hl(2).LineWidth=3;
    legend(hl,{'Example edge','Connected Edges'},'Location','SouthOutside');

    axisGeom;
    view(2);
    camlight headlight;
    axis off;
    clear hl;

    drawnow;

end

%% Example 2: Using the conTypes input to control output request

%%
% Request only the edge-edge connectivity
C=patchConnectivity(F,V,'ee')

%%
% Request both the edge-edge connectivity and the face-face connectivity
C=patchConnectivity(F,V,{'ee','ff'})

%%
% Loop over all types and compare computation speed

conTypeSet={'all','vv','vf','ve','ev','ef','ee','fv','ff','fe'};
for q=1:1:numel(conTypeSet)
    tic
    C=patchConnectivity(F,V,conTypeSet{q});
    t=toc;
    disp([conTypeSet{q},': ',num2str(t),' seconds'])
end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%Parametric ReLU function
x = -10:0.01:10;
a = input('Enter the value of the slope')
y = @(x) (x).*((x > 0)) + (a*x).*((x < 0)) ;
fplot(y);
xlabel('x');
ylabel('y');
grid on

author: TheAlgorithms

repo: MATLAB-Octave


text:
function javaFigs = mrvJavaFeature(javaFigs)
% Turn off javafigures in certain versions of Matlab
%
%    javaFigs = mrvJavaFeature(javaFigs)
%
% Rory turned off Java Figure features in versions 7 - 7.3 because java was
% buggy for Mathworks.  Starting in 7.4 (Matlab 2007a), Mathworks disabled
% this possibility. THey now think that java works great and they want you
% to report any bugs.  So, we only disable for the relevant versions.
%
% There are many places in the code where Rory did this.   Rather than
% finding them all and writing out this little bit, I wrote this routine
% which we should use to replace Rory's test code.  In the fullness of
% time, as everyone migrates beyond Matlab 7.4, we can eliminate these
% calls, too.
%
% Example - You need to alter code like this:
%
%  javaFigs = mrvJavaFeature;
%
% .... code here
%
%  mrvJavaFeature(javaFigs);
% 2007 by BAW.
warning('You are calling an obsolote function: mrvJavaFeature. Please remove it from your code.');
javaFigs = [];
return;

% If the person sent in a javaFigs variable, they want us to restore
if notDefined('javaFigs'), javaFigs = [];
else feature('javafigures',javaFigs); return;
end

% Otherwise, they want us to turn it off for certain versions of Matlab,
% specified here
version = ver('Matlab');
matlabVersion = version.Version;
mVersion = str2double(matlabVersion(1:3));
mMinorVersion = str2double(matlabVersion(3:end));

if ( (mVersion == 7) && (mMinorVersion < 4)) % version 7, < 7.4
    javaFigs = feature('javafigures');
    if ispref('VISTA', 'javaOn') % seems to work diff't for diff't machines
        feature('javafigures', getpref('VISTA', 'javaOn'));
    else
        feature('javafigures', 0);
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


function minus_log_dnsty = uc_opt_hyperpara(hyperpara,y,lags,options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computes the likelihood of the US system
% a(t)  = a1 a(t-1) + ... + ap a(t-p) + ea(t )    [transition 1]
% b(t)  = c(t-1)    + b(t-1) + eb(t);             [transition 2]
% c(t)  = c(t-1)             + ec(t);             [transition 3]
% y(t)  = a(t)      + b(t);                       [transistion 4]

% to a state space of the form
% x(t) = A x(t-1) + B Sigma' u(t) ~ N(0,I)
% y(t) = C*(cons + x(t-1))
% where A is the companion form of the lag struture

% Filippo Ferroni, 6/1/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% options.state_space_model = 2;
%
% options.only_logL    = 1;
% options.initialCond  = 2;
% options.aZero        = zeros(lags+2+1,1);
% options.aZero(end-2) = y(1);
% options.pZero        = 10*eye();

hyperpara(options.index_est)   = exp(hyperpara);
if isempty(options.index_fixed) == 0
    hyperpara(options.index_fixed) = exp(options.hyperpara_fidex);
end

phi   = (hyperpara(1:lags));
sigma = diag((hyperpara(1+lags:3+lags)));

log_dnsty       = kfilternan(phi,sigma,y,options);
minus_log_dnsty = -log_dnsty;


text:
function fem2d_pack_test20 ( )

%*****************************************************************************80
%
%% TEST20 tests SPHERE_GRID_Q9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  element_order = 9;
  nelemx = 3;
  nelemy = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST20\n' );
  fprintf ( 1, '  SPHERE_GRID_Q9_ELEMENT sets up a grid of\n' );
  fprintf ( 1, '    Q9 quadrilaterals on a sphere.\n' );
  fprintf ( 1, '  SPHERE_GRID_Q9_ELEMENT_NUM returns the number\n' );
  fprintf ( 1, '    of elements in the grid\n' );
  fprintf ( 1, '  SPHERE_GRID_Q9_NODE_NUM returns the number\n' );
  fprintf ( 1, '    of nodes in the grid.\n' );
  fprintf ( 1, '  SPHERE_GRID_Q9_NODE_XYZ returns the coordinates\n' );
  fprintf ( 1, '    of nodes in the grid.\n' );

  element_num = sphere_grid_q9_element_num ( nelemx, nelemy );
  node_num = sphere_grid_q9_node_num ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected number of nodes =    %d\n', node_num );
  fprintf ( 1, '  Expected number of elements = %d\n', element_num );

  element_node = sphere_grid_q9_element ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The elements and their nodes:\n' );
  fprintf ( 1, '\n' );

  for element = 1 : element_num
    fprintf ( 1, '%4d  ', element );
    for order = 1 : element_order
      fprintf ( 1, '%4d', element_node(order,element) );
    end
    fprintf ( 1, '\n' );
  end

  node_xyz = sphere_grid_q9_node_xyz ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The node coordinates:\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num
    fprintf ( 1, '  %4d  %12f  %12f  %12f\n', node, node_xyz(1:3,node) );
  end
%
%  Write the elements and nodes to files.
%
  r8mat_write ( 'sphere_q9_nodes.txt', 3, node_num, node_xyz );

  i4mat_write ( 'sphere_q9_elements.txt', element_order, element_num, ...
    element_node );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function clusters = xyz2clusters(xyz,P)
% Converts a 3-column x, y, z list of mm coordinates to a clusters
% structure given P, the filename of an analyze .img file to provide
% dimensions and voxel sizes.
%
% :Usage:
% ::
%
%    function cl = xyz2clusters(xyz,P)
%
% Uses this info from the image:
%   - VOL.M     - spm-style mat matrix
%   - VOL.VOX   - voxel sizes
%
%   - SPM.Z     - now 1s; could stores values in the original image in clusters.Z
%
% The following is created internally:
%   - SPM.XYZmm - mm coords, you input these
%   - SPM.XYZ   - voxel coords
%

V = spm_vol(P);

%wh = strmatch('Thalamus',L3); whos wh

VOL.M = V.mat;
VOL.VOX = diag(V.mat(1:3,1:3));

fprintf(1,'Converting mm to voxels -> ');
XYZ = mm2voxel(xyz,VOL,2);  % unique reordered
fprintf(1,'%3.0f voxels -> ',size(XYZ,1));
fprintf(1,'Converting voxels to mm -> ');
%XYZ = unique(XYZ,'rows');
XYZ = XYZ';

SPM.XYZmm = voxel2mm(XYZ,VOL.M);
SPM.XYZ = XYZ;
SPM.Z = ones(1,size(XYZ,2));

fprintf(1,'Making clusters.\n');

%cl = tor_extract_rois([],SPM,VOL);

% stuff from tor_extract_rois -- omit where it gets each region separately

    % ----------------------------------------------------------------------------------
    % define each cluster as cell in array.
    % ----------------------------------------------------------------------------------
    clusters = [];
cl_index = spm_clusters(SPM.XYZ);
%cl_index = SPM.Z;    %:size(SPM.XYZ,2);
cl.title = 'xyz';

    for i = 1:max(cl_index)
      %  if verbose, fprintf(1,['\nExtracting cluster ' num2str(i) ' of ' num2str(max(cl_index))]),end
        a = find(cl_index == i);

        %if isfield(SPM,'u'), cl.title = SPM.title;,else, cl.title = 'Untitled';,end
        %if isfield(SPM,'u'), cl.threshold = SPM.u;,else, cl.threshold = NaN;,end
        cl(i).voxSize = VOL.VOX;
        cl(i).M = VOL.M;
        cl(i).name = [cl.title '_' num2str(i) '_' mat2str(size(a,2)) '_voxels'];
        cl(i).numVox = size(a,2);
        cl(i).Z = SPM.Z(a);
        cl(i).XYZmm = SPM.XYZmm(:,a);
        cl(i).XYZ = SPM.XYZ(:,a);

        cl(i).center = mean(cl(i).XYZ');
        if size(cl(i).XYZmm,2) > 1, cl(i).mm_center = center_of_mass(cl(i).XYZmm,cl(i).Z);  % mean(cl.XYZmm');
        else cl(i).mm_center = cl(i).XYZmm';
        end
        clusters = [clusters, cl];

    end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ n_data_new, n, c ] = bernoulli_number_values ( n_data )

%*****************************************************************************80
%
%% BERNOULLI_NUMBER_VALUES returns some values of the Bernoulli numbers.
%
%  Discussion:
%
%    The Bernoulli numbers are rational.
%
%    If we define the sum of the M-th powers of the first N integers as:
%
%      SIGMA(M,N) = sum ( 0 <= I <= N ) I**M
%
%    and let C(I,J) be the combinatorial coefficient:
%
%      C(I,J) = I! / ( ( I - J )! * J! )
%
%    then the Bernoulli numbers B(J) satisfy:
%
%      SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) C(M+1,J) B(J) * (N+1)**(M+1-J)
%
%  First values:
%
%   B0  1                   =         1.00000000000
%   B1 -1/2                 =        -0.50000000000
%   B2  1/6                 =         1.66666666666
%   B3  0                   =         0
%   B4 -1/30                =        -0.03333333333
%   B5  0                   =         0
%   B6  1/42                =         0.02380952380
%   B7  0                   =         0
%   B8 -1/30                =        -0.03333333333
%   B9  0                   =         0
%  B10  5/66                =         0.07575757575
%  B11  0                   =         0
%  B12 -691/2730            =        -0.25311355311
%  B13  0                   =         0
%  B14  7/6                 =         1.16666666666
%  B15  0                   =         0
%  B16 -3617/510            =        -7.09215686274
%  B17  0                   =         0
%  B18  43867/798           =        54.97117794486
%  B19  0                   =         0
%  B20 -174611/330          =      -529.12424242424
%  B21  0                   =         0
%  B22  854,513/138         =      6192.123
%  B23  0                   =         0
%  B24 -236364091/2730      =    -86580.257
%  B25  0                   =         0
%  B26  8553103/6           =   1425517.16666
%  B27  0                   =         0
%  B28 -23749461029/870     = -27298231.0678
%  B29  0                   =         0
%  B30  8615841276005/14322 = 601580873.901
%
%  Recursion:
%
%    With C(N+1,K) denoting the standard binomial coefficient,
%
%    B(0) = 1.0
%    B(N) = - ( sum ( 0 <= K < N ) C(N+1,K) * B(K) ) / C(N+1,N)
%
%  Special Values:
%
%    Except for B(1), all Bernoulli numbers of odd index are 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA_NEW
%    from the previous call.
%
%    Output, integer N_DATA_NEW, the index of the test data.
%
%    Output, integer N, the order of the Bernoulli number.
%
%    Output, real C, the value of the Bernoulli number.
%
  n_max = 10;
  c_vec = [ ...
    1.0000000000E+00, -0.5000000000E+00,  0.1666666667E+00, ...
    0.0000000000E+00, -0.0333333333E+00, -0.02380952381E+00, ...
   -0.0333333333E+00,  0.0757575757E+00, -529.1242424E+00, ...
    601580873.9E+00 ];
  n_vec = [ ...
     0,  1,  2,  3,  4, 6,  8, 10, 20, 30 ];

  n_data_new = n_data;

  if ( n_data_new < 0 )
    n_data_new = 0;
  end

  n_data_new = n_data_new + 1;

  if ( n_max < n_data_new )
    n_data_new = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data_new);
    c = c_vec(n_data_new);
  end

  return
end


text:
function f = div( f )
%DIV   Divergence of a CHEBFUN2V.
%   DIV(F) returns the divergence of the CHEBFUN2V i.e.,
%       divergence(F) = F_x + F_y.
%
%  This is shorthand for the command DIVERGENCE.
%
% See also DIVERGENCE.

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

% Note that divergence of a 3-vector is the same, because the functions are
% of two variables.

f = divergence( f );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [overlap, only1, only2] = calculate_overlap(T1, T2, bounds)
% calculate_overlap  Calculates overlap for two trajectories
%
% The function calculates per-frame overlap between two trajectories. Besides the
% region overlap the function also returns cointainment of the second trajectory in the
% first and containment of the first trajectory in the second.
%
% If the trajectories are not of equal length, then the overlaps are calculated up
% to the end of the shorter one.
%
% Input:
% - T1 (cell): The first trajectory.
% - T2 (cell): The second trajectory.
% - bounds (vector): An optional bounds of valid region where the overlap is calculated.
%
% Output:
% - overlap: A vector of per-frame overlaps.
% - only1: A vector of per-frame containment of the second trajectory in the first.
% - only2: A vector of per-frame containment of the first trajectory in the second.

len = min(size(T1, 1), size(T2, 1));
T1 = T1(1:len, :);
T2 = T2(1:len, :);

if (~iscell(T1))
    T1 = num2cell(T1, 2);
end

if (~iscell(T2))
    T2 = num2cell(T2, 2);
end

if nargin < 3
    bounds = [];
end

if get_global_variable('legacy_rasterization', true)
    mode = 'legacy';
else
    mode = 'default';
end;

results = region_overlap(T1, T2, bounds, mode);

%results = cell2mat(cellfun(@(r1, r2) region_overlap(r1, r2), T1, T2, 'UniformOutput', false));

overlap = results(:, 1);
only1 = results(:, 2);
only2 = results(:, 3);

author: votchallenge

repo: toolkit-legacy


text:
function mec_weighting = lmsSR_computeWeightsMECv4(mec_confMapN, mec_thr_scaler, mec_weight_exp)
%
% LMSSR_COMPUTEWEIGHTSMECV4 Computes weights depending on the quality of the motion estimation/compensation stage. (Version 4: One global threshold, adapted by number of images)
%    mec_weighting = LMSSR_COMPUTEWEIGHTSMECV4(mec_confMapN, mec_thr_scaler, mec_weight_exp)
%
% Parameters: mec_confMapN     -   3-D matrix containing confidence information in form of SSD values for each pixel and each auxiliary frame
%             mec_thr_scaler   -   Scalar value responsible for scaling the standard deviation threshold which cuts off SSD values that are too large [default: 20]
%             mec_weight_exp   -   Exponent for intensifying the weighting [default: 2]
%
%
% Author: Michel Bätz (LMS)
%
% See also: lmsSR_framework
%

[height,width,num_imgs] = size(mec_confMapN);

% 1-SSD Weighting + Setting SSD values larger than twice the standard deviation to 0
mec_weighting   = ones(height,width,num_imgs+1); % Highest weight is set to 1 => Reference frame weights need no further adaptation
mec_conf_stddev = std(mec_confMapN(:));
mec_weighting(:,:,2:end) = 1 - (mec_confMapN / (mec_thr_scaler*(1/num_imgs)*mec_conf_stddev)); % SSD values larger than threshold_scaler times standard deviation are set to 0 [default: 20]
mec_weighting(mec_weighting < 0) = 0;

%mec_weighting(mec_weighting > 0) = 1; % QUICK-TEST
mec_weighting = mec_weighting.^mec_weight_exp; % Default: 2

author: thomas-koehler

repo: SupER


text:
% DEMOIL3 Oil data with deterministic training conditional.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'oil';
experimentNo = 3;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('dtc');
options.optimiser = 'scg';
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

if exist('printDiagram') & printDiagram
  fgplvmPrintPlot(model, lbls, capName, experimentNo);
end

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

% compute the nearest neighbours errors in latent space.
errors = fgplvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ n_data, n, c ] = phi_values ( n_data )

%*****************************************************************************80
%
%% PHI_VALUES returns some values of the PHI function.
%
%  Discussion:
%
%    PHI(N) is the number of integers between 1 and N which are
%    relatively prime to N.  I and J are relatively prime if they
%    have no common factors.  The function PHI(N) is known as
%    "Euler's totient function".
%
%    By convention, 1 and N are relatively prime.
%
%    In Mathematica, the function can be evaluated by:
%
%      EulerPhi[n]
%
%  First values:
%
%     N  PHI(N)
%
%     1    1
%     2    1
%     3    2
%     4    2
%     5    4
%     6    2
%     7    6
%     8    4
%     9    6
%    10    4
%    11   10
%    12    4
%    13   12
%    14    6
%    15    8
%    16    8
%    17   16
%    18    6
%    19   18
%    20    8
%
%  Formula:
%
%    PHI(U*V) = PHI(U) * PHI(V) if U and V are relatively prime.
%
%    PHI(P**K) = P**(K-1) * ( P - 1 ) if P is prime.
%
%    PHI(N) = N * Product ( P divides N ) ( 1 - 1 / P )
%
%    N = Sum ( D divides N ) PHI(D).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the argument of the PHI function.
%
%    Output, integer C, the value of the PHI function.
%
  n_max = 20;

  c_vec = [ ...
      1,   1,   2,   2,   4,   2,   6,   4,   6,   4, ...
      8,   8,  16,  20,  16,  40, 148, 200, 200, 648 ];

  n_vec = [ ...
      1,   2,   3,   4,   5,   6,   7,   8,   9,  10, ...
     20,  30,  40,  50,  60, 100, 149, 500, 750, 999 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data);
    c = c_vec(n_data);
  end

  return
end


text:
function [lb_e, sigma_e, tf_conv] = estimate_MP2(s, pout, gap,...
    niter, nx_target, flag_distance, tf_verbose)
% fitting with histogram, using L2 or dKL
if ~exist('pout', 'var')
    pout = 0.5;
end
if ~exist('gap', 'var')
    gap = 0.2;
end
if ~exist('niter', 'var')
    niter = 10;
end
if ~exist('tf_verbose', 'var')
    tf_verbose = false;
end
if ~exist('nx_target', 'var')
    nx_target = 3.2;
end

if ~exist('flag_distance', 'var')
    flag_distance = 'L2';
end

ep = 1e-4;
tol = 1e-3;
se = sort(s(s>ep),'descend');

% start with the moment estimator
[lb_e, sigma_e, ~] = estimate_MP(s, pout, gap, 50, false);
if tf_verbose; display(lb_e, sigma_e); end
tf_conv = false;
for i = 1: niter-1
    lb1_e = (1 + sqrt(lb_e))^2 * sigma_e^2;
    nbulk = sum(se < lb1_e);
    se_bulk = se(se < lb1_e);
    nbin = nbulk / nx_target;
    nbin = max(10, nbin);
    [Ns, bin_edges] = histcounts(se_bulk, ceil(nbin));

    bin_centers = (bin_edges(1:end-1) + bin_edges(2:end)) / 2;
    s_pdf = Ns ./ diff(bin_edges) / nbulk;

    switch flag_distance
        case 'L2'
            f_dist = @(p2) sum((s_pdf - MPdistr(bin_centers, p2(1), p2(2))).^2);
        case 'dKL'
            f_dist = @(p2) f_dist_dKL(p2, s_spdf, bin_centers);
    end;

    [p2_fit, ~] = fminsearch(f_dist, [lb_e, sigma_e], optimset('MaxFunEvals',1e5));
    lb_e_new = p2_fit(1);
    sigma_e_new = p2_fit(2);

    if abs(sigma_e_new - sigma_e) < tol && abs(lb_e_new-lb_e) < tol
        tf_conv = true;
        break;
    else
        lb_e = lb_e_new;
        sigma_e = sigma_e_new;
        if tf_verbose; display([lb_e, sigma_e]); end
    end
end
end

% Local functions
% not yet updated...
function d=f_dist_dKL(p2,nf,bin_centers)
nf_MP=prob_MP_bin(p2(1),p2(2),bin_centers);
d_list=nf_MP.*log(nf_MP./nf);
d_list(nf_MP==0)=0;
d=sum(d_list);
% % special treatment of nf=0 bins near 0
if isinf(d)
    d=1e5+sum((nf-prob_MP_bin(p2(1),p2(2),bin_centers)).^2);
end;
end

author: xiuyechen

repo: FishExplorer


text:
function out = PeubChannel(P, n, L)
%upper bound for the error probability of the AWGN channel
%P - channel SNR
%n - block length (scalar)
%L - length of messages transmitter between the source encoder and the channel encoder

%
%   Created in 2012 by Victoria Kostina (vkostina@caltech.edu)
%

%dP/dQ is bounded by a constant. Find that constant.
[~, gamma2] = fminbnd(@(x)-ncx2pdf(n*x, n, n*P)./chi2pdf(n*x/(1+P), n)*(1+P), 0, 10);
gamma2 = -gamma2;

A = n/2*log(1+P) + n/2 - L - log(gamma2);

Finner = @(v) ncx2pdf(n*v, n, n/P).*n.*exp(P/2/(1+P)*n*(v - thresv()));
tail = 1 - ncx2cdf(thresv()*n,n,n/P);
out = quad(Finner, 0, thresv())...                      %exponent < 1
    + tail;                                             %exponent = 1

    function out = thresv()
        %threshold for v as a function of v0
        out = max( 0, 2*(1+P)/P/n*A );
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function value = year_is_embolismic_hebrew ( y )

%*****************************************************************************80
%
%% YEAR_IS_EMBOLISMIC_HEBREW returns TRUE if the Hebrew year was embolismic.
%
%  Discussion:
%
%    In a 19 year cycle, there are 7 embolismic years.  During these years,
%    an extra month, "Adar II", (sometimes called "Veadar") is inserted after
%    the month of Adar.  Nonembolismic years are called "common" years.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, logical VALUE, TRUE if the year was embolismic.
%
  if ( 12 <= i4_modp ( 7 * y + 13, 19 ) )
    value = 1;
  else
    value = 0;
  end

  return
end


text:
function SyntaxAnalysis_Callback(hObject, eventdata, handles)

[filename, filepath] = uigetfile({'*.mat;*.xlsx'},'Select Detection Files OR Exported Detections','MultiSelect', 'on',char(handles.data.settings.detectionfolder));
if isnumeric(filename); return; end
filename = cellstr(filename);

settings = inputdlg({'Maximum bout seperation (s)','Exclude Classes with frequency below (0-1)'},'Syntax',[1 50],{'2','.01'});
boutlength = str2num(settings{1});
minfreq = str2num(settings{2});

[~,~,ext] = fileparts(filename{1});

h = waitbar(0,'Loading Files');
%% Load Files
AllCalls = table([],[],[],[],[],'VariableNames',{'BeginTime_s_','Label','File','Bout','Accepted'});

if strcmp(ext,'.mat')
    for i = 1:length(filename)
        Calls = loadCallfile([filepath filename{i}],handles);

        Calls = Calls(Calls.Accept == 1, :);
        begintime = Calls.Box(:,1); % Get Box Position
        CallClass = Calls.Type;
        dist = pdist2(begintime,begintime);
        dist(dist > boutlength) = 0;
        G = graph(dist);
        bout = conncomp(G)';

        try
            g = table(begintime, CallClass, repmat(i,height(Calls),1), bout, Calls.Accept,'VariableNames',{'BeginTime_s_','Label','File','Bout','Accepted'});
            AllCalls = [AllCalls; g];
            waitbar(i/length(filename),h)
        catch
            warning('No Accepted Calls in File');
        end

    end
elseif strcmp(ext,'.xlsx')
    for i = 1:length(filename)
        t = readtable([filepath filename{i}]);
        t = t(t.Accepted == 1,{'Accepted' 'BeginTime_s_','Label'});
        dist = pdist2(t.BeginTime_s_,t.BeginTime_s_);
        dist(dist > boutlength) = 0;
        G = graph(dist);
        bout = conncomp(G)';
        try
            g = table(bout,repmat(i,height(t),1),'VariableNames',{'Bout','File'});
            AllCalls = [AllCalls; g t];
            waitbar(i/length(filename),h)
        catch
            warning('No Accepted Calls in File');
        end

    end
end
AllCalls.Label = categorical(AllCalls.Label);
AllCalls = sortrows(AllCalls,{'File' 'BeginTime_s_'});

%% Exclude rare classes and rejected calls
catcounts = countcats(AllCalls.Label);
CommonClasses = (catcounts ./ sum(catcounts)) > minfreq;
cats = categories(AllCalls.Label);
goodCalls = any(AllCalls.Label == cats(CommonClasses)',2);
AllCalls = AllCalls(AllCalls.Accepted == 1 & goodCalls ,:);
AllCalls.Label = removecats(AllCalls.Label);
cats = categories(AllCalls.Label);

%% Create Transition Table
counts = [];
for i = 1:length(cats)
    for j = 1:length(cats)
        %         counts(j,i) = mean((AllCalls{1:end-1,'Label'} == cats{i}) & (AllCalls{2:end,'Label'} == cats{j}));
        x = sum(...
            ((AllCalls{1:end-1,'Label'} == cats{i}) & (AllCalls{2:end,'Label'} == cats{j})) & ...
            (AllCalls{1:end-1,'Bout'} == AllCalls{2:end,'Bout'}) &...
            (AllCalls{1:end-1,'File'} == AllCalls{2:end,'File'}));
        n = sum(AllCalls{2:end,'Label'} == cats{j});
        counts(j,i) = x / n;
    end
end
delete(h)

%% Create Figures
figure('position',[0 0 700 600],'color','w')
h = heatmap(cats,cats,counts,'ColorMethod','mean');
h.XLabel = 'Transition Probability';
h.YLabel = 'Syllable';
set(gcf,'Colormap',inferno);
set(gca,'GridVisible','off','FontSize',14);
h.CellLabelFormat = '%.3f';
colorbar off

figure('position',[300 0 700 600],'color','w')
counts(counts<.05)=0;
G = digraph(counts,cats);
idxg=find(G.Edges.Weight<.075);
G=rmedge(G,idxg);
g = plot(G,'EdgeCData',G.Edges.Weight);
axis off
g.ArrowSize = 15;
g.LineWidth = 15*G.Edges.Weight;
g.EdgeAlpha = .8;
if isfield(g, 'NodeFontSize')
    g.NodeFontSize = 12;
end
g.NodeColor = [.2 .2 .2];
layout(g,'circle');
set(gcf,'Colormap',inferno);
cb=colorbar;
cb.Label.String='Transition Probability';
cb.Label.Rotation=90;

%% Save Matrix
[file,path] = uiputfile('*.xlsx','Save Transion Matrix');

if ~isnumeric(file)
    output = [];
    counts = cell(length(cats));
    for i = 1:length(cats)
        for j = 1:length(cats)
            %         counts(j,i) = mean((AllCalls{1:end-1,'Label'} == cats{i}) & (AllCalls{2:end,'Label'} == cats{j}));
            x = sum(...
                ((AllCalls{1:end-1,'Label'} == cats{i}) & (AllCalls{2:end,'Label'} == cats{j})) & ...
                (AllCalls{1:end-1,'Bout'} == AllCalls{2:end,'Bout'}) &...
                (AllCalls{1:end-1,'File'} == AllCalls{2:end,'File'}));
            n = sum(AllCalls{2:end,'Label'} == cats{j});
            counts(j,i) = {x / n};
        end
    end
    output = [output; [{'Conditional Probability'} cell(1,length(cats))]];
    output = [output; [[{'Category'} cats'] ; [cats, counts]]];

    counts = cell(length(cats));
    for i = 1:length(cats)
        for j = 1:length(cats)
            %         counts(j,i) = mean((AllCalls{1:end-1,'Label'} == cats{i}) & (AllCalls{2:end,'Label'} == cats{j}));
            x = sum(...
                ((AllCalls{1:end-1,'Label'} == cats{i}) & (AllCalls{2:end,'Label'} == cats{j})) & ...
                (AllCalls{1:end-1,'Bout'} == AllCalls{2:end,'Bout'}) &...
                (AllCalls{1:end-1,'File'} == AllCalls{2:end,'File'}));
            counts(j,i) = {x };
        end
    end
    output = [output; cell(3,length(cats)+1)];
    output = [output; [{'Transition Count'} cell(1,length(cats))]];
    output = [output;[[{'Category'} cats'] ; [cats, counts]]];

    counts = cell(length(cats),1);
    for i = 1:length(cats)
        counts(i) = {sum(AllCalls{1:end,'Label'} == cats{i})};
    end

    output = [output; cell(3,length(cats)+1)];
    output = [output; [{'Total Count'} cell(1,length(cats))]];
    output = [output; cats, counts, cell(length(cats),length(cats)-1)];

    writetable(cell2table(output),[path file],'WriteVariableNames',0);
end

end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function C = times(A, B)
%*    CHEBOP multiplication.
%   C = A.*B, where ones of a A or B is a CHEBOP and the other a scalar
%   is equivalent to A*B.
%
%   All other instances returns an error.
%
% See also CHEBOP/MTIMES.

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

if ( isnumeric(A) || isnumeric(B) )
    C = mtimes(A, B);
else
    error('CHEBOP:TIMES:NotSupported', '%s.*%s is not supported.', ...
        class(A), class(B));
end

end


text:
function wevent = segment_event_waveforms(w, cobj, pretrigger, posttrigger, arrivalTimeCorrection)
%SEGMENT_EVENT_WAVEFORMS take a continuous vector of waveforms, and based
% arrival times in the infrasoundEvent vector, generate a cell vector where
% each element is a waveform vector extracted around the arrival times
% given
% Usage:
%     wevent = segment_event_waveforms(w, infrasoundEvent, pretrigger,
%                                      posttrigger,arrivalTimeCorrection)
% pretrigger and posttrigger are in seconds (e.g. 1)
if ~exist('arrivalTimeCorrection','var')
    arrivalTimeCorrection = 0.0;
end
disp('Segmenting event waveforms...')
numEvents = numel(infrasoundEvent);
for eventNumber=1:numEvents
    fprintf('- segmenting infrasound event %d of %d\n',eventNumber,numEvents );
    time1 = infrasoundEvent(eventNumber).FirstArrivalTime-pretrigger/86400-arrivalTimeCorrection/86400;
%     if arrivalTimeCorrection==0.0
%         time2 = infrasoundEvent(eventNumber).LastArrivalTime+posttrigger/86400-arrivalTimeCorrection/86400;
%     else
        time2 = infrasoundEvent(eventNumber).FirstArrivalTime+posttrigger/86400-arrivalTimeCorrection/86400;
%     end
    w2 = detrend(extract(w, 'time', time1, time2));
    wevent{eventNumber} = w2;
    clear w2
end

author: geoscience-community-codes

repo: GISMO


text:
function save_image(M, name, options)

% save_image - save an image or a set of images
%
%   save_image(M, name, options);
%
%   M and name can be cell arrays.
%
%   set options.clamp and options.rescale to clamp/rescale the image prior
%   to saving.
%
%   set options.format='bmp', 'jpg' ... to set the file format.
%
%   set options.base_str for the string that is appened to each filename.
%
%   Copyright (c) 2008 Gabriel Peyre

options.null = 0;
if iscell(M)
    if not(iscell(M)) || length(name)~=length(M)
        error('name should be a cell array as M');
    end
    for i=1:length(M)
        save_image(M{i}, name{i}, options);
    end
    return;
end

% trying to guess format
s = strfind(name, '.');
fmt = 'png';
if not(isempty(s))
    fmt = name(s(1)+1:end);
end
fmt = getoptions(options, 'format', fmt);

if iscell(fmt)
    for i=1:length(fmt)
        options.format = fmt{i};
        save_image(M, name, options);
    end
    return;
end

c = getoptions(options, 'clamp', []);
r = getoptions(options, 'rescale', []);
if c==1
    c = [0 1];
end
if r==1
    r = [0 1];
end
if not(isempty(c))
    M = clamp(M, c(1), c(2));
end
if not(isempty(r))
    M = rescale(M, r(1), r(2));
end
base_str = getoptions(options, 'base_str', '');

if strcmp(fmt, 'eps')
    warning('eps is unsuported for now');
    return;
    % detection of image magick
    global image_magick_path;
    if isempty(image_magick_path)
        image_magick_path = '/Applications/ImageMagick-6.4.3/bin/';
    end
    if not(exist([image_magick_path 'convert']))
        warning('Image magick not installed, you need to specify global image_magick_path');
    end

    warning off;
    imwrite(M, [base_str name '.png'],'png');
    warning on;
    system([ image_magick_path 'convert ./' base_str name '.png ./' base_str name '.eps' ]);
    delete( [base_str name '.png'] );
    return;
end

warning off;
imwrite(M, [base_str name '.' fmt], fmt);
warning on;

author: gpeyre

repo: matlab-toolboxes


text:
%DAGNN.UNIFORMSCALINGGRIDGENERATOR  Generate an iso-tropic scaling + translation
%   grid for bilinear resampling.
%   This layer maps 1 x 1 x 3 x N transformation parameters corresponding to
%   the scale and translation y and x respectively, to 2 x Ho x Wo x N
%   sampling grids compatible with dagnn.BlilinearSampler.

% (c) 2016 Ankush Gupta
classdef UniformScalingGridGenerator < dagnn.Layer

 properties
     Ho = 0;
     Wo = 0;
 end

  properties (Transient)
    % the grid --> this is cached
    % has the size: [2 x HoWo]
    xxyy ;
  end

  methods

    function outputs = forward(obj, inputs, ~)
      % input is a 1x1x3xN TENSOR corresponding to:
      % [  s 0 ty ]
      % [  0 s tx ]
      %
      % OUTPUT is a 2xHoxWoxN grid

      % reshape the tfm params into matrices:
      T = inputs{1};
      % check shape:
      sz_T = size(T);
      assert(all(sz_T(1:3) == [1 1 3]), 'transforms have incorrect shape');
      nbatch = size(T,4);
      S = reshape(T(1,1,1,:), 1,1,nbatch); % x,y scaling
      t = reshape(T(1,1,2:3,:), 2,1,nbatch); % translation
      % generate the grid coordinates:
      % use gpu?:
      useGPU = isa(T, 'gpuArray');
      if isempty(obj.xxyy)
        obj.initGrid(useGPU);
      end
      % transform the grid:
      g = bsxfun(@times, obj.xxyy, S); % scale
      g = bsxfun(@plus, g, t); % translate
      g = reshape(g, 2,obj.Ho,obj.Wo,nbatch);
      outputs = {g};
    end

    function [derInputs, derParams] = backward(obj, inputs, ~, derOutputs)
      %fprintf('restricted-affineGridGenerator backward\n');
      dY = derOutputs{1};
      useGPU = isa(dY, 'gpuArray');
      nbatch = size(dY,4);

      % create the gradient buffer:
      dA = zeros([1,1,3,nbatch], 'single');
      if useGPU, dA = gpuArray(dA); end

      dY  = reshape(dY, 2,obj.Ho*obj.Wo,nbatch);
      % gradient wrt the linear part:
      dA(1,1,1,:) = reshape(obj.xxyy,1,[]) * reshape(dY, [],nbatch);
      % gradient wrt translation (or bias):
      dA(1,1,2:3,:) = reshape(sum(dY,2),1,1,2,[]);

      derInputs = {dA};
      derParams = {};
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      nBatch = inputSizes{1}(4);
      outputSizes = {[2, obj.Ho, obj.Wo, nBatch]};
    end

    function obj = UniformScalingGridGenerator(varargin)
      obj.load(varargin);
      % get the output sizes:
      obj.Ho = obj.Ho;
      obj.Wo = obj.Wo;
      obj.xxyy = [];
    end

    function obj = reset(obj)
      reset@dagnn.Layer(obj) ;
      obj.xxyy = [] ;
    end

    function initGrid(obj, useGPU)
      % initialize the grid:
      % this is a constant
      xi = linspace(-1, 1, obj.Ho);
      yi = linspace(-1, 1, obj.Wo);
      [yy,xx] = meshgrid(yi,xi);
      xxyy = [xx(:), yy(:)]' ; % 2xM
      if useGPU
        xxyy = gpuArray(xxyy);
      end
      obj.xxyy = xxyy ;
    end

  end
end

author: jiangqy


text:
function [param] = get_parameters2(hdrf,t)
% Find model parameters
%
% Height - h
%
% Time to peak - p (in time units of TR seconds)
%
% Width (at half peak) - w
%
% Calculate Heights and Time to peak:
%
% delta = 1/(t(2)-t(1));

n = round(length(t)*0.8);
% n = round(t(end)*0.6*delta)

[~,p] = max(abs(hdrf(1:n)));
h = hdrf(p);

%if (p > t(end)*0.6*delta), warning('Late time to peak'), end;
if (p > t(end)*0.8), warning('Late time to peak'), end;

if (h >0)
    v = (hdrf >= h/2);
else
    v = (hdrf <= h/2);
end;

[~,b] = min(diff(v));
v(b+1:end) = 0;
w = sum(v);

cnt = p-1;
g =hdrf(2:end) - hdrf(1:(end-1));
while((cnt > 0) && (abs(g(cnt)) <0.001)),
    h = hdrf(cnt);
    p = cnt;
    cnt = cnt-1;
end;

param = zeros(3,1);
param(1) = h;
param(2) = p;
param(3) = w;

return;


text:
function child = children(tree,uid)
% XMLTREE/CHILDREN Return children's UIDs of node uid
% FORMAT child = children(tree,uid)
%
% tree   - a tree
% uid    - uid of the element
% child  - array of the UIDs of children of node uid
%__________________________________________________________________________
%
% Return UID's of children of node uid
%__________________________________________________________________________
% Copyright (C) 2002-2011  http://www.artefact.tk/

% Guillaume Flandin
% $Id: children.m 4460 2011-09-05 14:52:16Z guillaume $

%error(nargchk(2,2,nargin));

child = [];
uid = uid(:);
l = length(tree);
for i=1:length(uid)
    if uid(i) > 0 && uid(i) <= l
        if strcmp(tree.tree{uid(i)}.type,'element')
            child = [child tree.tree{uid(i)}.contents];
        end
    else
        error('[XMLTree] Invalid UID.');
    end
end
if isempty(child), child = []; end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ChannelMat = in_channel_eeglab_set(ChannelFile, isFixUnits)
% IN_CHANNEL_EEGLAB_SET:  Read 3D cartesian positions from an .SET EEGLAB file.
%
% USAGE:  ChannelMat = in_channel_eeglab_set(ChannelFile, isFixUnits=[ask])
%         ChannelMat = in_channel_eeglab_set(SetFileMat,  isFixUnits=[ask])
%
% INPUTS:
%     - ChannelFile  : Full path to the file (either .ela or .eps file)
%     - SetFileMat   : Structure contained in an EEGLAB .set file
%     - isFixUnit    : If 1: If the positions of the electrodes are normalized, try to convert them to something acceptable by brainstorm
%                            Otherwise, tries to convert the distance units to meters automatically
%                      If 0, does not apply  any correction
%                      If [], ask for use input

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

% Parse inputs
if (nargin < 2) || isempty(isFixUnits)
    isFixUnits = [];
end
% Load initial file (Matlab .mat format)
if ischar(ChannelFile)
    SetFileMat = load(ChannelFile, '-mat');
else
    SetFileMat = ChannelFile;
    ChannelFile = '';
end
% Convert Channel locations to Brainstorm format
if isfield(SetFileMat.EEG, 'chanlocs') && ~isempty(SetFileMat.EEG.chanlocs)
    % Check coordinates sytem
    isLoc = isfield(SetFileMat.EEG.chanlocs(1), 'X');
    if isLoc
        isNormalizedCs = (max([SetFileMat.EEG.chanlocs.X]) <= 1);
        if isempty(isNormalizedCs)
            isNormalizedCs = 0;
        elseif isNormalizedCs && isempty(isFixUnits)
            isFixUnits = java_dialog('confirm', ['The EEGLAB file you selected contains 3D electrodes positions, but they' 10 ...
                                               'seem to be spherical projections in a normalized coordinates system.' 10 ...
                                               'Brainstorm needs real 3D positions, you may need to import them separately.' 10 10 ...
                                               'Would you like Brainstorm to try to convert these positions ?'], 'EEGLAB electrodes positions');
        end
    end
    % Initialize returned structure
    nbChannels = SetFileMat.EEG.nbchan;
    ChannelMat = db_template('channelmat');
    ChannelMat.Comment = 'EEGLAB channels';

    % Copy channel data
    for iChan = 1:nbChannels
        % Electrode type and name
        if isfield(SetFileMat.EEG.chanlocs(iChan), 'type') && ~isempty(SetFileMat.EEG.chanlocs(iChan).type) && ischar(SetFileMat.EEG.chanlocs(iChan).type) && ismember(SetFileMat.EEG.chanlocs(iChan).type(1), 'abcdefghijklmnopqrstuvwyzABCDEFGHIJKLMNOPQRSTUVWXYZ')
            ChannelMat.Channel(iChan).Type = upper(SetFileMat.EEG.chanlocs(iChan).type);
        else
            ChannelMat.Channel(iChan).Type = 'EEG';
        end
        ChannelMat.Channel(iChan).Name = SetFileMat.EEG.chanlocs(iChan).labels;
        % Electrode location
        if isLoc
            if isempty(SetFileMat.EEG.chanlocs(iChan).X) || isempty(SetFileMat.EEG.chanlocs(iChan).Y) || isempty(SetFileMat.EEG.chanlocs(iChan).Z)
                ChannelMat.Channel(iChan).Loc = [];
                if strcmpi(ChannelMat.Channel(iChan).Type, 'EEG') && ~isempty([SetFileMat.EEG.chanlocs.X])
                    ChannelMat.Channel(iChan).Type = 'Misc';
                end
            elseif ~isNormalizedCs
                ChannelMat.Channel(iChan).Loc = [SetFileMat.EEG.chanlocs(iChan).X; ...
                                                 SetFileMat.EEG.chanlocs(iChan).Y; ...
                                                 SetFileMat.EEG.chanlocs(iChan).Z] ./ 1000;
            elseif isFixUnits
                ChannelMat.Channel(iChan).Loc = [SetFileMat.EEG.chanlocs(iChan).X / 9.4 + 0.007; ...
                                                 SetFileMat.EEG.chanlocs(iChan).Y / 11.5; ...
                                                 SetFileMat.EEG.chanlocs(iChan).Z / 8.7 + 0.042];
            else
                ChannelMat.Channel(iChan).Loc = [SetFileMat.EEG.chanlocs(iChan).X; ...
                                                 SetFileMat.EEG.chanlocs(iChan).Y; ...
                                                 SetFileMat.EEG.chanlocs(iChan).Z] / 10;
            end
        else
            ChannelMat.Channel(iChan).Loc = [];
        end
        ChannelMat.Channel(iChan).Orient  = [];
        ChannelMat.Channel(iChan).Comment = '';
        ChannelMat.Channel(iChan).Weight  = 1;
    end

    % Check distance units
    if isLoc && ~isNormalizedCs && ~isequal(isFixUnits, 0)
        if isempty(isFixUnits)
            isConfirmFix = 1;
        else
            isConfirmFix = 0;
        end
        ChannelMat = channel_fixunits(ChannelMat, 'mm', isConfirmFix);
    end
else
    ChannelMat = [];
end

author: brainstorm-tools


text:
% DEMOILFGPLVM9 Oil data with three dimensions and variational sparse approximation.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'oil';
experimentNo = 9;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('dtcvar');
options.kern = {'rbf', 'bias', 'whitefixed'};
options.optimiser = 'scg';
latentDim = 3;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);
model.kern.comp{3}.variance = 1e-4;
% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display them.
lvmScatterPlot(model, lbls,[],[1,2]); % default [1,2], change to plot
                                      % a different pair of lantent dims

% compute the nearest neighbours errors in latent space.
errors = fgplvmNearestNeighbour(model, lbls);


text:
% Simple algorithm to calculate the deflection, bending moment, shear force in a
% cantilever beam.
% load = loading distribution

% Example,
% 1.Uniform loading of 10 N/m: load=repmat(10,1,1001), l=20, dl=0.02
% 2.Point loading with P=400 N, load=[repmat(0,1,999) 10000 10000], l=20, dl=0.02

% length(load)=(l/dl)

function deflection=beam_deflection(load,ei,l,dl)
if length(load)~=(l/dl+1)
    error('Check inputs')
end
y=0:dl:l;
m=sum((y.*load))*dl;
v=sum(load)*dl;
u_4=load/ei;
u_3=v/ei;
for i=2:length(load)
    u_3(i)=u_3(i-1)-u_4(i-1)*dl;
end
u_2=m/ei;
for i=2:length(load)
    u_2(i)=u_2(i-1)-u_3(i-1)*dl;
end
u_1=0;
for i=2:length(load)
    u_1(i)=u_1(i-1)+u_2(i-1)*dl;
end
u=0;
for i=2:length(load)
    u(i)=u(i-1)+u_1(i-1)*dl;
end
deflection=u;
plot(y,u_2*ei)
hold on
plot(y,u_3*ei,'r')
legend('bending moment','shear force')
xlabel('length along the beam')
ylabel(' bending moment and shear force (SI units)')
grid
hold off
figure,plot(y,u,'r')
xlabel('length along the beam')
ylabel('deflection')
grid
title('Deflection')

author: Sable

repo: mcbench-benchmarks


text:
function value = p19_i4 ( action, name, value )

%*****************************************************************************80
%
%% P19_I4 sets or gets I4 parameters for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P19_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P19_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P19_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P19_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P19_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P19_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P19_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P19_I4 - Fatal error!' );

  end

  return
end


text:
%  Using 2D or 3D affine matrix to rotate, translate, scale, reflect and
%  shear a 2D image or 3D volume. 2D image is represented by a 2D matrix,
%  3D volume is represented by a 3D matrix, and data type can be real
%  integer or floating-point.
%
%  You may notice that MATLAB has a function called 'imtransform.m' for
%  2D spatial transformation. However, keep in mind that 'imtransform.m'
%  assumes y for the 1st dimension, and x for the 2nd dimension. They are
%  equivalent otherwise.
%
%  In addition, if you adjust the 'new_elem_size' parameter, this 'affine.m'
%  is equivalent to 'interp2.m' for 2D image, and equivalent to 'interp3.m'
%  for 3D volume.
%
%  Usage: [new_img new_M] = ...
%	affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);
%
%  old_img  -	original 2D image or 3D volume. We assume x for the 1st
%		dimension, y for the 2nd dimension, and z for the 3rd
%		dimension.
%
%  old_M  -	a 3x3 2D affine matrix for 2D image, or a 4x4 3D affine
%		matrix for 3D volume. We assume x for the 1st dimension,
%		y for the 2nd dimension, and z for the 3rd dimension.
%
%  new_elem_size (optional)  -  size of voxel along x y z direction for
%		a transformed 3D volume, or size of pixel along x y for
%		a transformed 2D image. We assume x for the 1st dimension
%		y for the 2nd dimension, and z for the 3rd dimension.
%		'new_elem_size' is 1 if it is default or empty.
%
%		You can increase its value to decrease the resampling rate,
%		and make the 2D image or 3D volume more coarse. It works
%		just like 'interp3'.
%
%  verbose (optional) - 1, 0
%		1:  show transforming progress in percentage
%		2:  progress will not be displayed
%		'verbose' is 1 if it is default or empty.
%
%  bg (optional)  -	background voxel intensity in any extra corner that
%		is caused by the interpolation. 0 in most cases. If it is
%		default or empty, 'bg' will be the average of two corner
%		voxel intensities in original data.
%
%  method (optional)  -	1, 2, or 3
%		1:  for Trilinear interpolation
%		2:  for Nearest Neighbor interpolation
%		3:  for Fischer's Bresenham interpolation
%		'method' is 1 if it is default or empty.
%
%  new_img  -	transformed 2D image or 3D volume
%
%  new_M  -	transformed affine matrix
%
%  Example 1 (3D rotation):
%	load mri.mat;   old_img = double(squeeze(D));
%	old_M = [0.88 0.5 3 -90; -0.5 0.88 3 -126; 0 0 2 -72; 0 0 0 1];
%	new_img = affine(old_img, old_M, 2);
%	[x y z] = meshgrid(1:128,1:128,1:27);
%	sz = size(new_img);
%	[x1 y1 z1] = meshgrid(1:sz(2),1:sz(1),1:sz(3));
%	figure; slice(x, y, z, old_img, 64, 64, 13.5);
%	shading flat; colormap(map); view(-66, 66);
%	figure; slice(x1, y1, z1, new_img, sz(1)/2, sz(2)/2, sz(3)/2);
%	shading flat; colormap(map); view(-66, 66);
%
%  Example 2 (2D interpolation):
%	load mri.mat;   old_img=D(:,:,1,13)';
%	old_M = [1 0 0; 0 1 0; 0 0 1];
%	new_img = affine(old_img, old_M, [.2 .4]);
%	figure; image(old_img); colormap(map);
%	figure; image(new_img); colormap(map);
%
%  This program is inspired by:
%  SPM5 Software from Wellcome Trust Centre for Neuroimaging
%	http://www.fil.ion.ucl.ac.uk/spm/software
%  Fischer, J., A. del Rio (2004). A Fast Method for Applying Rigid
%	Transformations to Volume Data, WSCG2004 Conference.
%	http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function [new_img, new_M] = affine(old_img, old_M, new_elem_size, verbose, bg, method)

   if ~exist('old_img','var') | ~exist('old_M','var')
      error('Usage: [new_img new_M] = affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);');
   end

   if ndims(old_img) == 3
      if ~isequal(size(old_M),[4 4])
         error('old_M should be a 4x4 affine matrix for 3D volume.');
      end
   elseif ndims(old_img) == 2
      if ~isequal(size(old_M),[3 3])
         error('old_M should be a 3x3 affine matrix for 2D image.');
      end
   else
      error('old_img should be either 2D image or 3D volume.');
   end

   if ~exist('new_elem_size','var') | isempty(new_elem_size)
      new_elem_size = [1 1 1];
   elseif length(new_elem_size) < 2
      new_elem_size = new_elem_size(1)*ones(1,3);
   elseif length(new_elem_size) < 3
      new_elem_size = [new_elem_size(:); 1]';
   end

   if ~exist('method','var') | isempty(method)
      method = 1;
   elseif ~exist('bresenham_line3d.m','file') & method == 3
      error([char(10) char(10) 'Please download 3D Bresenham''s line generation program from:' char(10) char(10) 'http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21057' char(10) char(10) 'to test Fischer''s Bresenham interpolation method.' char(10) char(10)]);
   end

   %  Make compatible to MATLAB earlier than version 7 (R14), which
   %  can only perform arithmetic on double data type
   %
   old_img = double(old_img);
   old_dim = size(old_img);

   if ~exist('bg','var') | isempty(bg)
      bg = mean([old_img(1) old_img(end)]);
   end

   if ~exist('verbose','var') | isempty(verbose)
      verbose = 1;
   end

   if ndims(old_img) == 2
      old_dim(3) = 1;
      old_M = old_M(:, [1 2 3 3]);
      old_M = old_M([1 2 3 3], :);
      old_M(3,:) = [0 0 1 0];
      old_M(:,3) = [0 0 1 0]';
   end

   %  Vertices of img in voxel
   %
   XYZvox = [	1		1		1
		1		1		old_dim(3)
		1		old_dim(2)	1
		1		old_dim(2)	old_dim(3)
		old_dim(1)	1		1
		old_dim(1)	1		old_dim(3)
		old_dim(1)	old_dim(2)	1
		old_dim(1)	old_dim(2)	old_dim(3)   ]';

   old_R = old_M(1:3,1:3);
   old_T = old_M(1:3,4);

   %  Vertices of img in millimeter
   %
   XYZmm = old_R*(XYZvox-1) + repmat(old_T, [1, 8]);

   %  Make scale of new_M according to new_elem_size
   %
   new_M = diag([new_elem_size 1]);

   %  Make translation so minimum vertex is moved to [1,1,1]
   %
   new_M(1:3,4) = round( min(XYZmm,[],2) );

   %  New dimensions will be the maximum vertices in XYZ direction (dim_vox)
   %  i.e. compute   dim_vox   via   dim_mm = R*(dim_vox-1)+T
   %  where, dim_mm = round(max(XYZmm,[],2));
   %
   new_dim = ceil(new_M(1:3,1:3) \ ( round(max(XYZmm,[],2))-new_M(1:3,4) )+1)';

   %  Initialize new_img with new_dim
   %
   new_img = zeros(new_dim(1:3));

   %  Mask out any changes from Z axis of transformed volume, since we
   %  will traverse it voxel by voxel below. We will only apply unit
   %  increment of mask_Z(3,4) to simulate the cursor movement
   %
   %  i.e. we will use   mask_Z * new_XYZvox   to replace   new_XYZvox
   %
   mask_Z = diag(ones(1,4));
   mask_Z(3,3) = 0;

   %  It will be easier to do the interpolation if we invert the process
   %  by not traversing the original volume. Instead, we traverse the
   %  transformed volume, and backproject each voxel in the transformed
   %  volume back into the original volume. If the backprojected voxel
   %  in original volume is within its boundary, the intensity of that
   %  voxel can be used by the cursor location in the transformed volume.
   %
   %  First, we traverse along Z axis of transformed volume voxel by voxel
   %
   for z = 1:new_dim(3)

      if verbose & ~mod(z,10)
         fprintf('%.2f percent is done.\n', 100*z/new_dim(3));
      end

      %  We need to find out the mapping from voxel in the transformed
      %  volume (new_XYZvox) to voxel in the original volume (old_XYZvox)
      %
      %  The following equation works, because they all equal to XYZmm:
      %  new_R*(new_XYZvox-1) + new_T  ==  old_R*(old_XYZvox-1) + old_T
      %
      %  We can use modified new_M1 & old_M1 to substitute new_M & old_M
      %      new_M1 * new_XYZvox       ==       old_M1 * old_XYZvox
      %
      %  where: M1 = M;   M1(:,4) = M(:,4) - sum(M(:,1:3),2);
      %  and:             M(:,4) == [T; 1] == sum(M1,2)
      %
      %  Therefore:   old_XYZvox = old_M1 \ new_M1 * new_XYZvox;
      %
      %  Since we are traverse Z axis, and   new_XYZvox   is replaced
      %  by   mask_Z * new_XYZvox, the above formula can be rewritten
      %  as:    old_XYZvox = old_M1 \ new_M1 * mask_Z * new_XYZvox;
      %
      %  i.e. we find the mapping from new_XYZvox to old_XYZvox:
      %  M = old_M1 \ new_M1 * mask_Z;
      %
      %  First, compute modified old_M1 & new_M1
      %
      old_M1 = old_M;   old_M1(:,4) = old_M(:,4) - sum(old_M(:,1:3),2);
      new_M1 = new_M;   new_M1(:,4) = new_M(:,4) - sum(new_M(:,1:3),2);

      %  Then, apply unit increment of mask_Z(3,4) to simulate the
      %  cursor movement
      %
      mask_Z(3,4) = z;

      %  Here is the mapping from new_XYZvox to old_XYZvox
      %
      M = old_M1 \ new_M1 * mask_Z;

      switch method
      case 1
         new_img(:,:,z) = trilinear(old_img, new_dim, old_dim, M, bg);
      case 2
         new_img(:,:,z) = nearest_neighbor(old_img, new_dim, old_dim, M, bg);
      case 3
         new_img(:,:,z) = bresenham(old_img, new_dim, old_dim, M, bg);
      end

   end;			% for z

   if ndims(old_img) == 2
      new_M(3,:) = [];
      new_M(:,3) = [];
   end

   return;					% affine

%--------------------------------------------------------------------
function img_slice = trilinear(img, dim1, dim2, M, bg)

   img_slice = zeros(dim1(1:2));
   TINY = 5e-2;					% tolerance

   %  Dimension of transformed 3D volume
   %
   xdim1 = dim1(1);
   ydim1 = dim1(2);

   %  Dimension of original 3D volume
   %
   xdim2 = dim2(1);
   ydim2 = dim2(2);
   zdim2 = dim2(3);

   %  initialize new_Y accumulation
   %
   Y2X = 0;
   Y2Y = 0;
   Y2Z = 0;

   for y = 1:ydim1

      %  increment of new_Y accumulation
      %
      Y2X = Y2X + M(1,2);		% new_Y to old_X
      Y2Y = Y2Y + M(2,2);		% new_Y to old_Y
      Y2Z = Y2Z + M(3,2);		% new_Y to old_Z

      %  backproject new_Y accumulation and translation to old_XYZ
      %
      old_X = Y2X + M(1,4);
      old_Y = Y2Y + M(2,4);
      old_Z = Y2Z + M(3,4);

      for x = 1:xdim1

         %  accumulate the increment of new_X, and apply it
         %  to the backprojected old_XYZ
         %
         old_X = M(1,1) + old_X  ;
         old_Y = M(2,1) + old_Y  ;
         old_Z = M(3,1) + old_Z  ;

         %  within boundary of original image
         %
         if (	old_X > 1-TINY & old_X < xdim2+TINY & ...
		old_Y > 1-TINY & old_Y < ydim2+TINY & ...
		old_Z > 1-TINY & old_Z < zdim2+TINY	)

            %  Calculate distance of old_XYZ to its neighbors for
            %  weighted intensity average
            %
            dx = old_X - floor(old_X);
            dy = old_Y - floor(old_Y);
            dz = old_Z - floor(old_Z);

            x000 = floor(old_X);
            x100 = x000 + 1;

            if floor(old_X) < 1
               x000 = 1;
               x100 = x000;
            elseif floor(old_X) > xdim2-1
               x000 = xdim2;
               x100 = x000;
            end

            x010 = x000;
            x001 = x000;
            x011 = x000;

            x110 = x100;
            x101 = x100;
            x111 = x100;

            y000 = floor(old_Y);
            y010 = y000 + 1;

            if floor(old_Y) < 1
               y000 = 1;
               y100 = y000;
            elseif floor(old_Y) > ydim2-1
               y000 = ydim2;
               y010 = y000;
            end

            y100 = y000;
            y001 = y000;
            y101 = y000;

            y110 = y010;
            y011 = y010;
            y111 = y010;

            z000 = floor(old_Z);
            z001 = z000 + 1;

            if floor(old_Z) < 1
               z000 = 1;
               z001 = z000;
            elseif floor(old_Z) > zdim2-1
               z000 = zdim2;
               z001 = z000;
            end

            z100 = z000;
            z010 = z000;
            z110 = z000;

            z101 = z001;
            z011 = z001;
            z111 = z001;

            x010 = x000;
            x001 = x000;
            x011 = x000;

            x110 = x100;
            x101 = x100;
            x111 = x100;

            v000 = double(img(x000, y000, z000));
            v010 = double(img(x010, y010, z010));
            v001 = double(img(x001, y001, z001));
            v011 = double(img(x011, y011, z011));

            v100 = double(img(x100, y100, z100));
            v110 = double(img(x110, y110, z110));
            v101 = double(img(x101, y101, z101));
            v111 = double(img(x111, y111, z111));

            img_slice(x,y) = v000*(1-dx)*(1-dy)*(1-dz) + ...
               v010*(1-dx)*dy*(1-dz) + ...
               v001*(1-dx)*(1-dy)*dz + ...
               v011*(1-dx)*dy*dz + ...
               v100*dx*(1-dy)*(1-dz) + ...
               v110*dx*dy*(1-dz) + ...
               v101*dx*(1-dy)*dz + ...
               v111*dx*dy*dz;

         else
            img_slice(x,y) = bg;

         end	% if boundary

      end	% for x
   end		% for y

   return;					% trilinear

%--------------------------------------------------------------------
function img_slice = nearest_neighbor(img, dim1, dim2, M, bg)

   img_slice = zeros(dim1(1:2));

   %  Dimension of transformed 3D volume
   %
   xdim1 = dim1(1);
   ydim1 = dim1(2);

   %  Dimension of original 3D volume
   %
   xdim2 = dim2(1);
   ydim2 = dim2(2);
   zdim2 = dim2(3);

   %  initialize new_Y accumulation
   %
   Y2X = 0;
   Y2Y = 0;
   Y2Z = 0;

   for y = 1:ydim1

      %  increment of new_Y accumulation
      %
      Y2X = Y2X + M(1,2);		% new_Y to old_X
      Y2Y = Y2Y + M(2,2);		% new_Y to old_Y
      Y2Z = Y2Z + M(3,2);		% new_Y to old_Z

      %  backproject new_Y accumulation and translation to old_XYZ
      %
      old_X = Y2X + M(1,4);
      old_Y = Y2Y + M(2,4);
      old_Z = Y2Z + M(3,4);

      for x = 1:xdim1

         %  accumulate the increment of new_X and apply it
         %  to the backprojected old_XYZ
         %
         old_X = M(1,1) + old_X  ;
         old_Y = M(2,1) + old_Y  ;
         old_Z = M(3,1) + old_Z  ;

         xi = round(old_X);
         yi = round(old_Y);
         zi = round(old_Z);

         %  within boundary of original image
         %
         if (	xi >= 1 & xi <= xdim2 & ...
		yi >= 1 & yi <= ydim2 & ...
		zi >= 1 & zi <= zdim2	)

            img_slice(x,y) = img(xi,yi,zi);

         else
            img_slice(x,y) = bg;

         end	% if boundary

      end	% for x
   end		% for y

   return;					% nearest_neighbor

%--------------------------------------------------------------------
function img_slice = bresenham(img, dim1, dim2, M, bg)

   img_slice = zeros(dim1(1:2));

   %  Dimension of transformed 3D volume
   %
   xdim1 = dim1(1);
   ydim1 = dim1(2);

   %  Dimension of original 3D volume
   %
   xdim2 = dim2(1);
   ydim2 = dim2(2);
   zdim2 = dim2(3);

   for y = 1:ydim1

      start_old_XYZ = round(M*[0     y 0 1]');
      end_old_XYZ   = round(M*[xdim1 y 0 1]');

      [X Y Z] = bresenham_line3d(start_old_XYZ, end_old_XYZ);

      %  line error correction
      %
%      del = end_old_XYZ - start_old_XYZ;
 %     del_dom = max(del);
  %    idx_dom = find(del==del_dom);
   %   idx_dom = idx_dom(1);
    %  idx_other = [1 2 3];
     % idx_other(idx_dom) = [];
      %del_x1 = del(idx_other(1));
%      del_x2 = del(idx_other(2));
 %     line_slope = sqrt((del_x1/del_dom)^2 + (del_x2/del_dom)^2 + 1);
  %    line_error = line_slope - 1;
% line error correction removed because it is too slow

      for x = 1:xdim1

         %  rescale ratio
         %
         i = round(x * length(X) / xdim1);

         if i < 1
            i = 1;
         elseif i > length(X)
            i = length(X);
         end

         xi = X(i);
         yi = Y(i);
         zi = Z(i);

         %  within boundary of the old XYZ space
         %
         if (	xi >= 1 & xi <= xdim2 & ...
		yi >= 1 & yi <= ydim2 & ...
		zi >= 1 & zi <= zdim2	)

            img_slice(x,y) = img(xi,yi,zi);

%            if line_error > 1
 %              x = x + 1;

%               if x <= xdim1
 %                 img_slice(x,y) = img(xi,yi,zi);
  %                line_error = line_slope - 1;
   %            end
    %        end		% if line_error
% line error correction removed because it is too slow

         else
            img_slice(x,y) = bg;

         end	% if boundary

      end	% for x
   end		% for y

   return;					% bresenham

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function crossfeed(varargin)

%CROSSFEED Checks for crossfeeding and duplicate data channels
% CROSSFEED Reads a short snippet of data from an entire network and cross
% correlates each trace against all others (without time shifting). Events
% are then clustered and a directory of png images is written out showing
% raw and filtered versions of the waveforms in each cluster. These images
% can be used to check for data problems across the network. The image
% directory is named FIG_yyyymmdd_HHMMSS according to the start time.
%
% CROSSFEED(STARTTIME) specify a start time in Matlab date format
%
% CROSSFEED(STARTTIME,DURATION) specify a trace duration in seconds.
%
% CROSSFEED(STARTTIME,DURATION,THRESHOLD) specify a cross correlation
% threshold for clustering.
%
% CROSSFEED(STARTTIME,DURATION,THRESHOLD,NETWORK) specify a network to
% examine.
%
% Using CROSSFEED with no arguments is equivalent to:
%    CROSSFEED( floor(now)-0.5 , 30 , 0.7 , 'AV' )
%
% Note that "now" is returned in local time. So the default data examined
% is 8 to 9 hours old, assuming AKST/AKDT.
%
% -- CAVEATS --
% This function is hardwired to run in the UAF Seis Lab assuming one day
% archive databases that contain waveforms and station-related tables (to
% get network affiliations)

% READ ARGUMENTS AND SET DEFAULT VALUES
if numel(varargin)>=1
    startTime = varargin{1};
else
     startTime = floor(now)-0.5;
end

if numel(varargin)>=2
    duration = varargin{2}/86400;
else
     duration = 30/86400;
end

if numel(varargin)>=3
    threshold = varargin{3};
else
     threshold = 0.7;
end

if numel(varargin)==4
    networkName = varargin{4};
else
     networkName = 'AV';
end

% GET CHANNEL LIST FROM DATABASE
dbWfdisc = ['/aerun/sum/db/archive/archive_' datestr(startTime,'yyyy') '/archive_' datestr(startTime,'yyyy') '_' datestr(startTime,'mm') '_' datestr(startTime,'dd')];
disp(['Reading STA_CHANs for network ' networkName ' from database: ']);
disp([ '   ' dbWfdisc ' ...']);

db = dbopen(dbWfdisc,'r');
db = dblookup(db,'','wfdisc','','');
db1 = dblookup(db,'','affiliation','','');
db = dbjoin(db,db1);
db = dbsubset(db,['net=="' networkName '"']);
[wfdisc.sta wfdisc.chan wfdisc.net] = dbgetv(db,'sta','chan','net');
dbclose(db);

% GET LIST OF UNIQUE STA CHANs
staChanNetList = strcat(wfdisc.sta,'#',wfdisc.chan,'#',wfdisc.net);
staChanNetList = unique(staChanNetList);
for n = 1:numel(staChanNetList)
    tmp = regexp(staChanNetList(n), '#', 'split');
    scnl(n) = scnlobject(tmp{1}{1},tmp{1}{2},tmp{1}{3},'');
end

% READ IN WAVEFORMS
disp(['Loading ' num2str(numel(scnl)) ' waveforms for this time period ...']);
disp(['   start time: ' datestr(startTime)]);
disp(['   duration:   ' num2str(duration*86400) ' seconds']);
ds = datasource('antelope',dbWfdisc);
w = waveform(ds,scnl,startTime,startTime+duration);
w = demean(w);
w = fillgaps(w,0);
w = detrend(w);

% PREP WAVEFORMS AND CROSS CORRELATE
c = correlation(w);
c = detrend(c);
c = demean(c);
c = taper(c);
c = butter(c,[2 12]);
%c = xcorr(c);

% DO ZERO LAG XCORR
D = get(c,'DATA');
wcoeff = 1./sqrt(sum(D.*D));
%wcoeff(find(wcoeff==Inf)) = 0;
corr = zeros(size(D,2));
lag = zeros(size(D,2));
for col = 1:size(D,2)
    for row = 1:size(D,2)
        corr(col,row) = sum(D(:,col).*D(:,row)) .* wcoeff(col) .* wcoeff(row);
    end
end
f = find(corr<0);
corr(f) = 0;
f = find(isnan(corr));
corr(f) = 0;
c = set(c,'corr',corr);
c = set(c,'lag',lag);
c = linkage(c);
c = cluster(c,threshold);

% MAKE FIGURE DIRECTORY
dirName = [networkName '_' datestr(startTime,'yyyy') datestr(startTime,'mm') datestr(startTime,'dd') '_' datestr(startTime,'HH') datestr(startTime,'MM') datestr(startTime,'SS')];
if ~exist(dirName,'dir')
    disp(['Creating new directory: ' dirName]);
    mkdir(dirName);
end

% PLOT TRACE CLUSTERS
index = find(c,'BIG',2);
clust = get(c,'CLUST');
disp('--> dbpick selection strings: ');
for n = 1:max(clust(index))
    f = find(clust==n);
    figure('Color','w','Position',[50 50 800 1100]);
    set(gcf,'DefaultAxesFontSize',14)

    % RAW TRACE FIGURE

    wOrig = set(w(f),'UNITS','normalized scale');
    wOrigCrop = extract(wOrig,'TIME',startTime+5/86400,startTime+25/86400);
    wOrigCrop = demean(detrend(wOrigCrop));
    for i = 1:numel(wOrigCrop)
       wOrigCrop(i) = 0.25 * wOrigCrop(i) ./ std(wOrigCrop(i)) + i;
    end
    wOrigCrop = set(wOrigCrop,'UNITS','normalized scale');

    subplot(2,1,1)
    plot(wOrigCrop);
    legend(wOrigCrop);
    xlim([0 20]);
    ylim([0 numel(wOrig)+1]);
    text(19.3,0.05 *(numel(wOrig)+1),['correlation threshold: ' num2str(threshold)],'HorizontalAlignment','Right','FontSize',12);
    text(19.3,0.10 *(numel(wOrig)+1),['trace duration: ' num2str(duration*86400) 's'],'HorizontalAlignment','Right','FontSize',12);
    text(19.3,0.15 *(numel(wOrig)+1),'Raw traces','HorizontalAlignment','Right','FontSize',12,'FontWeight','Bold');

    % PLOT FILTERED CROPPED TRACE
    subplot(2,1,2)
    c2 = subset(c,f);
    wFilt = waveform(c2);
    wFilt = set(wFilt,'UNITS','normalized scale');
    wFiltCrop = extract(wFilt,'TIME',startTime+10/86400,startTime+13/86400);
    for i = 1:numel(wOrig)
       wFiltCrop(i) = 0.25 * wFiltCrop(i) ./ std(wFiltCrop(i));
    end
    plot(wFiltCrop,'Linewidth',1);
    legend(wFiltCrop);
    xlim([0 3]);
    ylim([-0.8 0.8]);
    text(2.9,-0.75,'Filtered on 2-12 Hz','HorizontalAlignment','Right','FontSize',12,'FontWeight','Bold');

    % WRITE OUT FIGURES
    subSta = get(wFilt,'STATION');
    subChan = get(wFilt,'CHANNEL');
    scList = [];
    fileName = [dirName '/overlays_' subSta{1} '_' subChan{1} '-' subSta{2} '_' subChan{2} '.png'];
    set(gcf, 'paperorientation', 'portrait');
    set(gcf, 'paperposition', [.25 .25 8 11] );
    print('-dpng',fileName);
    %
    close all

    % WRITE DBPICK TEXT
    scList = [];
    for i = 1:numel(subSta)
        scList = [scList subSta{i} ':' subChan{i} ','];
    end
    scList = scList(1:end-1);
    disp(['sc ' scList]);
end
disp(['ts ' datestr(startTime,'yyyy/mm/dd HH:MM:SS')]);


text:
function data = plotData(f, g, h)
%PLOTDATA   Useful data values for plotting a CHEBTECH object.
%   DATA = PLOTDATA(F) returns a struct containing data that can be used for
%   plotting F. The struct DATA contains the following fields:
%
%       xLine: x-coordinates of for plotting smooth curves.
%       yLine: Function values of F at the coordinates stored in xLine.
%       xPoints: x-coordinates of the Chebyshev points used to represent F.
%       yPoints: Function value at the Chebyshev points used to represent F.
%
%   DATA.xLine and DATA.yLine are used for plotting smooth curves (usually
%   passed to PLOT() with the '-' option).
%
%   DATA.xPoints and DATA.yPoints contain the (x, F(x)) data at the Chebyshev
%   grid used to represent F, and are used for plots with marks (e.g.
%   PLOT(F,'-o').
%
%   DATA = PLOTDATA(F, G) is similar but for plot calls of the form PLOT(F, G),
%   where both F and G are CHEBTECH objects.
%
%   DATA = PLOTDATA(F, G, H) is for plots of the form PLOT3(F, G, H). In this
%   instance, DATA also contains fields zLine and zPoints for the data
%   corresponding to H.
%
% See also PLOT.

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

if ( nargin == 1 )
    g = [];
end
if ( nargin < 3 )
    h = [];
end

% Get the number of points: (Oversample the wavelength)
len = max([length(f), length(g), length(h)]);
npts = min(max(501, round(4*pi*len)), chebtech.techPref().maxLength);

% Initialise the output structure:
data = struct('xLine', [], 'yLine', [], 'xPoints', [], 'yPoints', [], ...
    'xLim', [], 'yLim', [], 'defaultXLim', 1, 'defaultYLim', 1);

if ( isempty(g) )
    % PLOT(F):

    % Values on oversampled Chebyshev grid (faster than evaluating a uniform grid!).
    data.xLine = f.chebpts(npts);
    data.yLine = get(prolong(f, npts), 'values');

    % Values on the Cheyshev grid tied to the CHEBTECH F:
    data.xPoints = f.points();
    data.yPoints = f.coeffs2vals(f.coeffs);

    % yLim:
    data.yLim = [min(data.yLine(:)) max(data.yLine(:))];

    % xLim:
    data.xLim = [-1 1];

elseif ( isa(g, 'chebtech') )
    % PLOT(F, G)

    % Also return the grid points used.
    % Grid data for f:
    data.fGrid.xLine = f.chebpts(npts);
    % Use the maximum of the lenghts of f, g and h to match the number of
    % values returned:
    data.fGrid.xPoints = f.chebpts(len);

    % Grid data for g:
    data.gGrid.xLine = g.chebpts(npts);
    % Use the maximum of the lenghts of f, g and h to match the number of
    % values returned:
    data.gGrid.xPoints = g.chebpts(len);

    % Values on oversampled Chebyshev grid (faster than evaluating a uniform grid!).
    data.xLine = get(prolong(f, npts), 'values');
    data.yLine = get(prolong(g, npts), 'values');

    % Values on the largest Cheyshev grid tied to the CHEBTECH objects F and G:
    data.xPoints = get(prolong(f, len), 'values');
    data.yPoints = get(prolong(g, len), 'values');

    % xLim:
    xdata = [get(f, 'lval'); data.xLine; get(f, 'rval')];
    data.xLim = [min(xdata(:)) max(xdata(:))];

    % yLim:
    ydata = [get(g, 'lval'); data.yLine; get(g, 'rval')];
    data.yLim = [min(ydata(:)) max(ydata(:))];

    if ( isa(h, 'chebtech') )
        % PLOT3(F, G, H)

        % Grid data for h:
        data.hGrid.xLine = h.chebpts(npts);
        % Use the maximum of the lenghts of f, g and h to match the number of
        % values returned:
        data.hGrid.xPoints = h.chebpts(len);

         % Values on oversampled Chebyshev grid:
        data.zLine = get(prolong(h, npts), 'values');
        data.zPoints = get(prolong(h, len), 'values');

    end

else

    error('CHEBFUN:CHEBTECH:plotData:dataType', 'Invalid data types.');

end

end


text:
function model_3d = train_nrsfm(cls,train_set)
    globals;
    params = get_params();
    md = isnan(train_set.points(1:end/2,:));
    points = train_set.points;
    points_bak = points;
    points(isnan(points))=0;
    try
        seg.poly_x = train_set.poly_x;
        seg.poly_y = train_set.poly_y;
    catch
        seg = train_set.mask;
    end

    % Setup variables for converting to PASCAL 3D frame
    load(fullfile(datadir,'partNames',cls));
    load(fullfile(datadir,'voc_kp_metadata'));
    cInd = find(ismember(metadata.categories,cls));
    [~,I] = sort(metadata.kp_names{cInd});
    invI(I) = 1:numel(I);
    rightCoordSys = metadata.right_coordinate_sys{cInd};
    rightCoordSys(1:6) = invI(rightCoordSys(:));
    names = partNames(rightCoordSys);

    % Do non-rigid SFM
    [P3, S_bar, V, RO, Tr, Z, sigma_sq, ~, ~, ~, ~,c]...
      = em_sfm(points, md, params.nrsfm.nBasis, seg, train_set.imsize, train_set.bbox,...
      params.nrsfm.tol, params.nrsfm.max_em_iter,rightCoordSys,train_set.rotP3d);

    model_3d.class          = cls;
    model_3d.S_bar          = S_bar;
    model_3d.defBasis       = V;
    model_3d.part_names     = train_set.labels;
    model_3d.sigma_sq       = sigma_sq;
    model_3d.points3        = P3;
    model_3d.rots           = RO;
    model_3d.trs            = Tr;
    model_3d.def_weights    = Z;
    model_3d.voc_image_id   = train_set.voc_image_id;
    if(isfield(train_set,'voc_rec_id'))
        model_3d.voc_rec_id     = train_set.voc_rec_id;
    end
    model_3d.bbox           = train_set.bbox;
    model_3d.points2        = points_bak;
    model_3d.params         = params;
    model_3d.c              = c;
    model_3d.seg            = seg;
    if(isfield(train_set,'flip'))
        model_3d.flip = train_set.flip;
    end
    if(isfield(train_set,'subtype'))
        model_3d.subtype = train_set.subtype;
    end
    if(isfield(train_set,'rotP3d'))
        model_3d.rotP3d = train_set.rotP3d;
    end

end

repo: CategoryShapes


text:
function [ f ] = DataDensityPlot( x, y, levels )
%DATADENSITYPLOT Plot the data density
%   Makes a contour map of data density
%   x, y - data x and y coordinates
%   levels - number of contours to show
%
% By Malcolm Mclean
%
    map = dataDensity(x, y, 256, 256);
    map = map - min(min(map));
    map = floor(map ./ max(max(map)) * (levels-1));
    f = figure();

    image(map);
    colormap(jet(levels));
    set(gca, 'XTick', [1 256]);
    set(gca, 'XTickLabel', [min(x) max(x)]);
    set(gca, 'YTick', [1 256]);
    set(gca, 'YTickLabel', [min(y) max(y)]);
    uiwait;
end


text:
function [xOpt,fval] = psorunhybridfcn(fitnessfcn,xOpt,...
    Aineq,bineq,Aeq,beq,LB,UB,nonlcon,options)
% Calls the hybrid function defined in options.HybridFcn, from PSO.

if iscell(options.HybridFcn)
    HybridFcn = options.HybridFcn{1} ;
    hybridOptions = options.HybridFcn{2} ;
else
    HybridFcn = options.HybridFcn ;
    hybridOptions = optimset('Display',options.Display,...
        'LargeScale','off') ;
    if isfield(hybridOptions,'Algorithm')
        hybridOptions.Algorithm = 'active-set' ;
    end
end

if options.Verbosity > 0
    fprintf('\nBest point before hybrid function: %s',...
        mat2str(xOpt,5))
    fprintf('\n\nTurning over to hybrid function %s...\n\n',...
        func2str(HybridFcn))
end

if exist(func2str(HybridFcn),'file') ~= 2
    warning('pso:hybridfcn:nofile',...
        'Hybrid function %s cannot be found. Check toolboxes.',...
        func2str(HybridFcn))
    fval = fitnessfcn(xOpt) ;
    return
end

% Check for constraints
if strcmp(func2str(HybridFcn),func2str(@fmincon)) && ...
        all([isempty([Aineq,bineq]), isempty([Aeq,beq]), ...
        isempty([LB;UB]),isempty(nonlcon)])
    if options.Verbosity > 0
        msg = sprintf('Warning: %s does not accept problems without',...
            func2str(HybridFcn)) ;
        fprintf('%s constraints. Switching to fminunc.\n\n',msg)
    end % if options.Verbosity
    HybridFcn = @fminunc ;
elseif (strcmp(func2str(HybridFcn),func2str(@fminunc)) || ...
        strcmp(func2str(HybridFcn),func2str(@fminsearch))) && ...
        ~all([isempty([Aineq,bineq]), isempty([Aeq,beq]), ...
        isempty([LB;UB]),isempty(nonlcon)])
    if options.Verbosity > 0
        msg = sprintf('Warning: %s does not accept problems with',...
            func2str(HybridFcn)) ;
        fprintf(0,'%s constraints. Switching to fmincon.\n\n',msg)
    end % if options.Verbosity
    HybridFcn = @fmincon ;
end

if strcmp(func2str(HybridFcn),func2str(@fmincon)) || ...
        strcmp(func2str(HybridFcn),func2str(@patternsearch))
    [xOpt,fval] = HybridFcn(fitnessfcn,xOpt,Aineq,bineq,...
        Aeq,beq,LB,UB,nonlcon,hybridOptions) ;
elseif strcmp(func2str(HybridFcn),func2str(@fminunc)) || ...
        strcmp(func2str(HybridFcn),func2str(@fminsearch))
    [xOpt,fval] = HybridFcn(fitnessfcn,xOpt,hybridOptions) ;
else
    warning('pso:hybridfcn:unrecognized',...
        'Unrecognized hybrid function. Ignoring for it now.')
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function y=absx2(x)

% y=absx2(x)
%
% Faster implementation of y=abs(x).^2.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% 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/>.

y=real(x).*real(x)+imag(x).*imag(x);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function rtk = detslip_gf_L1L5(rtk,obsr,obsb,nav)

sat=obsr.sat;

g1=gfobs_L1L5(obsr,obsb,nav.lam(sat,:));
if rtk.opt.nf<=2||g1==0,return;end

g0=rtk.sat(sat).gf2;
rtk.sat(sat).gf2=g1;

if g0~=0&&abs(g1-g0)>rtk.opt.csthres(1)
    rtk.sat(sat).slip(1)=bitor(rtk.sat(sat).slip(1),1);
    rtk.sat(sat).slip(3)=bitor(rtk.sat(sat).slip(3),1);
end

return

repo: GINav


text:
function [masked_image] = separate(working_image,background_dominant)
%Seperate function: Using a sum-of-gaussians model (2 gaussians) for the color at each pixel results in a simple way of seperating the background.
%Written by Alexander Farley, April 19 2011
%alexander.farley at utoronto.ca
%Updated October 8 2011: Added a flag for pictures where the foreground is dominant
%
%
%working_image: the input image to be counted
%background_dominant: set to 1 if the background occupies more than 50% of
%the image area, otherwise set to 0
%masked_image: binary background/cell separated image

[image_height image_width num_colors] = size(working_image);

%Seperate into color columns
Acolumn = working_image(:,:,1);
Bcolumn = working_image(:,:,2);
Ccolumn = working_image(:,:,3);

%Show original color data in 3D space
% figure
 Areduced = downsample(Acolumn(:),1000);
 Breduced = downsample(Bcolumn(:),1000);
 Creduced = downsample(Ccolumn(:),1000);
% scatter3(Areduced, Breduced, Creduced);

%Combine columns
sample_data = [Acolumn(:)'; Bcolumn(:)'; Ccolumn(:)';];
sample_data_reduced = [Areduced(:)'; Breduced(:)'; Creduced(:)'];

[W,M,V,L] = EM_GM(sample_data_reduced', 2);

%Calculate pdf using Gaussian parameters estimated using EM
membership = zeros(length(Acolumn(:)), 2);
for j=1:length(Acolumn(:))
    membership(j,1) = W(1)*mvnpdf(sample_data(:,j), M(:,1), V(:,:,1));
    membership(j,2) = W(2)*mvnpdf(sample_data(:,j), M(:,2), V(:,:,2));
end

is_bg = membership(:,2) > membership(:,1);
bg_mask = reshape(is_bg, image_height,image_width);
if background_dominant
if W(1) > W(2)
fg_mask = bg_mask;
else
    fg_mask = ~bg_mask;
end
else
if W(1) > W(2)
fg_mask = ~bg_mask;
else
    fg_mask = bg_mask;
end
end

masked_image = working_image;
masked_image(:,:,1) = fg_mask.*masked_image(:,:,1);
masked_image(:,:,2) = fg_mask.*masked_image(:,:,2);
masked_image(:,:,3) = fg_mask.*masked_image(:,:,3);

masked_image = masked_image > 1;