text:
function  c = fchcode(b, conn, dir)
%FCHCODE Computes the Freeman chain code of a boundary.
%   C = FCHCODE(B) computes the 8-connected Freeman chain code of a
%   set of 2-D coordinate pairs contained in B, an np-by-2 array. C
%   is a structure with the following fields:
%
%     c.fcc    = Freeman chain code (1-by-np)
%     c.diff   = First difference of code c.fcc (1-by-np)
%     c.mm     = Integer of minimum magnitude from c.fcc (1-by-np)
%     c.diffmm = First difference of code c.mm (1-by-np)
%     c.x0y0   = Coordinates where the code starts (1-by-2)
%
%   C = FCHCODE(B, CONN) produces the same outputs as above, but
%   with the code connectivity specified in CONN. CONN can be 8 for
%   an 8-connected chain code, or CONN can be 4 for a 4-connected
%   chain code. Specifying CONN=4 is valid only if the input
%   sequence, B, contains transitions with values 0, 2, 4, and 6,
%   exclusively.
%
%   C = FHCODE(B, CONN, DIR) produces the same outputs as above, but,
%   in addition, the desired code direction is specified. Values for
%   DIR can be:
%
%     'same'      Same as the order of the sequence of points in b.
%                 This is the default.
%
%     'reverse'   Outputs the code in the direction opposite to the
%                 direction of the points in B.  The starting point
%                 for each DIR is the same.
%
%   The elements of B are assumed to correspond to a 1-pixel-thick,
%   fully-connected, closed boundary. B cannot contain duplicate
%   coordinate pairs, except in the first and last positions, which
%   is a common feature of boundary tracing programs.
%
%   FREEMAN CHAIN CODE REPRESENTATION
%   The table on the left shows the 8-connected Freeman chain codes
%   corresponding to allowed deltax, deltay pairs. An 8-chain is
%   converted to a 4-chain if (1) if conn = 4; and (2) only
%   transitions 0, 2, 4, and 6 occur in the 8-code.  Note that
%   dividing 0, 2, 4, and 6 by 2 produce the 4-code.
%
%       -----------------------  ----------------
%       deltax | deltay | 8-code  corresp 4-code
%       -----------------------  ----------------
%         0        1       0            0
%        -1        1       1
%        -1        0       2            1
%        -1       -1       3
%         0       -1       4            2
%         1       -1       5
%         1        0       6            3
%         1        1       7
%       -----------------------  ----------------
%
%   The formula z = 4*(deltax + 2) + (deltay + 2) gives the following
%   sequence corresponding to rows 1-8 in the preceding table: z =
%   11,7,6,5,9,13,14,15. These values can be used as indices into the
%   table, improving the speed of computing the chain code. The
%   preceding formula is not unique, but it is based on the smallest
%   integers (4 and 2) that are powers of 2.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2003/11/21 14:34:49 $

% Preliminaries.
if nargin == 1
   dir = 'same';
   conn = 8;
elseif nargin == 2
   dir = 'same';
elseif nargin == 3
   % Nothing to do here.
else
   error('Incorrect number of inputs.')
end
[np, nc] = size(b);
if np < nc
   error('B must be of size np-by-2.');
end

% Some boundary tracing programs, such as boundaries.m, output a
% sequence in which the coordinates of the first and last points are
% the same. If this is the case, eliminate the last point.
if isequal(b(1, :), b(np, :))
   np = np - 1;
   b = b(1:np, :);
end

% Build the code table using the single indices from the formula
% for z given above:
C(11)=0; C(7)=1; C(6)=2; C(5)=3; C(9)=4;
C(13)=5; C(14)=6; C(15)=7;

% End of Preliminaries.

% Begin processing.
x0 = b(1, 1);
y0 = b(1, 2);
c.x0y0 = [x0, y0];

% Make sure the coordinates are organized sequentially:
% Get the deltax and deltay between successive points in b. The
% last row of a is the first row of b.
a = circshift(b, [-1, 0]);

% DEL = a - b is an nr-by-2 matrix in which the rows contain the
% deltax and deltay between successive points in b. The two
% components in the kth row of matrix DEL are deltax and deltay
% between point (xk, yk) and (xk+1, yk+1).  The last row of DEL
% contains the deltax and deltay between (xnr, ynr) and (x1, y1),
% (i.e., between the last and first points in b).
DEL = a - b;

% If the abs value of either (or both) components of a pair
% (deltax, deltay) is greater than 1, then by definition the curve
% is broken (or the points are out of order), and the program
% terminates.
if any(abs(DEL(:, 1)) > 1) | any(abs(DEL(:, 2)) > 1);
   error('The input curve is broken or points are out of order.')
end

% Create a single index vector using the formula described above.
z = 4*(DEL(:, 1) + 2) + (DEL(:, 2) + 2);

% Use the index to map into the table. The following are
% the Freeman 8-chain codes, organized in a 1-by-np array.
fcc = C(z);

% Check if direction of code sequence needs to be reversed.
if strcmp(dir, 'reverse')
   fcc = coderev(fcc); % See below for function coderev.
end

% If 4-connectivity is specified, check that all components
% of fcc are 0, 2, 4, or 6.
if conn == 4
   val = find(fcc == 1 | fcc == 3 | fcc == 5 | fcc ==7 );
   if isempty(val)
      fcc = fcc./2;
   else
      warning('The specified 4-connected code cannot be satisfied.')
   end
end

% Freeman chain code for structure output.
c.fcc = fcc;

% Obtain the first difference of fcc.
c.diff = codediff(fcc,conn); % See below for function codediff.

% Obtain code of the integer of minimum magnitude.
c.mm = minmag(fcc); % See below for function minmag.

% Obtain the first difference of fcc
c.diffmm = codediff(c.mm, conn);

%-------------------------------------------------------------------%
function cr = coderev(fcc)
%   Traverses the sequence of 8-connected Freeman chain code fcc in
%   the opposite direction, changing the values of each code
%   segment. The starting point is not changed. fcc is a 1-by-np
%   array.

% Flip the array left to right.  This redefines the starting point
% as the last point and reverses the order of "travel" through the
% code.
cr = fliplr(fcc);

% Next, obtain the new code values by traversing the code in the
% opposite direction. (0 becomes 4, 1 becomes 5, ... , 5 becomes 1,
% 6 becomes 2, and 7 becomes 3).
ind1 = find(0 <= cr & cr <= 3);
ind2 = find(4 <= cr & cr <= 7);
cr(ind1) = cr(ind1) + 4;
cr(ind2) = cr(ind2) - 4;

%-------------------------------------------------------------------%
function z = minmag(c)
%MINMAG Finds the integer of minimum magnitude in a chain code.
%   Z = MINMAG(C) finds the integer of minimum magnitude in a given
%   4- or 8-connected Freeman chain code, C. The code is assumed to
%   be a 1-by-np array.

% The integer of minimum magnitude starts with min(c), but there
% may be more than one such value. Find them all,
I = find(c == min(c));
% and shift each one left so that it starts with min(c).
J = 0;
A = zeros(length(I), length(c));
for k = I;
   J = J + 1;
   A(J, :) = circshift(c,[0 -(k-1)]);
end

% Matrix A contains all the possible candidates for the integer of
% minimum magnitude. Starting with the 2nd column, succesively find
% the minima in each column of A. The number of candidates decreases
% as the seach moves to the right on A.  This is reflected in the
% elements of J.  When length(J)=1, one candidate remains.  This is
% the integer of minimum magnitude.
[M, N] = size(A);
J = (1:M)';
for k = 2:N
   D(1:M, 1) = Inf;
   D(J, 1) = A(J, k);
   amin = min(A(J, k));
   J = find(D(:, 1) == amin);
   if length(J)==1
      z = A(J, :);
      return
   end
end

%-------------------------------------------------------------------%
function d = codediff(fcc, conn)
%CODEDIFF Computes the first difference of a chain code.
%   D = CODEDIFF(FCC) computes the first difference of code, FCC. The
%   code FCC is treated as a circular sequence, so the last element
%   of D is the difference between the last and first elements of
%   FCC.  The input code is a 1-by-np vector.
%
%   The first difference is found by counting the number of direction
%   changes (in a counter-clockwise direction) that separate two
%   adjacent elements of the code.

sr = circshift(fcc, [0, -1]); % Shift input left by 1 location.
delta = sr - fcc;
d = delta;
I = find(delta < 0);

type = conn;
switch type
case 4 % Code is 4-connected
   d(I) = d(I) + 4;
case 8 % Code is 8-connected
   d(I) = d(I) + 8;
end

author: HuangCongQing


text:
classdef LIRCMOP4 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP with large infeasible regions

%------------------------------- Reference --------------------------------
% Z. Fan, W. Li, X. Cai, H. Huang, Y. Fang, Y. You, J. Mo, C. Wei, and E.
% Goodman, An improved epsilon constraint-handling method in MOEA/D for
% CMOPs with large infeasible regions, Soft Computing, 2019, 23:
% 12491-12510.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenji Li

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            x_odd       = X(:,3:2:end);
            x_even      = X(:,2:2:end);
            len_odd     = size(x_odd,2);
            len_even    = size(x_even,2);
            g_1         = sum((x_odd - repmat(X(:,1),1,len_odd)).^2,2);
            g_2         = sum((x_even - repmat(X(:,1),1,len_even)).^2,2);
            PopObj(:,1) = X(:,1) + g_1;
            PopObj(:,2) = 1 - sqrt(X(:,1)) + g_2;
            PopCon(:,1) = (0.5 - g_1).*(0.51 - g_1);
            PopCon(:,2) = (0.5 - g_2).*(0.51 - g_2);
            PopCon(:,3) = 0.5 - sin(20 * pi * X(:,1));
            Population  = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
            R(sin(20*pi*R(:,1))<0.5,:) = [];
            R      = R + 0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R(:,1) = linspace(0,1,100)';
            R(:,2) = 1 - sqrt(R(:,1));
            R(sin(20*pi*R(:,1))<0.5,:) = nan;
            R      = R + 0.5;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [M1,MW,MW1] = perform_steerable_matching(M1,M,options)

% perform_steerable_matching - match multiscale histograms
%
% M1 = perform_steerable_matching(M1,M,options);
%
%   M1 is the image to synthesize.
%   M is the exemplar image.
%
%   This function match the histogram of the image and the histogram
%   of each sub-band of a steerable pyramid.
%
%   To do texture synthesis, one should apply several time this function.
%   You can do it by setting the value of options.niter_synthesis.
%   This leads to the synthesis as described in
%
%       Pyramid-Based Texture Analysis/Synthesis
%       D. Heeger, J. Bergen,
%       Siggraph 1995
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
niter_synthesis = getoptions(options, 'niter_synthesis', 1);
copy_lowpass = getoptions(options, 'copy_lowpass', 1);

if isfield(options, 'color_mode') && strcmp(options.color_mode, 'pca') && ~isfield(options, 'ColorP') && size(M,3)==3
    [tmp,options.ColorP] = change_color_mode(M,+1,options);
end
rgb_postmatching = getoptions(options, 'rgb_postmatching', 0);

if size(M,3)==3
    options.niter_synthesis = 1;
    if not(isfield(options, 'color_mode'))
        options.color_mode = 'hsv';
    end
    for iter=1:niter_synthesis
        % color images
        M  = change_color_mode(M, +1,options);
        M1 = change_color_mode(M1,+1,options);
        for i=1:size(M,3)
            M1(:,:,i) = perform_steerable_matching(M1(:,:,i),M(:,:,i), options);
        end
        M  = change_color_mode(M, -1,options);
        M1 = change_color_mode(M1,-1,options);
        if rgb_postmatching
        for i=1:size(M,3)
            M1(:,:,i) = perform_histogram_equalization(M1(:,:,i),M(:,:,i));
        end
        end
    end
    return;
end

if size(M,3)>1
    for i=1:size(M,3)
        [M1(:,:,i),MW,MW1] = perform_steerable_matching(M1(:,:,i),M(:,:,i),options);
    end
    return;
end

if not(isfield(options, 'nb_orientations'))
    options.nb_orientations = 4;
end

nb_orientations = getoptions(options, 'nb_orientations', 1, 1);
Jmin = getoptions(options, 'JminSynthesis', 4);

for iter=1:niter_synthesis
    % spatial equalization
    M1 = perform_histogram_equalization(M1,M);
    % forward transforms
    MW1 = perform_steerable_transform(M1, Jmin, options);
    if iscell(M)
        MW = M;
    else
        MW = perform_steerable_transform(M, Jmin, options);
    end
    % wavelet domain equalization
    MW1 = perform_histogram_equalization(MW1,MW);
    if copy_lowpass
        % copy verbatim the low pass coefficients
        MW1{end} = MW{end};
    end
    % backward transform
    M1 = perform_steerable_transform(MW1, Jmin, options);
    % spatial equalization
    M1 = perform_histogram_equalization(M1,M);
end

author: gpeyre

repo: matlab-toolboxes


text:
function p = nodes_brick27 ( )

%*****************************************************************************80
%
%% NODES_BRICK27 returns the natural coordinates of the BRICK27 element.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real P(3,27), the coordinates.
%
  p = [ ...
    -1.0, -1.0, -1.0; ...
    +1.0, -1.0, -1.0; ...
    +1.0, +1.0, -1.0; ...
    -1.0, +1.0, -1.0; ...
    -1.0, -1.0, +1.0; ...
    +1.0, -1.0, +1.0; ...
    +1.0, +1.0, +1.0; ...
    -1.0, +1.0, +1.0; ...
     0.0, -1.0, -1.0; ...
    +1.0,  0.0, -1.0; ...
     0.0, +1.0, -1.0; ...
    -1.0,  0.0, -1.0; ...
    -1.0, -1.0,  0.0; ...
    +1.0, -1.0,  0.0; ...
    +1.0, +1.0,  0.0; ...
    -1.0, +1.0,  0.0; ...
     0.0, -1.0, +1.0; ...
    +1.0,  0.0, +1.0; ...
     0.0, +1.0, +1.0; ...
    -1.0,  0.0, +1.0; ...
     0.0,  0.0, -1.0; ...
     0.0, -1.0,  0.0; ...
    +1.0,  0.0,  0.0; ...
     0.0, +1.0,  0.0; ...
    -1.0,  0.0,  0.0; ...
     0.0,  0.0, +1.0; ...
     0.0,  0.0,  0.0 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Computation of the S_Method
% Usage:
%
%       output = S_method( input1, L, wl, win, param_overlap, fftl);
%
%
% Parameters:
%
%    signal
%
%         Input one dimensional signal to be analysed. An analytic signal
%         is required for this function, however, if signal is real, a
%         default analytic transformer routine will be called from this
%         function before computing spectrogram.
%
%    L
%
%        Define the length of frequency window in the S-Method (2*L+1)
%
%    window_length
%
%         Length of choosen window.
%
%
%    window_type
%
%         One of  'rect', 'hann', 'hamm', 'bart'
%
%   overlap
%
%         Number of samples in common between two consecutives windows.
%         0 < overlap < window_length.
%
%    fft_length
%
%         Zero-padding at the fft stage of the analysis may be specified by
%         giving an fft_length larger than normal. If fft_length is not
%         specified then the next highest power of two above the
%         signal length will be used. If fft_length is not a power
%         of two, then the next highest power of two is used.
%
%
%    tfd
%
%        The computed time-frequency distribution. size(tfd) will
%        return [a/2+1, b], where a is the next largest power of two above
%        window_length, and b is ceil((length(signal)-overlap)/(window_length-overlap)) .
%
%
%
%  See Also: spec
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
function y = htmlcolors(x)

% HTMLCOLORS looks up the RGB value for a named color (string), or the name for a given RGB value
%
% Use as
%   rgb = htmlcolors(name)
% or
%   name = htmlcolors(rgb)
% or
%   list = htmlcolors
%
% See https://www.rapidtables.com/web/color/html-color-codes.html
% and https://www.color-hex.com/color-palettes/
%
% See also FT_COLORMAP, COLORMAP, COLORMAPEDITOR, BREWERMAP, MATPLOTLIB, CMOCEAN

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

mapping = {
  'lightsalmon',           [255,160,122]
  'salmon',                [250,128,114]
  'darksalmon',            [233,150,122]
  'lightcoral',            [240,128,128]
  'indianred',             [205,92,92]
  'crimson',               [220,20,60]
  'firebrick',             [178,34,34]
  'red',                   [255,0,0]
  'darkred',               [139,0,0]
  'coral',                 [255,127,80]
  'tomato',                [255,99,71]
  'orangered',             [255,69,0]
  'gold',                  [255,215,0]
  'orange',                [255,165,0]
  'darkorange',            [255,140,0]
  'lightyellow',           [255,255,224]
  'lemonchiffon',          [255,250,205]
  'lightgoldenrodyellow',  [250,250,210]
  'papayawhip',            [255,239,213]
  'moccasin',              [255,228,181]
  'peachpuff',             [255,218,185]
  'palegoldenrod',         [238,232,170]
  'khaki',                 [240,230,140]
  'darkkhaki',             [189,183,107]
  'yellow',                [255,255,0]
  'lawngreen',             [124,252,0]
  'chartreuse',            [127,255,0]
  'limegreen',             [50,205,50]
  'lime',                  [0,255,0]
  'forestgreen',           [34,139,34]
  'green',                 [0,128,0]
  'darkgreen',             [0,100,0]
  'greenyellow',           [173,255,47]
  'yellowgreen',           [154,205,50]
  'springgreen',           [0,255,127]
  'mediumspringgreen',     [0,250,154]
  'lightgreen',            [144,238,144]
  'palegreen',             [152,251,152]
  'darkseagreen',          [143,188,143]
  'mediumseagreen',        [60,179,113]
  'seagreen',              [46,139,87]
  'olive',                 [128,128,0]
  'darkolivegreen',        [85,107,47]
  'olivedrab',             [107,142,35]
  'lightcyan',             [224,255,255]
  'cyan',                  [0,255,255]
  'aqua',                  [0,255,255]
  'aquamarine',            [127,255,212]
  'mediumaquamarine',      [102,205,170]
  'paleturquoise',         [175,238,238]
  'turquoise',             [64,224,208]
  'mediumturquoise',       [72,209,204]
  'darkturquoise',         [0,206,209]
  'lightseagreen',         [32,178,170]
  'cadetblue',             [95,158,160]
  'darkcyan',              [0,139,139]
  'teal',                  [0,128,128]
  'powderblue',            [176,224,230]
  'lightblue',             [173,216,230]
  'lightskyblue',          [135,206,250]
  'skyblue',               [135,206,235]
  'deepskyblue',           [0,191,255]
  'lightsteelblue',        [176,196,222]
  'dodgerblue',            [30,144,255]
  'cornflowerblue',        [100,149,237]
  'steelblue',             [70,130,180]
  'royalblue',             [65,105,225]
  'blue',                  [0,0,255]
  'mediumblue',            [0,0,205]
  'darkblue',              [0,0,139]
  'navy',                  [0,0,128]
  'midnightblue',          [25,25,112]
  'mediumslateblue',       [123,104,238]
  'slateblue',             [106,90,205]
  'darkslateblue',         [72,61,139]
  'lavender',              [230,230,250]
  'thistle',               [216,191,216]
  'plum',                  [221,160,221]
  'violet',                [238,130,238]
  'orchid',                [218,112,214]
  'fuchsia',               [255,0,255]
  'magenta',               [255,0,255]
  'mediumorchid',          [186,85,211]
  'mediumpurple',          [147,112,219]
  'blueviolet',            [138,43,226]
  'darkviolet',            [148,0,211]
  'darkorchid',            [153,50,204]
  'darkmagenta',           [139,0,139]
  'purple',                [128,0,128]
  'indigo',                [75,0,130]
  'pink',                  [255,192,203]
  'lightpink',             [255,182,193]
  'hotpink',               [255,105,180]
  'deeppink',              [255,20,147]
  'palevioletred',         [219,112,147]
  'mediumvioletred',       [199,21,133]
  'white',                 [255,255,255]
  'snow',                  [255,250,250]
  'honeydew',              [240,255,240]
  'mintcream',             [245,255,250]
  'azure',                 [240,255,255]
  'aliceblue',             [240,248,255]
  'ghostwhite',            [248,248,255]
  'whitesmoke',            [245,245,245]
  'seashell',              [255,245,238]
  'beige',                 [245,245,220]
  'oldlace',               [253,245,230]
  'floralwhite',           [255,250,240]
  'ivory',                 [255,255,240]
  'antiquewhite',          [250,235,215]
  'linen',                 [250,240,230]
  'lavenderblush',         [255,240,245]
  'lightgray',             [211,211,211]
  'silver',                [192,192,192]
  'darkgray',              [169,169,169]
  'gray',                  [128,128,128]
  'dimgray',               [105,105,105]
  'lightslategray',        [119,136,153]
  'slategray',             [112,128,144]
  'darkslategray',         [47,79,79]
  'black',                 [0,0,0]
  'cornsilk',              [255,248,220]
  'blanchedalmond',        [255,235,205]
  'bisque',                [255,228,196]
  'navajowhite',           [255,222,173]
  'wheat',                 [245,222,179]
  'burlywood',             [222,184,135]
  'tan',                   [210,180,140]
  'rosybrown',             [188,143,143]
  'sandybrown',            [244,164,96]
  'goldenrod',             [218,165,32]
  'peru',                  [205,133,63]
  'chocolate',             [210,105,30]
  'saddlebrown',           [139,69,19]
  'sienna',                [160,82,45]
  'brown',                 [165,42,42]
  'maroon',                [128,0,0]
  };

if nargin==0
  % return the list with names
  y = mapping(:,1);
elseif ischar(x)
  % look up the corresponding RGB values
  sel = strcmp(mapping(:,1), x);
  y = mapping{sel, 2};
  y = y/255;
elseif isnumeric(x)
  x = x*255;
  % look up the corresponding name
  rgb = vertcat(mapping{:,2});
  rgb(:,1) = rgb(:,1) - x(1);
  rgb(:,2) = rgb(:,2) - x(2);
  rgb(:,3) = rgb(:,3) - x(3);
  [dum, sel] = min(sum(rgb.^2, 2));
  y = mapping{sel,1};
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [mesh, meshCurvature, layerCurvature, gLocs2d, gLocs3d, numNodes] ...
    = mfmAssignGray(mesh, params, insideNodes)
% Assign gray matter locations to the unfolded positions in a flattened
% map.
%
% [mesh, meshCurvature, layerCurvature, gLocs2d, gLocs3d, numNodes] ...
%    = mfmAssignGray(mesh, params, insideNodes)
%
% Author:  Winawer
%
%   Purpose: Assign gray matter locations to the unfolded positions in a flattened
%   map.
%
%   Sub-routine derived from Alex's unfoldMeshFromGUI code.
%
% See Also:  unfoldMeshFromGUI
% Get all the variables we may need
meshFileName     = params.meshFileName;
grayFileName     = params.grayFileName;
flatFileName     = params.flatFileName;
startCoords      = params.startCoords;
scaleFactor      = params.scaleFactor;
perimDist        = params.perimDist;
statusHandle     = params.statusHandle;
busyHandle       = params.busyHandle;
spacingMethod    = params.spacingMethod;
adjustSpacing    = params.adjustSpacing;
gridSpacing      = params.gridSpacing;
showFigures      = params.showFigures;
saveExtra        = params.saveExtra;
truePerimDist    = params.truePerimDist;
hemi             = params.hemi;
nperims          = params.NPERIMS;
saveIntermeidate = params.SAVE_INTERMEDIATE;
numberOfSteps    = params.NUMBEROFSTEPS;

% Update status
statusStringAdd(statusHandle,'Reading gray graph...');
[gNodes, gEdges] = readGrayGraph(grayFileName, hemi);

% Do all the gray node stages in a loop so that we can have arbitrary
% numbers of layers.
numGrayLayers=max(gNodes(6,:));

% Get the indices for all the gnodes (all layers)
for t=1:numGrayLayers
    grayNodeIndices{t}=find(gNodes(6,:)==t);
end

% Extract the layer 1 nodes. These are special because they are used to map
% down to the unfolded boundary mesh.
l1gNodes=gNodes(:,grayNodeIndices{1});
l1mesh.vertices=l1gNodes(1:3,:);
l1mesh.indices=grayNodeIndices{1};

% How many gNodes are there?
nGnodes=length(gNodes);
% How many gEdges are there?
nGedges=length(gEdges);

% We want to make a grey node connection matrix - which grey nodes are connected to which other gnodes?
statusStringAdd(statusHandle,'Finding grey connection matrix (slow)');

grayConMat=makeGrayConMat(gNodes,gEdges,busyHandle);

% We can assign layer 1 grey nodes to the white matter mesh using
% assignToNearest.dll (see assignToNearest.c) Can't do this for higher
% levels of grey matter because they might be mis-assigned. (Also, potential
% problem near  very crinkly edges. - Could we accidentally assign a l1
% grey matter node to the wrong WM point?) (I think the answer is 'yes,
% rarely'. If a single layer or gray is sandwiched between two sides of a
% sulcus (say) : It's grown from one side but mrFlatMesh has no way of
% telling which one.  Going deeper into mrGray's source code to determine
% the parentage of l1 gray nodes might be possible...)

% So for higher grey matter points, we have to restrict the possible sub
% node search space by assigning them >only< to  points they are connected
% to. Note that a single layer2 grey node may be connected to several l1
% nodes

% The gray may be defined over the entire mesh but we only want to deal with gray points over the
% unfolded part. The strategy should be....
% 1) do assignToNearest for each mesh point to find the nearest connected l1 node
% 2) Use the set of l1 nodes found in 1) to build up a list of other connected gray nodes
% 3) Repeat stage 2 for l3,l4

statusStringAdd(statusHandle,'Mapping L1 to mesh.');

% Find 3D coords of all the l1 gnodes
l1GNodeCoords=l1gNodes(1:3,:)';

% Find 3D coords of all the mesh points (not just the unfolded ones)
meshCoords=mesh.uniqueVertices;

% Mesh coords are in voxels, gnode coords are in voxels. Bring them into alignment
l1GNodeCoordsN=l1GNodeCoords-1; % Use this -1 to bring them into perfect alignment. This was the source of an error in V1.0

% Now restrict the mesh coords to just those points in the unfold
meshCoordsN=meshCoords(mesh.insideNodes,:);

% And now restrict the set of l1 gray nodes so that only those that are relatively near the
% mesh are included in the search - this is done first as a simple bounds check
boundedL1NodeIndices=boundsCheck3D(min(meshCoordsN)-3,max(meshCoordsN)+3,l1GNodeCoordsN);
boundedL1GNodes=l1GNodeCoordsN(boundedL1NodeIndices,:);
boundedL1NodeIndices=grayNodeIndices{1}(boundedL1NodeIndices); % This is now a list of indices into the full gNode array

statusStringAdd(statusHandle,'Finding nearest Ll gray points to mesh (very slow)');

% Find mesh points near the l1 gNodes
[boundedL1toMeshIndices,sqrDist]=nearpoints(boundedL1GNodes',(meshCoordsN)');
% This returns a list of indices into the meshCoordsN array that links a single 3D mesh point to each l1Gnode)
% and then eliminate any l1gNodes that are more than a set distance away from the mesh - here sqrt(0.25)

% *************************
closeEnough=find(sqrDist<=0.2501);
% *************************

boundedL1toMeshNodes  = boundedL1GNodes(closeEnough,:); % remember, assignToNearest returns the squared distance
boundedL1toMeshIndices= boundedL1toMeshIndices(closeEnough);

% For each member of the bounded l1Gnodes, this tells us the index of the full mesh point that it maps to.
fullBoundedL1toMeshIndices=insideNodes(boundedL1toMeshIndices);

restNodeIndices{1}=boundedL1NodeIndices(closeEnough);

statusStringAdd(statusHandle,'Setting L1 glocs');
% We can start setting gLocs
% Note the +1 here! Because we subtracted 1 from the original node coords
% to align them with the mesh.
layerGlocs3d{1}=boundedL1toMeshNodes+1;
layerGlocs2d{1}=mesh.X(fullBoundedL1toMeshIndices,:);

numNodes(1)=length(layerGlocs3d{1});

layerCurvature{1}=mesh.uniqueCols(fullBoundedL1toMeshIndices,1);

if (showFigures)
    h = unfoldPlotL1Mesh(grayConMat(restNodeIndices{1},restNodeIndices{1}), layerGlocs2d{1});
    statusStringAdd(statusHandle,'Displaying L1 gray mesh.');

	[p f ext] = fileparts(flatFileName);
	savePath = fullfile(p, [f ' Gray Mesh.png']);
	saveas(h, savePath);
	fprintf('[%s]: Saved %s.\n', mfilename, savePath);
end

% Now we have to find l2tol1Indices, l3tol2Indices and l4tol3Indices. This
% is faster since for each point, we restrict its  potential nearest
% neighbours to points that it is connected to in the previous layer.  We
% also restrict the l2 nodes to just those that are connected to the
% restricted l1 nodes and the l3 nodes to those connected to the l2 nodes.
% Use the full connection matrix to find which l2 Gnodes are connected to the restricted l1Gnodes
statusStringAdd(statusHandle,'Mapping other gray layers (2,3,4..)');

% Set up the 3 critical arrays curvature, gLocs2d, gLocs3d. These are what
% eventually get written to the flat.
meshCurvature=layerCurvature{1};
gLocs2d=layerGlocs2d{1};
gLocs3d=layerGlocs3d{1};

for t=2:numGrayLayers

    layerInterconnect=grayConMat(:, restNodeIndices{t-1}); % Finds all the things connected to the t-1th layer...
    [restNodeIndices{t},dummy]=find(layerInterconnect);
    restNodeIndices{t}=unique(restNodeIndices{t});
    restNodeIndices{t}=intersect(restNodeIndices{t},grayNodeIndices{t}); % Only take potential l2 node indices
    finalConnectIndices{t-1}=findNearestConnected(gNodes', restNodeIndices{t},restNodeIndices{t-1},grayConMat);

    % Set glocs3d, glocs2d, meshCurvature
    layerGlocs3d{t}=gNodes(1:3,restNodeIndices{t})';
    numNodes(t)=length(layerGlocs3d{t});
    layerGlocs2d{t}=layerGlocs2d{t-1}(finalConnectIndices{t-1},:);
    layerCurvature{t}=zeros(length(finalConnectIndices{t-1}),1);
    gLocs2d=[gLocs2d;layerGlocs2d{t}];
    gLocs3d=[gLocs3d;layerGlocs3d{t}];
    meshCurvature=[meshCurvature;layerCurvature{t}];

end

if (showFigures),
	unfoldPlotgLocs2d(numGrayLayers,layerGlocs2d);
	[p f ext] = fileparts(flatFileName);
	savePath = fullfile(p, [f ' Grid Locations.png']);
	saveas(gcf, savePath);
	fprintf('[%s]: Saved %s.\n', mfilename, savePath);
end

statusStringAdd(statusHandle,'Creating flat.mat structure');

% We want to transform  the curvature to run between 1 to 64 for display
% purposes.  We want the mean, however, to be preserved so that we can
% calculate the average curvature and also we want to keep the values
% accurate in terms of real curvature (which is 1/radius in pixels but
% should be mm; though mrGray may do some unwanted scaling).
% So, instead of this:
%     meshCurvature=normalize((meshCurvature))*63+1;%
% we now do this to (a) map to 1 to 64,
%
% Note that the curvature values here are pulled from the rgba field of the
% mesh sturct. This is created by meGray and is the actual RGB (& alpha)
% value used to color each triangle when the mesh is rendered. If the mesh
% has no color overlays, then this will be an 8-bit grayscale value (R=G=B)
% and thus range from 0-255, with curvature=0 at 127. (RFD/BW).
meshCurvature = (layerCurvature{1}/255)*63 ;

return

%----------------------------------
% ---- SUBROUTINES
%----------------------------------
%----------------------------------
function h = unfoldPlotL1Mesh(x,layer1Glocs2d)
h = figure;
gplot(x,layer1Glocs2d);
title('L1 gray mesh'); zoom on;

return;

%----------------------------------
function unfoldPlotgLocs2d(numGrayLayers,layerGlocs2d)
%
nFigsAcross=ceil(sqrt(numGrayLayers));
nFigsDown=ceil(sqrt(numGrayLayers));
for t=1:numGrayLayers
    if (t<=numGrayLayers)
        subplot(nFigsDown,nFigsAcross,t);

        plot(layerGlocs2d{t}(:,1),layerGlocs2d{t}(:,2),'.');
        axis equal;
    end % end check on current layer index
end % Next image
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function SpharmMatExpandAlignPDM(confs, objs)
numSbj = length(objs);
if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end
if (~exist([confs.OutDirectory '/Logs'],'dir'))
    mkdir([confs.OutDirectory '/Logs']);
end

h = waitbar(0,'Please wait...');
for i = 1:numSbj
    file1 = objs{i};
    [pa, name, ext] = fileparts(file1);
    file2 = [pa '/' name(1:end-4) 'para' ext];
    diary(fullfile([confs.OutDirectory '/Logs'],[date, '_', num2str(round(cputime)), '_' name '_ExpPDM.log']));

    if ~exist(file2,'file')
        disp(sprintf('Corresponding parameter file to %s doesnot exist',name));
    else

        optionSTR = '';
        for j = 1:length(confs.vars)
            vals = eval(sprintf('confs.%s', confs.vars{j}));
            if ~isempty(vals) & (confs.args(j) < 10)
                vals = num2str(vals);
                optionSTR = [optionSTR ' ' sprintf('-%s ', confs.vars{j}) vals];
            elseif (confs.args(j) < 200) & (confs.args(j) >= 10) & ~isempty(deblank(vals)) & ~strcmp(confs.vars{j}, 'others')
                optionSTR = [optionSTR ' ' sprintf('-%s ', confs.vars{j}) vals];
            elseif (confs.args(j) < 200) & (confs.args(j) >= 10) & ~isempty(deblank(vals)) & strcmp(confs.vars{j}, 'others')
                optionSTR = [optionSTR ' ' vals];
            end
        end

        if ispc
            cSTR = sprintf('!%s %s %s %s', confs.command, file1, file2, optionSTR);
            disp('This function is not working on MS Windows. Please use Linux version.');
        elseif isunix
            preload = '/lib64/libgcc_s.so.1';
            if exist(preload, 'file')
                setenv('LD_PRELOAD', preload);
            else
                setenv('LD_PRELOAD', 'lib/libcc_s_so.1');
            end
            cSTR = sprintf('!./%s %s %s %s', confs.command, file1, file2, optionSTR);
            curPath = pwd;
            cd(confs.path);
                eval(cSTR);
            cd(curPath);
            movefile([pa '/' name 'SPHARM*.meta'],confs.OutDirectory);
            movefile([pa '/' name 'SPHARM*.coef'],confs.OutDirectory);
        elseif ismac
            disp('Not implemented yet.');
        end
    end
    diary('off');

    waitbar(i/numSbj)
end
close(h);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function ellipse_grid_display ( ng, xy )

%*****************************************************************************80
%
%% ELLIPSE_GRID_DISPLAY displays grid points inside an ellipse.
%
%  Discussion:
%
%    The ellipse is specified as
%
%      ( ( X - C1 ) / R1 )^2 + ( ( Y - C2 ) / R2 )^2 = 1
%
%    The user supplies a number N.  There will be N+1 grid points along
%    the shorter axis.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NG, the number of grid points inside the ellipse.
%
%    Input, real XY(2,NG), the grid points.
%
  scatter ( xy(1,:), xy(2,:), 'b.' );
  axis equal
  title ( sprintf ( '%d grid points inside an ellipse', ng ) )
  grid on

  return
end


text:
% plot using mapping toolbox from shapefile
% download the shapefile from http://www.gadm.org (recommende is the one with 6 dissolved layers
%S=shaperead('/Users/reyesc/Downloads/gadm28_levels/gadm28_adm1.shp','UseGeoCoords',true); %1st administrative level
S=shaperead('/Users/reyesc/Downloads/gadm28_levels/gadm28_adm0.shp','UseGeoCoords',true); %country administrative level
figure
axesm lambert
framem
plotm([S.Lat],[S.Lon]);

%plot(ax,[S.Lon],[S.Lat],'color',[.5 .5 .5],'Tag','borders')

author: CelsoReyes


text:
function [stats,Y,X] = spm_mci_mvnpost (post,method,verbose,max_lag)
% Are MCMC samples consistent with Gaussian posterior ?
% FORMAT [stats,Y,X] = spm_mci_mvnpost (post,method,verbose,max_lag)
%
% post      posterior data structure from spm_mci_post
% method    'ESS' or 'thinning'
% verbose   create plots
% max_lag   maximum potential lag for MAR model
%
% stats     (multivariate) normal test statistics
%           See spm_mci_mvntest.m
% Y         uncorrelated posterior samples
% X         original posterior samples
%
% Run Gaussianity test on Markov chain samples
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_mvnpost.m 6697 2016-01-27 14:57:28Z spm $

try, pmax=max_lag; catch, pmax=10; end
try, meth=method; catch, meth='ESS'; end
try, ver=verbose; catch, ver=1; end

j=post.ind;
X=post.P(:,j)';
Nj=length(j);
Np=size(post.P,1);

if any(std(X)==0)
    disp('Warning from spm_mci_mvnpost: some parameters have zero variance');
    stats=[];
    return
end

switch meth
    case 'ESS',
        %C=cov(post.P(:,post.ind)');
        for p=1:Np,
            ess(p)=spm_mci_ess(post.P(p,j));
            %v(p)=C(p,p)*ess(p)/Nj;
        end
        %post.Cp=diag(v);
        %mess=mean(ess);
        mess=ceil(min(ess));
        stats = spm_mci_mvntest(X,mess);

        if ver
            disp('ESS method');
            disp('Effective Sample Sizes:');
            disp(ess);
        end

    case 'thinning',
        for p=1:Np,
            [tmp,m(p)]=spm_mci_ess(post.P(p,j));
        end
        lag=max(m)+1;
        Y=X(1:lag:end,:);
        stats = spm_mci_mvntest(Y);
        if ver
            disp('Thinning method');
            disp(sprintf('Using only every %d-th sample',lag));
        end

    otherwise
        disp('Unknown method in spm_mci_mvnpost.m');
        return
end

author: spm


text:
function V = mttkrp(X,U,n)
%MTTKRP Matricized tensor times Khatri-Rao product for sparse tensor.
%
%   V = MTTKRP(X,U,n) efficiently calculates the matrix product of the
%   n-mode matricization of X with the Khatri-Rao product of all
%   entries in U, a cell array of matrices, except the nth.  How to
%   most efficiently do this computation depends on the type of tensor
%   involved.
%
%   See also SPTENSOR, TENSOR/MTTKRP, SPTENSOR/TTV
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

% In the sparse case, it is most efficient to do a series of TTV operations
% rather than forming the Khatri-Rao product.

N = ndims(X);

if (n == 1)
    R = size(U{2},2);
else
    R = size(U{1},2);
end

V = zeros(size(X,n),R);
for r = 1:R
    % Set up cell array with appropriate vectors for ttv multiplication
    Z = cell(N,1);
    for i = [1:n-1,n+1:N]
        Z{i} = U{i}(:,r);
    end
    % Perform ttv multiplication
    V(:,r) = double(ttv(X, Z, -n));
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function smaps =  msCreateMultipleSegmentations(pE, adjlist, nsp, nsegall)
% 1) Randomly select superpixel s1, then randomly selects different superpixel
% s2 within same segment (if one exists); remove s1,s2 from s
% 2) Then, for randomly ordered i:
%     if si is adjacent to s1, assign si to s1 with probability pE(si, s1)
%     if si is adjacent to s2, assign si to s2 with probability pE(si, s2)
%     remove si from s
% 3) Repeat (2) until s is empty
%
% Note: this differs from the published version of the algorithm in that
% duplicate segments are removed (number is replaced with 0 in smaps)

% randomly split each segment further into connected components
%adjmat2 = curradjmat;
nmaps = numel(nsegall);

smaps = zeros(nsp, nmaps);

adjmat = zeros(nsp, nsp);
for k = 1:size(adjlist, 1)
    s1 = adjlist(k, 1);
    s2 = adjlist(k, 2);
    adjmat(s1, s2) = k;
    adjmat(s2, s1) = k;
end

nadj = zeros(nsp, 1);
for k = 1:nsp
    adj{k} = find(adjmat(k, :));
    nadj(k) = numel(adj{k});
end
for k = 1:nsp
    normalization = 1;%2./(nadj(k)+nadj(adj{k}));
    pOffAll{k} = (1-pE(adjmat(k, adj{k}))).^normalization + 1E-10;
    pOnAll{k} = pE(adjmat(k, adj{k})).^normalization + 1E-10;
end

for m = 1:nmaps
    rind = randperm(nsp);
    smap = zeros(nsp, 1);
    nseg = nsegall(m);

    nseg = min(nseg, nsp);
    smap(rind(1:nseg)) = (1:nseg);

    if nseg==nsp % more segments than superpixels: return on sp per segment
        smaps(:, m) = smap;
        continue;
    end

    rind(1:nseg) = [];

    nleft = nsp;
    %tic
    while sum(nadj(rind))>0 % do until all possible sp are assigned
        for r = rind
            %p = -Inf*ones(nseg, 1);

            if all(smap(adj{r})>0)
                p = -Inf*ones(nseg, 1);
                for k = unique(smap(adj{r}))'
                    %if any(smap(adj{r})==k)
                        sOn = (smap(adj{r})==k);
                        sOff = ~sOn; %setdiff(find(smap(adj{r})>0), sOn);
                        pOn = prod(pOnAll{r}(sOn));
                        pOff = 1;
                        if ~isempty(sOff)
                            pOff = prod(pOffAll{r}(sOff));
                        end
                        p(k) = log(pOn)+log(pOff);
                    %end
                end
                [kval, kmax] = max(p);
                smap(r) = kmax;
            else
                for k = unique(smap(adj{r}))' %1:nseg
                    if k > 0
                        kadj = (smap(adj{r})==k);
                    %if sum(kadj)>0
                        pOn = prod(pOnAll{r}(kadj));
                        pOff = prod(pOffAll{r}(kadj));
                        if rand(1) < pOn / (pOn+pOff) + 0.05
                            smap(r) = k;
                        end
                    %end
                    end
                end
            end
        end

        rind(smap(rind)>0) = [];

        % break if it seems stuck (can happen when pE is very low,
        % nsegments is small)
        if numel(nleft)>50 && all(nleft(end-49:end)==numel(rind))
            %disp(['break init early: ' num2str(numel(rind)) ' left'])
            %disp(sort(rind))
            break;
        end
        nleft(end+1) = numel(rind);

    end

    count = zeros(nseg, 1);
    for k = 1:nseg
        count(k) = sum(smap==k);
    end
    %toc
    %disp(['0: ' num2str(evaluateEdgeProb(adjlist, pE, smap))])
    %tic
    for t = 1:10
        lastmap = smap;
        rind = randperm(nsp);
        for r = rind
            if smap(r) > 0 && count(smap(r)) > 1
                p = -Inf*ones(nseg, 1);

                for k = unique(smap(adj{r}))' %1:nseg
                    if k > 0
                    %if sum(smap(adj{r})==k)>0
                        sOn = smap(adj{r})==k;
                        sOff = ~sOn; %setdiff([1:nadj(r)], sOn);
                        pOn = sum(log(pOnAll{r}(sOn)));
                        pOff = sum(log(pOffAll{r}(sOff)));
                        p(k) = pOn + pOff;
                    end
                    %end
                end
                [kval, kmax] = max(p);
                if kmax~=smap(r)
                    count(smap(r)) = count(smap(r))-1;
                    smap(r) = kmax;
                end

            end
        end
        if all(lastmap==smap)
            break;
        end
    end
    %disp(['1: ' num2str(evaluateEdgeProb(adjlist, pE, smap))])
    %toc
    %disp(num2str(t))
    smaps(:, m) = smap;
end

smaps = msPruneSegments(smaps);

%%

function p = evaluateEdgeProb(adjlist, pE, smap)
hasedge = zeros(size(adjlist, 1), 1);
for k = 1:size(adjlist, 1)
    if smap(adjlist(k, 1))==smap(adjlist(k, 2))
        hasedge(k)=1;
    end
end
p = sum(log(pE(hasedge==1)));
p = p+sum(log(1-pE(hasedge==0)));

author: zouchuhang

repo: LayoutNet


text:
function z = languageModelHead(z,p,word_embeddings)
% languageModelHead   The standard BERT operations for masked
% language modeling.
%
%   Z = languageModelHead(X,languageModelWeights,wordEmbeddingMatrix)
%   applies the language model operations to an input X. The input X must
%   be an unformatted dlarray of size
%   hiddenSize-by-numInputTokens-by-numObs. The languageModelWeights must
%   be a struct with fields 'transform' and 'LayerNorm' such as the
%   mdl.Parameters.Weights.masked_LM struct where mdl = bert(). The
%   wordEmbeddingMatrix must be the word embedding matrix used by the bert
%   model such as mdl.Parameters.Weights.embeddings.word_embeddings where
%   mdl = bert().

% Copyright 2021 The MathWorks, Inc.
z = transformer.layer.convolution1d(z,p.transform.kernel,p.transform.bias);
z = transformer.layer.gelu(z);
z = transformer.layer.normalization(z,p.LayerNorm.gamma,p.LayerNorm.beta);
z = dlmtimes(word_embeddings.',z);
z = z + p.output.bias;
z = softmax(z,'DataFormat','CTB');
end

author: matlab-deep-learning

repo: transformer-models


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

for c = 1:num_labels
    init_theta = zeros(n + 1, 1);
    options = optimset('GradObj', 'on', 'MaxIter', 50);
    [theta] = fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), init_theta, options);
    all_theta(c,:) = theta';
end;

% =========================================================================

end

repo: MachineLearning


text:
%#codegen
function [pixel_val, pixel_valid] = aMediantFilter_2D(c_data, c_idx)

smax = 9;
persistent window;
if isempty(window)
    window = zeros(smax, smax);
end

cp = ceil(smax/2); % center pixel;

w3 = -1:1;
w5 = -2:2;
w7 = -3:3;
w9 = -4:4;

r3 = cp + w3;      % 3x3 window
r5 = cp + w5;      % 5x5 window
r7 = cp + w7;      % 7x7 window
r9 = cp + w9;      % 9x9 window

d3x3 = window(r3, r3);
d5x5 = window(r5, r5);
d7x7 = window(r7, r7);
d9x9 = window(r9, r9);

center_pixel = window(cp, cp);

% use 1D filter for 3x3 region
outbuf = get_median_1d(d3x3(:)');
[min3, med3, max3] = getMinMaxMed_1d(outbuf);

% use 2D filter for 5x5 region
outbuf = get_median_2d(d5x5);
[min5, med5, max5] = getMinMaxMed_2d(outbuf);

% use 2D filter for 7x7 region
outbuf = get_median_2d(d7x7);
[min7, med7, max7] = getMinMaxMed_2d(outbuf);

% use 2D filter for 9x9 region
outbuf = get_median_2d(d9x9);
[min9, med9, max9] = getMinMaxMed_2d(outbuf);

pixel_val = get_new_pixel(min3, med3, max3, ...
    min5, med5, max5, ...
    min7, med7, max7, ...
    min9, med9, max9, ...
    center_pixel);

% we need to wait until 9 cycles for the buffer to fill up
% output is not valid every time we start from col1 for 9 cycles.
persistent datavalid
if isempty(datavalid)
    datavalid = false;
end
pixel_valid = datavalid;
datavalid = (c_idx >= smax);

% build the 9x9 buffer
window(:,2:smax) = window(:,1:smax-1);
window(:,1) = c_data;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [min, med, max] = getMinMaxMed_1d(inbuf)

max = inbuf(1);
med = inbuf(ceil(numel(inbuf)/2));
min = inbuf(numel(inbuf));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [min, med, max] = getMinMaxMed_2d(inbuf)

[nrows, ncols] = size(inbuf);
max = inbuf(1, 1);
med = inbuf(ceil(nrows/2), ceil(ncols/2));
min = inbuf(nrows, ncols);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function new_pixel  = get_new_pixel(...
    min3, med3, max3, ...
    min5, med5, max5, ...
    min7, med7, max7, ...
    min9, med9, max9, ...
    center_data)

if (med3 > min3 && med3 < max3)
    new_pixel = get_center_data(min3, med3, max3,center_data);
elseif (med5 > min5 && med5 < max5)
    new_pixel = get_center_data(min5, med5, max5,center_data);
elseif (med7 > min7 && med7 < max7)
    new_pixel = get_center_data(min7, med7, max7,center_data);
elseif (med9 > min9 && med9 < max9)
    new_pixel = get_center_data(min9, med9, max9,center_data);
else
    new_pixel = center_data;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [new_data] = get_center_data(min,med,max,center_data)
if center_data <= min || center_data >= max
    new_data = med;
else
    new_data = center_data;
end
end

author: Sable

repo: mcbench-benchmarks


text:
function [net,stats] = cnn_train_dag(net, imdb, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;
opts.saveMomentum = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = true;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end
if isnan(opts.val), opts.val = [] ; end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
    if isempty(opts.derOutputs)
        error('DEROUTPUTS must be specified when training.\n') ;
    end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
    fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
    [net, state, stats] = loadState(modelPath(start)) ;
else
    state = [] ;
end

for epoch=start+1:opts.numEpochs

    % Set the random seed based on the epoch and opts.randomSeed.
    % This is important for reproducibility, including when training
    % is restarted from a checkpoint.

    rng(epoch + opts.randomSeed) ;
    prepareGPUs(opts, epoch == start+1) ;

    % Train for one epoch.
    params = opts ;
    params.epoch = epoch ;
    params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
    params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
    params.val = opts.val(randperm(numel(opts.val))) ;
    params.imdb = imdb ;
    params.getBatch = getBatch ;

    if numel(opts.gpus) <= 1
        [net, state] = processEpoch(net, state, params, 'val') ;

        [net, state] = processEpoch(net, state, params, 'train') ;
        [net, state] = processEpoch(net, state, params, 'val') ;
        if ~evaluateMode
            saveState(modelPath(epoch), net, state) ;
        end
        lastStats = state.stats ;
    else
        spmd
            [net, state] = processEpoch(net, state, params, 'train') ;
            [net, state] = processEpoch(net, state, params, 'val') ;
            if labindex == 1 && ~evaluateMode
                saveState(modelPath(epoch), net, state) ;
            end
            lastStats = state.stats ;
        end
        lastStats = accumulateStats(lastStats) ;
    end

    stats.train(epoch) = lastStats.train ;
    stats.val(epoch) = lastStats.val ;
    clear lastStats ;
    saveStats(modelPath(epoch), stats) ;

    if opts.plotStatistics
        switchFigure(1) ; clf ;
        plots = setdiff(...
            cat(2,...
            fieldnames(stats.train)', ...
            fieldnames(stats.val)'), {'num', 'time'}) ;
        for p = plots
            p = char(p) ;
            values = zeros(0, epoch) ;
            leg = {} ;
            for f = {'train', 'val'}
                f = char(f) ;
                if isfield(stats.(f), p)
                    tmp = [stats.(f).(p)] ;
                    values(end+1,:) = tmp(1,:)' ;
                    leg{end+1} = f ;
                end
            end
            subplot(1,numel(plots),find(strcmp(p,plots))) ;
            plot(1:epoch, values','o-') ;
            xlabel('epoch') ;
            title(p) ;
            legend(leg{:}) ;
            grid on ;
        end
        drawnow ;
        print(1, modelFigPath, '-dpdf') ;
    end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.momentum)
    state.momentum = num2cell(zeros(1, numel(net.params))) ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
    net.move('gpu') ;
    state.momentum = cellfun(@gpuArray, state.momentum, 'uniformoutput', false) ;
end
if numGpus > 1
    parserv = ParameterServer(params.parameterServer) ;
    net.setParameterServer(parserv) ;
else
    parserv = [] ;
end

% profile
if params.profile
    if numGpus <= 1
        profile clear ;
        profile on ;
    else
        mpiprofile reset ;
        mpiprofile on ;
    end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
    fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
        fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
    batchSize = min(params.batchSize, numel(subset) - t + 1) ;

    for s=1:params.numSubBatches
        % get this image batch and prefetch the next
        batchStart = t + (labindex-1) + (s-1) * numlabs ;
        batchEnd = min(t+params.batchSize-1, numel(subset)) ;
        batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
        num = num + numel(batch) ;
        if numel(batch) == 0, continue ; end

        inputs = params.getBatch(params.imdb, batch) ;

        if params.prefetch
            if s == params.numSubBatches
                batchStart = t + (labindex-1) + params.batchSize ;
                batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
            else
                batchStart = batchStart + numlabs ;
            end
            nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
            params.getBatch(params.imdb, nextBatch) ;
        end

        if strcmp(mode, 'train')
            net.mode = 'normal' ;
            net.accumulateParamDers = (s ~= 1) ;
            net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
        else
            net.mode = 'test' ;
            net.vars(11).precious = 1;
            net.eval(inputs) ;
        end

        % test in batch
%         imt = inputs{2};
%         labels = inputs{end};
%         score_batch = squeeze(net.vars(11).value);
%         [~, predLabel_batch] = max(score_batch, [], 1);
%         acc_batch = 1-mean(labels(:)==predLabel_batch(:));
%         % test individually
%         predLabel_single = zeros(1, length(labels));
%         score_single = zeros(10, length(labels));
%         for iii = 1:length(labels)
%             curLabel = labels(iii);
%             curInput = {inputs{1}, imt(:,:,:,iii), inputs{3}, curLabel};
%
%             net.vars(11).precious = 1;
%             net.eval(curInput);
%
%             curA = squeeze(net.vars(11).value);
%             score_single(:, iii) = curA(:);
%             [~, curPredLabel] = max(curA, [], 1);
%             predLabel_single(iii) = curPredLabel;
%         end
%         acc_single = 1-mean(labels(:)==predLabel_single(:));
%
%         fprintf('\nacc in batch form: %.4f\n', acc_batch);
%         fprintf('acc in single form: %.4f\n', acc_single);
%         fprintf('score difference: %.4f\n', norm(score_single(:)-score_batch(:)));

    end

    % Accumulate gradient.
    if strcmp(mode, 'train')
        if ~isempty(parserv), parserv.sync() ; end
        state = accumulateGradients(net, state, params, batchSize, parserv) ;
    end

    % Get statistics.
    time = toc(start) + adjustTime ;
    batchTime = time - stats.time ;
    stats.num = num ;
    stats.time = time ;
    stats = params.extractStatsFn(stats,net) ;
    currentSpeed = batchSize / batchTime ;
    averageSpeed = (t + batchSize - 1) / time ;
    if t == 3*params.batchSize + 1
        % compensate for the first three iterations, which are outliers
        adjustTime = 4*batchTime - time ;
        stats.time = time + adjustTime ;
    end

    fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
    for f = setdiff(fieldnames(stats)', {'num', 'time'})
        f = char(f) ;
        fprintf(' %s: %.3f', f, stats.(f)) ;
    end
    fprintf('\n') ;
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
    if numGpus <= 1
        state.prof.(mode) = profile('info') ;
        profile off ;
    else
        state.prof.(mode) = mpiprofile('info');
        mpiprofile off ;
    end
end
if ~params.saveMomentum
    state.momentum = [] ;
else
    state.momentum = cellfun(@gather, state.momentum, 'uniformoutput', false) ;
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

    if ~isempty(parserv)
        parDer = parserv.pullWithIndex(p) ;
    else
        parDer = net.params(p).der ;
    end

    switch net.params(p).trainMethod

        case 'average' % mainly for batch normalization
            thisLR = net.params(p).learningRate ;
            net.params(p).value = vl_taccum(...
                1 - thisLR, net.params(p).value, ...
                (thisLR/batchSize/net.params(p).fanout),  parDer) ;

        case 'gradient'
            thisDecay = params.weightDecay * net.params(p).weightDecay ;
            thisLR = params.learningRate * net.params(p).learningRate ;

            if thisLR>0 || thisDecay>0
                % Normalize gradient and incorporate weight decay.
                parDer = vl_taccum(1/batchSize, parDer, ...
                    thisDecay, net.params(p).value) ;

                % Update momentum.
                state.momentum{p} = vl_taccum(...
                    params.momentum, state.momentum{p}, ...
                    -1, parDer) ;

                % Nesterov update (aka one step ahead).
                if params.nesterovUpdate
                    delta = vl_taccum(...
                        params.momentum, state.momentum{p}, ...
                        -1, parDer) ;
                else
                    delta = state.momentum{p} ;
                end

                % Update parameters.
                net.params(p).value = vl_taccum(...
                    1,  net.params(p).value, thisLR, delta) ;
            end
        otherwise
            error('Unknown training method ''%s'' for parameter ''%s''.', ...
                net.params(p).trainMethod, ...
                net.params(p).name) ;
    end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
    s = char(s) ;
    total = 0 ;

    % initialize stats stucture with same fields and same order as
    % stats_{1}
    stats__ = stats_{1} ;
    names = fieldnames(stats__.(s))' ;
    values = zeros(1, numel(names)) ;
    fields = cat(1, names, num2cell(values)) ;
    stats.(s) = struct(fields{:}) ;

    for g = 1:numel(stats_)
        stats__ = stats_{g} ;
        num__ = stats__.(s).num ;
        total = total + num__ ;

        for f = setdiff(fieldnames(stats__.(s))', 'num')
            f = char(f) ;
            stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

            if g == numel(stats_)
                stats.(s).(f) = stats.(s).(f) / total ;
            end
        end
    end
    stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
    stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
    save(fileName, 'stats', '-append') ;
else
    save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
    error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
    try
        set(0,'CurrentFigure',n) ;
    catch
        figure(n) ;
    end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
    % check parallel pool integrity as it could have timed out
    pool = gcp('nocreate') ;
    if ~isempty(pool) && pool.NumWorkers ~= numGpus
        delete(pool) ;
    end
    pool = gcp('nocreate') ;
    if isempty(pool)
        parpool('local', numGpus) ;
        cold = true ;
    end

end
if numGpus >= 1 && cold
    fprintf('%s: resetting GPU\n', mfilename)
    clearMex() ;
    if numGpus == 1
        gpuDevice(opts.gpus)
    else
        spmd
            clearMex() ;
            gpuDevice(opts.gpus(labindex))
        end
    end
end

author: aimerykong


text:
function [fnames, defaults] = mysubs_fields

% function [fnames, defaults] = mysubs_fields
% Additional fields for class cfg_branch. See help of
% @cfg_item/subs_fields for general help about this function.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: mysubs_fields.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

fnames={};
defaults={};


text:
function options = spset(varargin)
% SPSET   Create/alter sparse interpolation OPTIONS structure
%   OPTIONS = SPSET('NAME1',VALUE1,'NAME2',VALUE2,...) creates an
%   options structure OPTIONS in which the named properties have
%   the specified values. Any unspecified properties have default
%   values. It is sufficient to type only the leading characters
%   that uniquely identify the property. Case is ignored for
%   property names.
%
%   OPTIONS = SPSET(OLDOPTS,'NAME1',VALUE1,...) alters an existing
%   options structure OLDOPTS.
%
%   OPTIONS = SPSET(OLDOPTS,NEWOPTS) combines an existing options
%   structure OLDOPTS with a new options structure NEWOPTS. Any new
%   properties overwrite corresponding old properties.
%
%   SPSET with no input arguments displays all property names and
%   their possible values.
%
%   See also SPGET, SPINTERP, SPVALS.
%
%
% SPSET PROPERTIES
%
% GridType - Sparse grid type and basis functions to use
%   [ {Clenshaw-Curtis} | Maximum | NoBoundary | Chebyshev |
%   Gauss-Patterson ]. For an illustration of  the grid types,
%   run CMPGRIDS.
%
% RelTol - Relative error tolerance [ positive scalar {1e-2} ]
%   A relative error tolerance that applies to all hierarchical
%   surpluses of the sparse grid representation. The grid is
%   further refined until all hierarchical surpluses are less than
%   max(RelTol*(max(fevalRange)-min(fevalRange)),AbsTol), with
%   fevalRange containing all results evaluating FUN up to
%   that point.
%
% AbsTol - Absolute error tolerance [ positive scalar {1e-6} ]
%   Used by the error control stated under RelTol.
%
% Vectorized - Vectorized function FUN [ on | {off} ]
%   Indicates if FUN is available for vectorized evaluation.
%   Vectorized coding of FUN can significantly reduce the
%   computation time used by SPVALS. For an example using a
%   vectorized function, please see SPDEMO.
%
% MinDepth - Minimum interpolation depth [ integer {2} ]
%   Minimum number of hierarchical interpolation levels n to
%   compute.
%
% MaxDepth - Maximum interpolation depth [ integer {8} ]
%   Maximum number of hierarchical interpolation levels n to
%   compute. The maximum supported depth (and the default) is 6
%   for the Gauss-Patterson grid.
%
% VariablePositions - Position of the ranges to interpolate in the
%   argument list when FUN is evaluated. [ 1xD vector {[]} ]
%   By setting VariablePositions, SPVALS will evaluate FUN with respect
%   to some of its input parameters, but not necessarily the first
%   D ones. With VARPOS, the actual position of each variable may
%   be set. VARPOS must be a 1xD array.
%
% NumberOfOutputs - Number of outputs [ integer {1} ]
%   If FUN produces multiple outputs (all must be scalar), indicate
%   this here to perform the sparse grid computation for many
%   output variables at once. Also see the example SPDEMOVAROUT.
%
% PrevResults - Previous sparse grid representation [ struct {[]} ]
%   An existing result structure obtained from SPVALS may be
%   provided to further refine an existing sparse grid.
%
% FunctionArgType - Function argument type [ {list} | vector ]
%   Indicates whether the objective function takes the input
%   parameters as a comma-separated list (default) or as a vector.
%
% KeepFunctionValues - Keep the function evaluations [ {off} | on ]
%   If this parameter is set, a structure field 'fvals' is
%   returned, containing a cell array with the function values at
%   the sparse grid points.
%
% KeepGrid - Keep the sparse grid points [ {off} | on ]
%   If this parameter is set, a structure field 'grid' is
%   returned, containing a cell array with the the sparse grid
%   points.
%
% DimensionAdaptive - Set dimension adaptivity option [ {off} | on ]
%   Dimension-adaptive grids try to adaptively find important
%   dimensions and adjust the sparse grid structure
%   accordingly. The implementation of dimension-adaptive sparse
%   grids follows the approach by Gerstner/Griebel.
%
% MinPoints - Minimum number of support nodes [ integer {100} ]
%   This parameter only applies to dimension-adaptive sparse grids,
%   and indicates the minimum number of support nodes
%   (i.e. function evaluations to perform).
%
% MaxPoints - Maximum number of support nodes [ integer {10000} ]
%   This parameter only applies to dimension-adaptive sparse grids.
%   The dimension-adaptive algorithm is aborted once the function
%   evaluation count exceeds this number.
%
% DimadaptDegree - Fine-tuning parameter to alter the degree of
%   dimension-adaptivity [ positive scalar {0.9} ]. A value of 1
%   places complete emphasis on the error estimates, and thus leads
%   to "greedy" dimension-adaptivity. A value of 0 disregards the
%   error estimates, and constructs a conventional sparse grid
%   based on the amount of work involved.
%
% DegreeStrategy - [ {balancing} | depth ] Strategy for the degree
%   of dimensional adaptivity.
%   The 'balancing' strategy balances the number of grid points
%   generated according to the greedy, error estimate-based refine-
%   ment rule compared to the number of points generated by the
%   conventional (regular) sparse grid refinement rule. I.e., a
%   DimadaptDegree value of 0.9 would mean that around 90% of the
%   grid points are generated by the error estimate-based rule, and
%   the remaining points are selected according to the regular
%   rule.
%   The 'depth' strategy makes sure that the maximum level depth
%   reached by the error estimate-based refinement in one dimension
%   does not get too deep compared to the depth reached in the other
%   dimensions. This measurement degree strategy is the one used
%   prior to version 5.1.0 of the toolbox, and is described on page
%   61 of A. Klimke, "Uncertainty modeling using fuzzy arithmetic
%   and sparse grids". This approach is still supported but no
%   the default strategy.
%
% SparseIndices - [ {auto} | on | off ] Manually turn the efficient
%   sparse storage scheme (new feature since version 3.0) of the
%   multi-index arrays on or off. The default switch auto uses the
%   new scheme for the ClenshawCurtis and the Chebyshev grid, and
%   the old (full) storage scheme from version 2.x for the Maximum
%   and the NoBoundary grid (the sparse grid storage scheme is not
%   supported for these two grid types).
%
% DropTol - Drop tolerance [ {auto} | off | 1x2 double vector ]
%   During the sparse grid construction progress, the spvals
%   algorithm may add subgrids with hierarchical surpluses that are
%   all close to zero or of negligible magnitude compared to the
%   surpluses of other sub-grids. In particular, this occurs when
%   additive structure is present in the objective function.
%   To increase the performance of the spinterp algorithm,
%   subgrids where all (absolute) hierarchical surpluses are less
%   than max(relDropTol*(max(fevalRange)-min(fevalRange),absDropTol))
%   can be omitted from the interpolation process by running the
%   sppurge algoritm (see help sppurge).
%   You may specify the absolute and the relative drop tolerance
%   as a vector [ absDropTol, relDropTol ], or turn it off completely
%   (= behavior of version 3.0 and earlier). The switch 'auto' uses
%   the values absDropTol = 0, relDropTol = 100*eps, that is, by
%   default, only a relative drop tolerance is used.
%
% EnableDCT - [ {on} | off ]  Use the DCT when constructing the
%   Chebyshev-Gauss-Lobatto type sparse grid.

% Author : Andreas Klimke, Universitaet Stuttgart
% Version: 1.9
% Date   : November 18, 2007

% Change log:
% V1.0   : Sep 23, 2003
%          Initial release
% V1.1   : Dec 23, 2003
%          Added FunctionArgType handling for more flexibility in
%          defining the objective function.
% V1.2   : Jan 07, 2004
%          Added KeepFunctionValues property to allow the user to
%          access the function values at the sparse grid points.
% V1.3   : Jan 08, 2004
%          Added KeepGrid property to allow the user to access the
%          sparse grid points.
% V1.4   : Apr 22, 2004
%          Added options for dimension-adaptive sparse grids.
% V1.5   : Feb 22, 2005
%          Added handling of sparse sets of indices (useful for
%          high-dimensional interpolation)
% V1.6   : January 13, 2006
%          Added documentation of SparseIndices property
% V1.7   : February 2, 2006
%          Added DropTol property
% V1.8   : March 1, 2006
%          Added DCT feature
% V1.9   : November 18, 2007
%          Added new grid type : Gauss-Patterson
% V2.0   : January 20, 2007
%          Added new parameter for dimensional adaptivity degree
%          strategy.

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

% Note: SPSET is similar in syntax and code the the options
% handling with ODEGET and ODESET of the MATLAB ODE suite by Marc
% Reichelt and Lawrence Shampine.

% Print out possible values of properties.
if (nargin == 0) && (nargout == 0)
  fprintf(['          GridType: [ {Clenshaw-Curtis} | Maximum |' ...
					 ' NoBoundary | Chebyshev | Gauss-Patterson ]\n']);
  fprintf('            RelTol: [ positive scalar {1e-2} ]\n');
  fprintf('            AbsTol: [ positive scalar {1e-6} ]\n');
  fprintf('        Vectorized: [ {off} | on ]\n');
	fprintf('          MinDepth: [ positive integer {2} ]\n');
	fprintf('          MaxDepth: [ positive integer {8} ] (6 for Gauss-Patt.)\n');
  fprintf(' VariablePositions: [ double matrix {[]} ]\n');
  fprintf('   NumberOfOutputs: [ positive integer {1} ]\n');
  fprintf('       PrevResults: [ struct {[]} ]\n');
	fprintf('   FunctionArgType: [ {list} | vector ]\n');
	fprintf('KeepFunctionValues: [ {off} | on ]\n');
	fprintf('          KeepGrid: [ {off} | on ]\n');
	fprintf(' DimensionAdaptive: [ {off} | on ]\n');
	fprintf('         MinPoints: [ positive integer {100} ]\n');
	fprintf('         MaxPoints: [ positive integer {10000} ]\n');
	fprintf('    DimadaptDegree: [ positive scalar in [0,1] {0.9} ]\n');
	fprintf('    DegreeStrategy: [ {balancing} | depth ]\n');
	fprintf('     SparseIndices: [ {auto} | on | off ]\n');
	fprintf(['           DropTol: [ {auto} | off |' ...
					 ' [absTol, relTol] (pos. scalars)]\n']);
	fprintf('         EnableDCT: [ {on} | off ]\n');
  fprintf('\n');
  return;
end

Names = {'GridType', 'RelTol', 'AbsTol', 'Vectorized', 'MinDepth', ...
				 'MaxDepth', 'VariablePositions', 'NumberOfOutputs', ...
				 'PrevResults', 'FunctionArgType', 'KeepFunctionValues', ...
				 'KeepGrid', 'DimensionAdaptive', 'MinPoints', ...
				 'MaxPoints', 'DimadaptDegree', 'DegreeStrategy', ...
				 'SparseIndices', 'DropTol', 'EnableDCT'};

m = length(Names);
names = cell(m,1);
for k = 1:m
	names{k} = lower(Names{k});
end

% Combine all leading options structures o1, o2, ... in odeset(o1,o2,...).
options = [];
for k = 1:m
  options.(Names{k}) = [];
end

k = 1;
while k <= nargin
  arg = varargin{k};
  if ischar(arg)                        % arg is an option name
    break;
  end
  if ~isempty(arg)                      % [] is a valid options argument
    if ~isa(arg,'struct')
      msg = sprintf(['Expected argument %d to be a string property' ...
										 ' name or an options structure\ncreated with' ...
										 ' SPSET.'], k);
      error(msg);
    end
    for l = 1:m
      if any(strcmp(fieldnames(arg),Names{l}))
        val = arg.(Names{l});
      else
        val = [];
      end
      if ~isempty(val)
        options.(Names{l}) = val;
      end
    end
		k = k + 1;
  else
		k = k + 1;
		break;
	end
end

% now process the name-value pairs.
if rem(nargin-k+1,2) ~= 0
  error('Arguments must occur in name-value pairs.');
end

expectval = 0;                      % start expecting a name, not a
                                    % value
while k <= nargin
  arg = varargin{k};

  if ~expectval
    if ~ischar(arg)
      msg = sprintf(['Expected argument %d to be a string property' ...
										 ' name.'], k);
      error(msg);
    end
    lowArg = lower(arg);
    matched = strmatch(lowArg,names);
    if isempty(matched)
      msg = sprintf('Unrecognized property name ''%s''.', arg);
      error(msg);
    elseif length(matched) > 1
      msg = sprintf('Ambiguous property name ''%s'' ', arg);
      msg = [msg '(' Names{matched(1)}];
      for l = matched(2:length(matched))'
        msg = [msg ', ' Names{l}];
      end
      msg = sprintf('%s).', msg);
      error(msg);
    end
    expectval = 1;                      % we expect a value next

  else
    options.(Names{matched}) = arg;
    expectval = 0;

  end
  k = k + 1;
end

if expectval
  msg = sprintf('Expected value for property ''%s''.', arg);
  error(msg);
end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% KERNEL/DISTANCES
%
% Contents of FAIR's Distance Toolbox
%
% This folder provides several distance measures, wrappers and utilities.
%
% general purpose tools:
%   contents        - this file
%   distance        - the specific distance measure used in FAIR
%     initialize:   distance('reset','distance','SSD');
%     usage:        D = distance(Tc,Rc,omega,m);
%   testDistances   - test the files in this folder
%   see also E9_Hands_MLIR_SSD_mbElas
%
% measures and wrappers:
%   MI              - wrapper for MIcc
%   MImex           - Mutual Information based on rhoSplineC (C-implementation)
%   MIspline        - Mutual Information based on rhoSpline  (MATLAB-implementation)
%   rhoSpline       - MATLAB implementation of joint entropy estimator (spline kernel)
%   rhoSplineC      - C based implementation of joint entropy estimator (spline kernel)
%   NCC             - Normalized Cross Correlation
%   NCCmex          - Wrapper for NCCmexC
%   NCCmexC         - C version of NCC
%   NGF             - wrapper for NGFdot
%   NGFdot          - Normalized Gradient Fields, dot product based (THE implementation)
%   NGFcross        - Normalized Gradient Fields, cross product based (to be used with care)
%   NGFmex          - Wrapper for NGFdotMexC
%   NGFdotMexC      - C version of NGFdot
%   SSD             - Sum of Squared Differences
%   SSDmex          - Wrapper for SSDmexC
%   SSDmexC         - C version of SSD
%   parseDistance   - parse distances for derivative checks
%  see also E9_Hands_MLIR_SSD_mbElas
%
%==============================================================================

function debit = contents
if nargout == 0,
    help(mfilename);
    return;
end;

debit = {
  'contents.m'
  'distance.m'
  'testDistances.m'

  'SSD.m'
  'SSDmex.m'
  'SSDmexC.cpp'

  'NCC.m'
  'NCCmex.m'
  'NCCmexC.cpp'

  'NGF.m'
  'NGFcross.m'
  'NGFdot.m'
  'NGFmex.m'
  'NGFdotMexC.cpp'

  'rhoSpline.m'
  'rhoSplineC.cpp'

  'MI.m'
  'MImex.m'
  'MIspline.m'

  'parseDistance.m'
  };
  %==============================================================================


text:
function D = EuDist2(fea_a,fea_b,bSqrt)
% Euclidean Distance matrix
%   D = EuDist(fea_a,fea_b)
%   fea_a:    nSample_a * nFeature
%   fea_b:    nSample_b * nFeature
%   D:      nSample_a * nSample_a
%       or  nSample_a * nSample_b

if ~exist('bSqrt','var')
    bSqrt = 1;
end

if (~exist('fea_b','var')) | isempty(fea_b)
    [nSmp, nFea] = size(fea_a);

    aa = sum(fea_a.*fea_a,2);
    ab = fea_a*fea_a';

    aa = full(aa);
    ab = full(ab);

    if bSqrt
        D = sqrt(repmat(aa, 1, nSmp) + repmat(aa', nSmp, 1) - 2*ab);
        D = real(D);
    else
        D = repmat(aa, 1, nSmp) + repmat(aa', nSmp, 1) - 2*ab;
    end

    D = max(D,D');
    D = D - diag(diag(D));
    D = abs(D);
else
    [nSmp_a, nFea] = size(fea_a);
    [nSmp_b, nFea] = size(fea_b);

    aa = sum(fea_a.*fea_a,2);
    bb = sum(fea_b.*fea_b,2);
    ab = fea_a*fea_b';

    aa = full(aa);
    bb = full(bb);
    ab = full(ab);

    if bSqrt
        D = sqrt(repmat(aa, 1, nSmp_b) + repmat(bb', nSmp_a, 1) - 2*ab);
        D = real(D);
    else
        D = repmat(aa, 1, nSmp_b) + repmat(bb', nSmp_a, 1) - 2*ab;
    end

    D = abs(D);
end

repo: lrslibrary


text:
function [xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk,dxpdalpha] = project_points2(X,om,T,f,c,k,alpha)

%project_points2.m
%
%[xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk] = project_points2(X,om,T,f,c,k,alpha)
%
%Projects a 3D structure onto the image plane.
%
%INPUT: X: 3D structure in the world coordinate frame (3xN matrix for N points)
%       (om,T): Rigid motion parameters between world coordinate frame and camera reference frame
%               om: rotation vector (3x1 vector); T: translation vector (3x1 vector)
%       f: camera focal length in units of horizontal and vertical pixel units (2x1 vector)
%       c: principal point location in pixel units (2x1 vector)
%       k: Distortion coefficients (radial and tangential) (4x1 vector)
%       alpha: Skew coefficient between x and y pixel (alpha = 0 <=> square pixels)
%
%OUTPUT: xp: Projected pixel coordinates (2xN matrix for N points)
%        dxpdom: Derivative of xp with respect to om ((2N)x3 matrix)
%        dxpdT: Derivative of xp with respect to T ((2N)x3 matrix)
%        dxpdf: Derivative of xp with respect to f ((2N)x2 matrix if f is 2x1, or (2N)x1 matrix is f is a scalar)
%        dxpdc: Derivative of xp with respect to c ((2N)x2 matrix)
%        dxpdk: Derivative of xp with respect to k ((2N)x4 matrix)
%
%Definitions:
%Let P be a point in 3D of coordinates X in the world reference frame (stored in the matrix X)
%The coordinate vector of P in the camera reference frame is: Xc = R*X + T
%where R is the rotation matrix corresponding to the rotation vector om: R = rodrigues(om);
%call x, y and z the 3 coordinates of Xc: x = Xc(1); y = Xc(2); z = Xc(3);
%The pinehole projection coordinates of P is [a;b] where a=x/z and b=y/z.
%call r^2 = a^2 + b^2.
%The distorted point coordinates are: xd = [xx;yy] where:
%
%xx = a * (1 + kc(1)*r^2 + kc(2)*r^4 + kc(5)*r^6)      +      2*kc(3)*a*b + kc(4)*(r^2 + 2*a^2);
%yy = b * (1 + kc(1)*r^2 + kc(2)*r^4 + kc(5)*r^6)      +      kc(3)*(r^2 + 2*b^2) + 2*kc(4)*a*b;
%
%The left terms correspond to radial distortion (6th degree), the right terms correspond to tangential distortion
%
%Finally, convertion into pixel coordinates: The final pixel coordinates vector xp=[xxp;yyp] where:
%
%xxp = f(1)*(xx + alpha*yy) + c(1)
%yyp = f(2)*yy + c(2)
%
%
%NOTE: About 90 percent of the code takes care fo computing the Jacobian matrices
%
%
%Important function called within that program:
%
%rodrigues.m: Computes the rotation matrix corresponding to a rotation vector
%
%rigid_motion.m: Computes the rigid motion transformation of a given structure

if nargin < 7,
   alpha = 0;
   if nargin < 6,
      k = zeros(5,1);
      if nargin < 5,
         c = zeros(2,1);
         if nargin < 4,
            f = ones(2,1);
            if nargin < 3,
               T = zeros(3,1);
               if nargin < 2,
                  om = zeros(3,1);
                  if nargin < 1,
                     error('Need at least a 3D structure to project (in project_points.m)');
                     return;
                  end;
               end;
            end;
         end;
      end;
   end;
end;

[m,n] = size(X);

[Y,dYdom,dYdT] = rigid_motion(X,om,T);

inv_Z = 1./Y(3,:);

x = (Y(1:2,:) .* (ones(2,1) * inv_Z)) ;

bb = (-x(1,:) .* inv_Z)'*ones(1,3);
cc = (-x(2,:) .* inv_Z)'*ones(1,3);

dxdom = zeros(2*n,3);
dxdom(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(1:3:end,:) + bb .* dYdom(3:3:end,:);
dxdom(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(2:3:end,:) + cc .* dYdom(3:3:end,:);

dxdT = zeros(2*n,3);
dxdT(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(1:3:end,:) + bb .* dYdT(3:3:end,:);
dxdT(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(2:3:end,:) + cc .* dYdT(3:3:end,:);

% Add distortion:

r2 = x(1,:).^2 + x(2,:).^2;

dr2dom = 2*((x(1,:)')*ones(1,3)) .* dxdom(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdom(2:2:end,:);
dr2dT = 2*((x(1,:)')*ones(1,3)) .* dxdT(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdT(2:2:end,:);

r4 = r2.^2;

dr4dom = 2*((r2')*ones(1,3)) .* dr2dom;
dr4dT = 2*((r2')*ones(1,3)) .* dr2dT;

r6 = r2.^3;

dr6dom = 3*((r2'.^2)*ones(1,3)) .* dr2dom;
dr6dT = 3*((r2'.^2)*ones(1,3)) .* dr2dT;

% Radial distortion:

cdist = 1 + k(1) * r2 + k(2) * r4 + k(5) * r6;

dcdistdom = k(1) * dr2dom + k(2) * dr4dom + k(5) * dr6dom;
dcdistdT = k(1) * dr2dT + k(2) * dr4dT + k(5) * dr6dT;
dcdistdk = [ r2' r4' zeros(n,2) r6'];

xd1 = x .* (ones(2,1)*cdist);

dxd1dom = zeros(2*n,3);
dxd1dom(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdom;
dxd1dom(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdom;
coeff = (reshape([cdist;cdist],2*n,1)*ones(1,3));
dxd1dom = dxd1dom + coeff.* dxdom;

dxd1dT = zeros(2*n,3);
dxd1dT(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdT;
dxd1dT(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdT;
dxd1dT = dxd1dT + coeff.* dxdT;

dxd1dk = zeros(2*n,5);
dxd1dk(1:2:end,:) = (x(1,:)'*ones(1,5)) .* dcdistdk;
dxd1dk(2:2:end,:) = (x(2,:)'*ones(1,5)) .* dcdistdk;

% tangential distortion:

a1 = 2.*x(1,:).*x(2,:);
a2 = r2 + 2*x(1,:).^2;
a3 = r2 + 2*x(2,:).^2;

delta_x = [k(3)*a1 + k(4)*a2 ;
   k(3) * a3 + k(4)*a1];

%ddelta_xdx = zeros(2*n,2*n);
aa = (2*k(3)*x(2,:)+6*k(4)*x(1,:))'*ones(1,3);
bb = (2*k(3)*x(1,:)+2*k(4)*x(2,:))'*ones(1,3);
cc = (6*k(3)*x(2,:)+2*k(4)*x(1,:))'*ones(1,3);

ddelta_xdom = zeros(2*n,3);
ddelta_xdom(1:2:end,:) = aa .* dxdom(1:2:end,:) + bb .* dxdom(2:2:end,:);
ddelta_xdom(2:2:end,:) = bb .* dxdom(1:2:end,:) + cc .* dxdom(2:2:end,:);

ddelta_xdT = zeros(2*n,3);
ddelta_xdT(1:2:end,:) = aa .* dxdT(1:2:end,:) + bb .* dxdT(2:2:end,:);
ddelta_xdT(2:2:end,:) = bb .* dxdT(1:2:end,:) + cc .* dxdT(2:2:end,:);

ddelta_xdk = zeros(2*n,5);
ddelta_xdk(1:2:end,3) = a1';
ddelta_xdk(1:2:end,4) = a2';
ddelta_xdk(2:2:end,3) = a3';
ddelta_xdk(2:2:end,4) = a1';

xd2 = xd1 + delta_x;

dxd2dom = dxd1dom + ddelta_xdom ;
dxd2dT = dxd1dT + ddelta_xdT;
dxd2dk = dxd1dk + ddelta_xdk ;

% Add Skew:

xd3 = [xd2(1,:) + alpha*xd2(2,:);xd2(2,:)];

% Compute: dxd3dom, dxd3dT, dxd3dk, dxd3dalpha

dxd3dom = zeros(2*n,3);
dxd3dom(1:2:2*n,:) = dxd2dom(1:2:2*n,:) + alpha*dxd2dom(2:2:2*n,:);
dxd3dom(2:2:2*n,:) = dxd2dom(2:2:2*n,:);
dxd3dT = zeros(2*n,3);
dxd3dT(1:2:2*n,:) = dxd2dT(1:2:2*n,:) + alpha*dxd2dT(2:2:2*n,:);
dxd3dT(2:2:2*n,:) = dxd2dT(2:2:2*n,:);
dxd3dk = zeros(2*n,5);
dxd3dk(1:2:2*n,:) = dxd2dk(1:2:2*n,:) + alpha*dxd2dk(2:2:2*n,:);
dxd3dk(2:2:2*n,:) = dxd2dk(2:2:2*n,:);
dxd3dalpha = zeros(2*n,1);
dxd3dalpha(1:2:2*n,:) = xd2(2,:)';

% Pixel coordinates:
if length(f)>1,
    xp = xd3 .* (f * ones(1,n))  +  c*ones(1,n);
    coeff = reshape(f*ones(1,n),2*n,1);
    dxpdom = (coeff*ones(1,3)) .* dxd3dom;
    dxpdT = (coeff*ones(1,3)) .* dxd3dT;
    dxpdk = (coeff*ones(1,5)) .* dxd3dk;
    dxpdalpha = (coeff) .* dxd3dalpha;
    dxpdf = zeros(2*n,2);
    dxpdf(1:2:end,1) = xd3(1,:)';
    dxpdf(2:2:end,2) = xd3(2,:)';
else
    xp = f * xd3 + c*ones(1,n);
    dxpdom = f  * dxd3dom;
    dxpdT = f * dxd3dT;
    dxpdk = f  * dxd3dk;
    dxpdalpha = f .* dxd3dalpha;
    dxpdf = xd3(:);
end;

dxpdc = zeros(2*n,2);
dxpdc(1:2:end,1) = ones(n,1);
dxpdc(2:2:end,2) = ones(n,1);

return;

% Test of the Jacobians:

n = 10;

X = 10*randn(3,n);
om = randn(3,1);
T = [10*randn(2,1);40];
f = 1000*rand(2,1);
c = 1000*randn(2,1);
k = 0.5*randn(5,1);
alpha = 0.01*randn(1,1);

[x,dxdom,dxdT,dxdf,dxdc,dxdk,dxdalpha] = project_points2(X,om,T,f,c,k,alpha);

% Test on om: OK

dom = 0.000000001 * norm(om)*randn(3,1);
om2 = om + dom;

[x2] = project_points2(X,om2,T,f,c,k,alpha);

x_pred = x + reshape(dxdom * dom,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on T: OK!!

dT = 0.0001 * norm(T)*randn(3,1);
T2 = T + dT;

[x2] = project_points2(X,om,T2,f,c,k,alpha);

x_pred = x + reshape(dxdT * dT,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on f: OK!!

df = 0.001 * norm(f)*randn(2,1);
f2 = f + df;

[x2] = project_points2(X,om,T,f2,c,k,alpha);

x_pred = x + reshape(dxdf * df,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on c: OK!!

dc = 0.01 * norm(c)*randn(2,1);
c2 = c + dc;

[x2] = project_points2(X,om,T,f,c2,k,alpha);

x_pred = x + reshape(dxdc * dc,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on k: OK!!

dk = 0.001 * norm(k)*randn(5,1);
k2 = k + dk;

[x2] = project_points2(X,om,T,f,c,k2,alpha);

x_pred = x + reshape(dxdk * dk,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on alpha: OK!!

dalpha = 0.001 * norm(k)*randn(1,1);
alpha2 = alpha + dalpha;

[x2] = project_points2(X,om,T,f,c,k,alpha2);

x_pred = x + reshape(dxdalpha * dalpha,2,n);

norm(x2-x)/norm(x2 - x_pred)

repo: ekfmonoslam


text:
function hRD = plugin_rd_rxtx
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevices.fmcomms5.common.plugin_rd('ZC702', 'Rx & Tx');
AnalogDevices.fmcomms5.zc702.add_io(hRD,'Rx & Tx');

author: analogdevicesinc


text:
function number = comp_enum ( n, k )

%*****************************************************************************80
%
%% COMP_ENUM returns the number of compositions of the integer N into K parts.
%
%  Discussion:
%
%    A composition of the integer N into K parts is an ordered sequence
%    of K nonnegative integers which sum to N.  The compositions (1,2,1)
%    and (1,1,2) are considered to be distinct.
%
%    The 28 compositions of 6 into three parts are:
%
%      6 0 0,  5 1 0,  5 0 1,  4 2 0,  4 1 1,  4 0 2,
%      3 3 0,  3 2 1,  3 1 2,  3 0 3,  2 4 0,  2 3 1,
%      2 2 2,  2 1 3,  2 0 4,  1 5 0,  1 4 1,  1 3 2,
%      1 2 3,  1 1 4,  1 0 5,  0 6 0,  0 5 1,  0 4 2,
%      0 3 3,  0 2 4,  0 1 5,  0 0 6.
%
%    The formula for the number of compositions of N into K parts is
%
%      Number = ( N + K - 1 )! / ( N! * ( K - 1 )! )
%
%    Describe the composition using N '1's and K-1 dividing lines '|'.
%    The number of distinct permutations of these symbols is the number
%    of compositions.  This is equal to the number of permutations of
%    N+K-1 things, with N identical of one kind and K-1 identical of another.
%
%    Thus, for the above example, we have:
%
%      Number = ( 6 + 3 - 1 )! / ( 6! * (3-1)! ) = 28
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms for Computers and Calculators,
%    Second Edition,
%    Academic Press, 1978,
%    ISBN: 0-12-519260-6,
%    LC: QA164.N54.
%
%  Parameters:
%
%    Input, integer N, the integer whose compositions are desired.
%
%    Input, integer K, the number of parts in the composition.
%
%    Output, integer NUMBER, the number of compositions of N
%    into K parts.
%
  number = i4_choose ( n + k - 1, n );

  return
end


text:
function dropout_forward()
    global config mem;
    curr_layer_idx = config.misc.current_layer - 1;
    if config.misc.training;
        dm = config.NEW_MEM(single(uint8(rand(size(mem.activations{curr_layer_idx})))));
        mem.activations{curr_layer_idx} = mem.activations{curr_layer_idx} .* dm;
    else
        mem.activations{curr_layer_idx} = mem.activations{curr_layer_idx} ./ 2;
    end
end

author: jimmy-ren


text:
function [Dictionary, data, coefs] = gererateSyntheticDictionaryAndData(N, L, dim, K, SNRdB)

randn('state',sum(100*clock));
rand('state',sum(100*clock));

Dictionary = randn(dim,K);
Dictionary = Dictionary*diag(1./sqrt(sum(Dictionary.*Dictionary)));

[data,coefs] = CreateDataFromDictionarySimple(Dictionary, N, L);

if (SNRdB==0) | (SNRdB == 80)
    return
else
    noise = randn(size(data));
    actualNoise = calcNoiseFromSNR(SNRdB,data, noise);
    SNR = calcSNR(data, data+actualNoise);
    data =  data + actualNoise*SNR/SNRdB;
end

function [D,xOrig] = CreateDataFromDictionarySimple(dictionary, numElements, numCoef)
maxRangeOfCoef = 1;
resolution = 0.0001;

xOrig = zeros(size(dictionary,2),numElements);
%vecOfValues = -1*maxRangeOfCoef:resolution:maxRangeOfCoef;
%coefs = randsrc(numCoef,numElements,vecOfValues);
coefs = randn(numCoef,numElements)*maxRangeOfCoef;
xOrig(1:numCoef,:) = coefs;
for i=1:size(xOrig,2)
    xOrig(:,i) = xOrig(randperm(size(xOrig,1)),i);
end
%dictionaryElementIndices = randsrc(numCoef*numElements,1,[1:size(dictionary,2)])   ;
%matrixOfIndices = repmat([1:numElements],numCoef,1);
%xOrig(sub2ind(size(xOrig),dictionaryElementIndices,matrixOfIndices(:))) = coefs;
D = dictionary*xOrig;

function  actualNoise = calcNoiseFromSNR(TargerSNR, signal, randomNoise)
signal = signal(:);
randomNoiseRow = randomNoise(:);
signal_2 = sum(signal.^2);
ActualNoise_2 = signal_2/(10^(TargerSNR/10));
noise_2 = sum(randomNoiseRow.^2);
ratio = ActualNoise_2./noise_2;
actualNoise = randomNoiseRow.*repmat(sqrt(ratio),size(randomNoiseRow,1),1);
actualNoise = reshape(actualNoise,size(randomNoise));

function SNR = calcSNR(origSignal, noisySignal)
errorSignal = origSignal-noisySignal;
signal_2 = sum(origSignal.^2);
noise_2 = sum(errorSignal.^2);

SNRValues = 10*log10(signal_2./noise_2);
SNR = mean(SNRValues);

author: lbasek

repo: image-denoising-benchmark


text:
function [d,a] =  training(a,d)

  for i=1:length(a.child)
    [d a.child{i}]=train(a.child{i},d);
    if i<length(a.child) & (isa(a.child{i},'group') | isa(a.child{i},'param'))
      a.child{i}.group='one_for_each';
    end;
   end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore