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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym acsc (@var{x})
%% Symbolic acsc function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = acsc (x)
%%   @result{} y = (sym) acsc(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = acsc(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('acsc', x);
end

%!error acsc (sym(1), 2)
%!assert (isequaln (acsc (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = acsc(x);
%! f2 = acsc(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = acsc(A);
%! f2 = acsc(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = acsc (d);
%! f = acsc (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)


text:
function area = triangleArea(pt1, pt2, pt3)
%TRIANGLEAREA Signed area of a triangle.
%
%   AREA = triangleArea(P1, P2, P3)
%   Computes area of the triangle whose vertices are given by P1, P2 and
%   P3. Each vertex is a 1-by-2 row vector.
%
%   AREA = triangleArea(PTS)
%   Concatenates vertex coordinates in a 3-by-2 array. Each row of the
%   array contains coordinates of one vertex.
%
%
%   Example
%   % Compute area of a Counter-Clockwise (CCW) oriented triangle
%     triangleArea([10 10], [30 10], [10 40])
%     ans =
%         300
%
%   % Compute area of a Clockwise (CW) oriented triangle
%     triangleArea([10 40], [30 10], [10 10])
%     ans =
%         -300
%
%   See also
%   polygonArea, triangleArea3d

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-08-23, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% if data is given as one array, split vertices
if nargin == 1
    pt2 = pt1(2,:);
    pt3 = pt1(3,:);
    pt1 = pt1(1,:);
end

% compute individual vectors
v12 = bsxfun(@minus, pt2, pt1);
v13 = bsxfun(@minus, pt3, pt1);

% compute area from cross product
area = (v13(:,2) .* v12(:,1) - v13(:,1) .* v12(:,2)) / 2;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function res = bf_inverse_champagne(BF, S)
% Computes Champagne filters
% See Owen et al. Performance evaluation of the Champagne source
% reconstruction algorithm on simulated and real M/EEG data. Neuroimage. 2012 Mar;60(1):305-23
% Code contributed by Sri Nagarajan
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: bf_inverse_champagne.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
if nargin == 0
    nem = cfg_entry;
    nem.tag = 'nem';
    nem.name = 'Number of EM iterations';
    nem.strtype = 'n';
    nem.num = [1 1];
    nem.val = {100};

    vcs         = cfg_menu;
    vcs.tag     = 'vcs';
    vcs.name    = 'Voxel covariance structure';
    vcs.labels  = {
        'scalar'
        'diagonal'
        'general'
        }';
    vcs.values  = {0, 1, 2};
    vcs.val = {2};

    nupd         = cfg_menu;
    nupd.tag     = 'nupd';
    nupd.name    = 'Noise covariance';
    nupd.labels  = {
        'use provided'
        'learn scalar'
        'learn diagonal'
        }';
    nupd.values  = {0, 1, 2};
    nupd.val = {0};

    champagne      = cfg_branch;
    champagne.tag  = 'champagne';
    champagne.name = 'Champagne';
    champagne.val  = {nem, vcs, nupd};
    res = champagne;

    return
elseif nargin < 2
    error('Two input arguments are required');
end

C     =  BF.features.(S.modality).C;
Y     =  BF.features.(S.modality).Y;

L = S.L;

W = cell(size(L));
nvert = numel(W);
nd = size(L{1}, 2);

LL = zeros(size(L{1}, 1), nvert*nd);
ind = 1;

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', nvert, ['Preparing ' S.modality ' leadfields']); drawnow;
if nvert > 100, Ibar = floor(linspace(1, nvert,100));
else Ibar = 1:nvert; end

for i = 1:nvert

    for j = 1:nd
        lf = L{i}(:, j);
        LL(:, ind) = lf./norm(lf);
        ind = ind+1;
    end
    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

Fgraph  = spm_figure('GetWin', 'Champagne'); figure(Fgraph); clf

%****************************************************************************************
[gamma,x,w,sigu,like]=champagne_aug2015(Y, LL, C, S.nem, nd, S.vcs, S.nupd, [], 0, Fgraph);
%****************************************************************************************

spm_progress_bar('Init', nvert, ['Preparing ' S.modality ' filters']); drawnow;
if nvert > 100, Ibar = floor(linspace(1, nvert,100));
else Ibar = 1:nvert; end

for i = 1:nvert
    W{i} = w((i-1)*nd+(1:3), :);

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

spm('Pointer', 'Arrow');

res.W = W;

author: spm


text:
%{
 * Copyright (C) 2013-2020, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function cost = costIoU(theta_x, theta_y, theta_z, T, X, Y, intrinsic)
    H = eye(4);
    H(1:3,1:3) = rotx(theta_x) * roty(theta_y) * rotz(theta_z);
    H(1:3,4) = T';
    P = intrinsic * [eye(3) zeros(3,1)] * H;
    L_X_transformed = P * X; % transformed points in LiDAR frame
    C_X_transformed = L_X_transformed ./ L_X_transformed(3,:);
    cost = 0;
    for i=1:size(C_X_transformed,2)/4
        k = 4*(i-1) + 1;
        vertices_X = [C_X_transformed(1,k) C_X_transformed(1,k+1) C_X_transformed(1,k+3) C_X_transformed(1,k+2);
                      C_X_transformed(2,k) C_X_transformed(2,k+1) C_X_transformed(2,k+3) C_X_transformed(2,k+2)];
        vertices_Y = [Y(1,k) Y(1,k+1) Y(1,k+3) Y(1,k+2);
                      Y(2,k) Y(2,k+1) Y(2,k+3) Y(2,k+2)];
%         vertices_X = C_X_transformed(1:2,k:k+3);
%         vertices_Y = Y(1:2,k:k+3);

%         conv_X = convhull(vertices_X(1,:)', vertices_X(2,:)');
%         conv_Y = convhull(vertices_Y(1,:)', vertices_Y(2,:)');
%         [vertices_X(1,:), vertices_X(2,:)] = poly2cw(vertices_X(1,conv_X(1:4)), vertices_X(2,conv_X(1:4)));
%         [vertices_Y(1,:), vertices_Y(2,:)] = poly2cw(vertices_Y(1,conv_Y(1:4)), vertices_Y(2,conv_Y(1:4)));
%         pgon_X =  polyshape(vertices_X(1,conv_X), vertices_X(2,conv_X));
%         pgon_Y =  polyshape(vertices_Y(1,conv_Y), vertices_Y(2,conv_Y));

        [vertices_X(1,:), vertices_X(2,:)] = poly2cw(vertices_X(1, :), vertices_X(2, :));
        [vertices_Y(1,:), vertices_Y(2,:)] = poly2cw(vertices_Y(1, :), vertices_Y(2, :));

        pgon_X =  polyshape(vertices_X(1, :), vertices_X(2, :));
        pgon_Y =  polyshape(vertices_Y(1, :), vertices_Y(2, :));
        if all(issimplified([pgon_X, pgon_Y]))
            [I_polyout,~,~] = intersect(pgon_X, pgon_Y);

            [U_polyout,~,~] = union(pgon_X ,pgon_Y);
            I_area = max(area(I_polyout), 1e-5);
            U_area = area(U_polyout);
            cost = cost + I_area/U_area;
        else
            cost = -1e3;
        end
    end
    cost = -cost;
end

author: UMich-BipedLab


text:
% ObserveEvidence Modify a vector of factors given some evidence.
%   F = ObserveEvidence(F, E) sets all entries in the vector of factors, F,
%   that are not consistent with the evidence, E, to zero. F is a vector of
%   factors, each a data structure with the following fields:
%     .var    Vector of variables in the factor, e.g. [1 2 3]
%     .card   Vector of cardinalities corresponding to .var, e.g. [2 2 2]
%     .val    Value table of size prod(.card)
%   E is an N-by-2 matrix, where each row consists of a variable/value pair.
%     Variables are in the first column and values are in the second column.

function F = ObserveEvidence(F, E)

% Iterate through all evidence

for i = 1:size(E, 1),
    v = E(i, 1); % variable
    x = E(i, 2); % value

    % Check validity of evidence
    if (x == 0),
        warning(['Evidence not set for variable ', int2str(v)]);
        continue;
    end;

    for j = 1:length(F),
		  % Does factor contain variable?
        indx = find(F(j).var == v);

        if (~isempty(indx)),

		  	   % Check validity of evidence
            if (x > F(j).card(indx) || x < 0 ),
                error(['Invalid evidence, X_', int2str(v), ' = ', int2str(x)]);
            end;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % YOUR CODE HERE
            % Adjust the factor F(j) to account for observed evidence
            % Hint: You might find it helpful to use IndexToAssignment
            %       and SetValueOfAssignment
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           for i = 1:length(F(j).val)
		  assign = IndexToAssignment(i,F(j).card);
		  if assign(indx)!=x
			  F(j) = SetValueOfAssignment(F(j),assign,0);
		  end
	  end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

				% Check validity of evidence / resulting factor
            if (all(F(j).val == 0)),
                warning(['Factor ', int2str(j), ' makes variable assignment impossible']);
            end;

        end;
    end;
end;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
%% "Spirograph" Code
%  Jonathan Jamieson - September 2013
%  unigamer@gmail.com
%  www.jonathanjamieson.com

% This is the main script. Run it in Matlab and then mess around with
% the "User Configurable Parameters". By default arm at the bottom
% is arm "a" and the left arm is arm "b". The "c" coordinates denote
% the position of the rotating platform.

%% User Configurable Parameters -------------------------------------------
% a - arm a parameters
% b - arm b parameters
% c - parameters of the rotating disc where the paper would be

% Time-span
t=1:0.01:900;

% The positions of the local coordinate systems in the global one
xa = 18;
ya = 2.5;

xb = 6.5;
yb = 12;

xc = 18;
yc = 12;

% The positions of where the arms are mounted to the rotating discs
A_x_local =-1.8;
A_y_local = 0.2;

B_x_local = 3;
B_y_local = 1;

% The angular speeds of the rotating discs
A_omega = -3*5;
B_omega = 1*5;
C_omega = 0.22*5;

% The lengths of the arms
La = 8;
Lb = 13;

% The distance between where the two rods join and the pen
extension = 1;

% Radii of the discs
ra = 2; % The "r" parameters are only
rb = 4; % used for the purpose of
rc = 6; % visualisation, not calculation

% End of User Configurable Parameters -------------------------------------

%% Graphical Stuff
scrsz = get(0,'ScreenSize');  % Get screen size
figure('Name','The drawing','Position',[scrsz(1)+scrsz(3)/4 scrsz(2)+scrsz(4)/6 scrsz(3)/2 scrsz(4)/1.5]);
hold on
title('The Plot');
axis equal

field_size = 30;
xlim(gca, [0 field_size]);
ylim(gca, [0 field_size]);

L1_handle = line([0 0],[0 0],'Color','red','LineWidth',2);
L2_handle = line([0 0],[0 0],'Color','red','LineWidth',2);
P_handle = line([0 0],[0 0],'Color','red','LineWidth',2);

PlotCircle(xa,ya,ra);
PlotCircle(xb,yb,rb);
PlotCircle(xc,yc,rc);

pen_log = zeros(numel(t),2); % Pre-allocate the matrix, this speeds stuff up and is good practice.

%% Main Loop
% This loop is executed for each time-step

try
  for loop_counter=1:1:numel(t)

    %% Calculate the position of the pen at the current time-step
    current_time=t(loop_counter);
    [A_x_global,A_y_global] = Local2Global(A_x_local,A_y_local,xa,ya,current_time*A_omega);
    [B_x_global,B_y_global] = Local2Global(B_x_local,B_y_local,xb,yb,current_time*B_omega);
    [C_x_global,C_y_global] = Calculate_Point_C(A_x_global,B_x_global,A_y_global,B_y_global,La,Lb); % This is the intersection of two circles
    [P_x_global,P_y_global ] = CalculatePenPosition(A_x_global,A_y_global,C_x_global,C_y_global,extension);
    [C_x_local,C_y_local] = Global2Local(P_x_global,P_y_global,xc,yc,current_time*C_omega);

    pen_log((loop_counter),:) = [C_x_local,C_y_local]; % Update the pen_log with the new pen position

    %% Create a temporary matrix with the global position of each pen point generated so far
    % at the current time. Since the local coordinate system is rotating a new set of temporary coordinates need to be generated each time step.
    pen_log_temp = zeros(loop_counter,2);

    for loop_plot=1:1:loop_counter;

        [temp_x temp_y] = Local2Global(pen_log(loop_plot,1),pen_log(loop_plot,2),xc,yc,current_time*C_omega);

        pen_log_temp(loop_plot,:) = [temp_x temp_y];
    end

    if loop_counter ~= 1; % Don't try and delete on the first iteration because nothing will have been plotted yet
            delete(h_plot);
    end

    %% Redraw everything for the current time-step
    h_plot = plot(pen_log_temp(:,1),pen_log_temp(:,2),'Color','magenta');

    set(L1_handle,'XData',[A_x_global  C_x_global],'YData',[A_y_global C_y_global]);
    set(L2_handle,'XData',[B_x_global  C_x_global],'YData',[B_y_global C_y_global]);
    set(P_handle,'XData',[C_x_global  P_x_global],'YData',[C_y_global P_y_global]);

    refresh
    pause(0.000001);

  end

    %% Check figure hasn't been closed
    catch E
    if ~strcmp(E.identifier, 'MATLAB:class:InvalidHandle')
        pen_log = pen_log(any(pen_log,2),:); % Delete unused pen_log entries
        x = pen_log(:,1); % Make a vector from the first column
        y = pen_log(:,2); % Make a vector from the second column
        clearvars -except x y % Delete all variables except the coordinates generated in the drawing process. Uncomment this if you wish to debug something.
        break

    end
end

% comet(x,y) % Uncomment this function if you want an animated plot at the end of the drawing
% plot (x,y) % Uncomment this function if you want a plot at the end of the drawing

author: Sable

repo: mcbench-benchmarks


text:
function lab = rgb2lab(rgb_im)

    if ~isa(rgb_im,'uint8'),
        rgb_im = im2uint8(rgb_im);
    end

    cform = makecform('srgb2lab');
    lab = applycform(rgb_im,cform);

end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
function plotci3()
    %  plot a circle containing ni events
    %  around each grid point
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    figure
    orient landscape
    axes('position',[ 0.1 0.1 0.8 0.8])
    set(gca,'NextPlot','add')
    axis([ min(gx) max(gx) min(gy) max(gy)])
    axis image
    xlabel('Distance in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('Depth in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)

    if exist('maex', 'var')
        pl = plot(maex,-maey,'xw');
        set(pl,'MarkerSize',10,'LineWidth',2)
    end
    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','TickDir','out')
    h1 = gca;

    set(gca,'NextPlot','add');
    pl = plot(newgri(:,1),newgri(:,2),'+k','Markersize',3)
    st = 3;
    [X,Y] = meshgrid(gx,gy);
    [m,n]= size(r);
    set(gca,'NextPlot','add')
    x = -pi-0.1:0.1:pi;
    for i = 1:st:m
        for k = 1:st:n
            if r(i,k) <= ZG.tresh_km
                plot(X(i,k)+r(i,k)*1*sin(x),Y(i,k)+r(i,k)*1*cos(x),'color',[0.5 0.5 0.5])
                %   plot(X(i,k)+4.0*sin(x),Y(i,k)+4.0*cos(x),'r')
                set(gca,'NextPlot','add')
                plovo = plot(X(i,k),Y(i,k),'^k','Markersize',5)
                set(plovo,'LineWidth',1,'MarkerSize',5,...
                    'MarkerFaceColor','w','MarkerEdgeColor','k');

            end
        end
    end

    set(gca,'Color',color_bg)

end

author: CelsoReyes


text:
% This is a simplest implementation of the proposed RIFT algorithm. In this implementation,...
% rotation invariance part and corner point detection are not included.

clc;clear;close all;
warning('off')

addpath sar-optical   % type of multi-modal data

str1='pair1.jpg';   % image pair
str2='pair2.jpg';
im1 = im2uint8(imread(str1));
im2 = im2uint8(imread(str2));

if size(im1,3)==1
    temp=im1;
    im1(:,:,1)=temp;
    im1(:,:,2)=temp;
    im1(:,:,3)=temp;
end

if size(im2,3)==1
    temp=im2;
    im2(:,:,1)=temp;
    im2(:,:,2)=temp;
    im2(:,:,3)=temp;
end

disp('RIFT feature detection and description')
% RIFT feature detection and description
[des_m1,des_m2] = RIFT_no_rotation_invariance(im1,im2,4,6,96);

disp('nearest matching')
% nearest matching
[indexPairs,matchmetric] = matchFeatures(des_m1.des,des_m2.des,'MaxRatio',1,'MatchThreshold', 100);
matchedPoints1 = des_m1.kps(indexPairs(:, 1), :);
matchedPoints2 = des_m2.kps(indexPairs(:, 2), :);
[matchedPoints2,IA]=unique(matchedPoints2,'rows');
matchedPoints1=matchedPoints1(IA,:);

disp('outlier removal')
%outlier removal
H=FSC(matchedPoints1,matchedPoints2,'affine',2);
Y_=H*[matchedPoints1';ones(1,size(matchedPoints1,1))];
Y_(1,:)=Y_(1,:)./Y_(3,:);
Y_(2,:)=Y_(2,:)./Y_(3,:);
E=sqrt(sum((Y_(1:2,:)-matchedPoints2').^2));
inliersIndex=E<3;
cleanedPoints1 = matchedPoints1(inliersIndex, :);
cleanedPoints2 = matchedPoints2(inliersIndex, :);

disp('Show matches')
% Show results
figure; showMatchedFeatures(im1, im2, cleanedPoints1, cleanedPoints2, 'montage');

disp('registration result')
% registration
image_fusion(im2,im1,double(H));

author: LJY-RS

repo: RIFT-multimodal-image-matching


text:
% SOSDEMO7 --- Chebyshev polynomials
% Section 3.7 of SOSTOOLS User's Manual

clear; echo on;

ndeg = 8;   % Degree of Chebyshev polynomial

syms x gam;

% =============================================
% First, initialize the sum of squares program
prog = sosprogram([x],[gam]);

% Create the polynomial P
Z = monomials(x,[0:ndeg-1]);
[prog,P1] = sospolyvar(prog,Z);
P = P1 + gam * x^ndeg;           % The leading coeff of P is gam

% Imposing the inequalities
prog = sosineq(prog, 1 - P, [-1, 1]);
prog = sosineq(prog, 1 + P, [-1, 1]);

% And setting objective
prog = sossetobj(prog, -gam);

% Then solve the program
prog = sossolve(prog);

% =============================================
% Finally, get solution
SOLV = sosgetsol(prog, P)
GAM = sosgetsol(prog, gam)

echo off

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [nodes] = idtf_model_nodes(n_meshes, n_lines, n_pointsets)
% File:      idtf_model_nodes.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.06.10 - 2012.06.24
% Language:  MATLAB R2012a
% Purpose:   create IDTF model nodes
% Copyright: Ioannis Filippidis, 2012-

% depends
%   verbatim

%% mesh nodes
mesh_node = mesh_node_str;
temp_nodes = cell(1, n_meshes);
for i=1:n_meshes
    cur_mesh_node = sprintf(mesh_node, i, i);
    temp_nodes{1, i} = cur_mesh_node;
end
mesh_nodes = [temp_nodes{:} ];
%disp(mesh_nodes)

%% line nodes
line_node = line_node_str;
temp_nodes = cell(1, n_lines);
for i=1:n_lines
    cur_line_node = sprintf(line_node, i, i);
    temp_nodes{1, i} = cur_line_node;
end
line_nodes = [temp_nodes{:} ];
%disp(line_nodes)

%% pointset nodes
point_node = point_node_str;
temp_nodes = cell(1, n_pointsets);
for i=1:n_pointsets
    cur_point_node = sprintf(point_node, i, i);
    temp_nodes{1, i} = cur_point_node;
end
point_nodes = [temp_nodes{:} ];
%disp(point_nodes)

%% output
nodes = [mesh_nodes, line_nodes, point_nodes];

function [str] = mesh_node_str
%mesh_id_number, mesh_id_number
str = verbatim;
%{

NODE "MODEL" {
     NODE_NAME "Mesh%d"
     PARENT_LIST {
          PARENT_COUNT 1
          PARENT 0 {
               PARENT_NAME "<NULL>"
               PARENT_TM {
                    1.000000 0.000000 0.000000 0.000000
                    0.000000 1.000000 0.000000 0.000000
                    0.000000 0.000000 1.000000 0.000000
                    0.000000 0.000000 0.000000 1.000000
               }
          }
     }
     RESOURCE_NAME "MyMesh%d"
     MODEL_VISIBILITY "BOTH"
}

%}

function [str] = line_node_str
% lineset_id_number, lineset_id_number
str = verbatim;
%{

NODE "MODEL" {
     NODE_NAME "Line%d"
     PARENT_LIST {
          PARENT_COUNT 1
          PARENT 0 {
               PARENT_NAME "<NULL>"
               PARENT_TM {
                    1.000000 0.000000 0.000000 0.000000
                    0.000000 1.000000 0.000000 0.000000
                    0.000000 0.000000 1.000000 0.000000
                    0.000000 0.000000 0.000000 1.000000
               }
          }
     }
     RESOURCE_NAME "MyLine%d"
}

%}

function [str] = point_node_str
% pointset_id_number, pointset_id_number
str = verbatim;
%{

NODE "MODEL" {
     NODE_NAME "Points%d"
     PARENT_LIST {
          PARENT_COUNT 1
          PARENT 0 {
               PARENT_NAME "<NULL>"
               PARENT_TM {
                    1.000000 0.000000 0.000000 0.000000
                    0.000000 1.000000 0.000000 0.000000
                    0.000000 0.000000 1.000000 0.000000
                    0.000000 0.000000 0.000000 1.000000
               }
          }
     }
     RESOURCE_NAME "MyPoints%d"
}

%}


text:
function [KKTError,costL,timeElapsed] = NMPC_KKTError(x0,p,rho,lambda,mu,u,x)

    timerStart = Timer();
    % Global variables
    global ParNMPCGlobalVariable
    lambdaDim                  = ParNMPCGlobalVariable.dim.x;
    xDim                       = ParNMPCGlobalVariable.dim.x;

    % init
    KKTError.stateEquation   = 0;
    KKTError.C               = 0;
    KKTError.Hu              = 0;
    KKTError.costateEquation = 0;

    [~,sizeSeg,DoP]  = size(x);
    if coder.target('MATLAB') % Normal excution
        % in serial
        numThreads = 0;
    else % Code generation
        numThreads = DoP;
    end

    % Local
    KKTxEquation      = zeros(sizeSeg,DoP);
    KKTC              = zeros(sizeSeg,DoP);
    KKTHu             = zeros(sizeSeg,DoP);
    KKTlambdaEquation = zeros(sizeSeg,DoP);
    L                 = zeros(sizeSeg,DoP);

    % Coupling variable for each segment
    lambdaNext      = zeros(lambdaDim,sizeSeg,DoP);
    xPrev           = zeros(xDim,sizeSeg,DoP);
    xPrev(:,1,1)    = x0;
    for i=2:1:DoP
        xPrev(:,1,i) = x(:,sizeSeg,i-1);
        lambdaNext(:,sizeSeg,i-1) = lambda(:,1,i);
    end

    parfor (i=1:1:DoP,numThreads)
%     for i=1:1:DoP
        lambda_i = lambda(:,:,i);
        mu_i = mu(:,:,i);
        u_i = u(:,:,i);
        x_i = x(:,:,i);
        p_i = p(:,:,i);

        xPrev_i      = xPrev(:,:,i);
        lambdaNext_i = lambdaNext(:,:,i);

        [KKTxEquation_i,KKTC_i,KKTHu_i,KKTlambdaEquation_i,L_i,LB_i] =...
            KKT_error_func(lambda_i,mu_i,u_i,x_i,p_i,xPrev_i,lambdaNext_i,rho,i);

        KKTxEquation(:,i)      = KKTxEquation_i;
        KKTC(:,i)              = KKTC_i;
        KKTHu(:,i)             = KKTHu_i;
        KKTlambdaEquation(:,i) = KKTlambdaEquation_i;
        L(:,i)                 = L_i;
    end

    KKTError.stateEquation   = max(KKTxEquation(:));
    KKTError.C               = max(KKTC(:));
    KKTError.Hu              = max(KKTHu(:));
    KKTError.costateEquation = max(KKTlambdaEquation(:));
    costL = sum(L(:));
    timerEnd = Timer();
    timeElapsed = timerEnd  - timerStart;
end

author: deng-haoyang

repo: ParNMPC


text:
classdef ProductComputerFactory < handle

    properties (Access = private)
        secondOrder
        fourthOrder
        secondOrderOut
        productComputer
    end

    methods (Access = public)

        function pComp = create(obj,C,e)
            obj.init(C,e)
            obj.createProductComputer()
            pComp = obj.getProductComputer();
        end
    end

    methods (Access = private)

        function init(obj,C,e)
           obj.fourthOrder = C;
           obj.secondOrder = e;
        end

        function  createProductComputer(obj)
            C = obj.fourthOrder;
            e = obj.secondOrder;

            if obj.isVoigt()
                p = ProductComputerForVoigt(C,e);
            elseif obj.isTensor()
                p = ProductComputerForTensor(C,e);
            else
                 error('Not possible to compute Product with these classes')
            end

            obj.productComputer = p;

        end

        function itIs = isVoigt(obj)
            isFourthVoigt = strcmp(obj.fourthOrder.getRepresentation(),'voigt');
            isSecondVoigt = strcmp(obj.secondOrder.getRepresentation(),'voigt');
            itIs = isFourthVoigt && isSecondVoigt;
        end

        function itIs = isTensor(obj)
            isFourthVoigt = strcmp(obj.fourthOrder.getRepresentation(),'tensor');
            isSecondVoigt = strcmp(obj.secondOrder.getRepresentation(),'tensor');
            itIs = isFourthVoigt && isSecondVoigt;
        end

        function p = getProductComputer(obj)
            p = obj.productComputer;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function adj = meshAdjacencyMatrix(faces, varargin)
%MESHADJACENCYMATRIX Compute adjacency matrix of a mesh from set of faces
%
%   ADJMAT = meshAdjacencyMatrix(FACES)
%   Returns a sparse NV-by-NV matrix (NV being the maximum vertex index)
%   containing vertex adjacency of the mesh represented by FACES.
%   FACES is either a NF-by-3, a NF-by-4 index array, or a Nf-by-1 cell
%   array.
%
%   Example
%     [v f] = createCube;
%     adj = meshAdjacencyMatrix(f);
%
%   See also
%     meshes3d, triangulateFaces, smoothMesh

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

% Ensures faces is a N-by-3 or N-by-4 array
if iscell(faces) || (isnumeric(faces) && size(faces, 2) > 4)
    faces = triangulateFaces(faces);
end

% forces faces to b efloating point array, for sparse function
if ~isfloat(faces)
    faces = double(faces);
end

% populate a sparse matrix
if size(faces, 2) == 3
    adj = sparse(...
        [faces(:,1); faces(:,1); faces(:,2); faces(:,2); faces(:,3); faces(:,3)], ...
        [faces(:,3); faces(:,2); faces(:,1); faces(:,3); faces(:,2); faces(:,1)], ...
        1.0);
elseif size(faces, 2) == 4
    adj = sparse(...
        [faces(:,1); faces(:,1); faces(:,2); faces(:,2); faces(:,3); faces(:,3); faces(:,4); faces(:,4)], ...
        [faces(:,4); faces(:,2); faces(:,1); faces(:,3); faces(:,2); faces(:,4); faces(:,3); faces(:,1)], ...
        1.0);
end

% remove double adjacencies
adj = min(adj, 1);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%% KUKA sunrise toolbox example.
% moving end-effector of the robot on an ellipse by utilizing the point to
% point elliptical motion functions of the KST

% Copy right: Mohammad SAFEEA
% 30-April-2018

global t_Kuka;
ip='172.31.1.147';
t_Kuka=net_establishConnection(ip);

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  disp('Connection could not be establised, script aborted');
  return;
end

%% Go to some initial configuration
disp('moving in joint space to initial configuration');
jPos={0., pi / 180 * 30, 0, -pi / 180 * 60, 0,...
                        pi / 180 * 90, 0};
disp(jPos);
relVel=0.15;
movePTPJointSpace( t_Kuka , jPos, relVel); % move to initial configuration

%% move a little bit back on the X direction
disp('moving -60 mm in the X direction')
deltaX=-60;deltaY=0;deltaZ=0.;
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
vel=50;
movePTPLineEefRelBase( t_Kuka , Pos, vel);

%% put the pen on the level of the page
disp('moving -85 mm in the Z direction')
deltaX=0;deltaY=0;deltaZ=-85.;
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
vel=50;
movePTPLineEefRelBase( t_Kuka , Pos, vel);
pause(1);

%% Define the ellipse,
disp('Drawing an ellipse in a plane parallel to XY axes of the base')
c=[0; 50]; % this is the displacement of the center of the ellipse with respect to the current position of EEF,
% taken in the XY plane of the robot base
ratio=0.5; % the radious ratio (a/b) of the ellipse
velocity=40;
accel=25;
theta=2*pi;
TefTool=eye(4);
movePTPEllipse_XY(t_Kuka,c,ratio,theta,velocity,accel,TefTool);

%% Define the ellipse, dimentsions are in (meter)
disp('Drawing an ellipse in a plane parallel to XY axes of the base')
c=[0; 50]; % this is the displacement of the center of the ellipse with respect to the current position of EEF,
% taken in the XZ plane of the robot base
ratio=0.5; % the radious ratio (a/b) of the ellipse
velocity=40;
accel=25;
theta=2*pi;
TefTool=eye(4);
movePTPEllipse_XZ(t_Kuka,c,ratio,theta,velocity,accel,TefTool);

%% Define the ellipse, dimentsions are in (meter)
disp('Drawing an ellipse in a plane parallel to XY axes of the base')
c=[0; 50]; % this is the displacement of the center of the ellipse with respect to the current position of EEF,
% taken in the YZ plane of the robot base
ratio=0.5; % the radious ratio (a/b) of the ellipse
velocity=40;
accel=25;
theta=2*pi;
TefTool=eye(4);
movePTPEllipse_YZ(t_Kuka,c,ratio,theta,velocity,accel,TefTool);

net_turnOffServer(t_Kuka);

repo: KST-Kuka-Sunrise-Toolbox


text:
function determ = chow_determinant ( alpha, beta, n )

%*****************************************************************************80
%
%% CHOW_DETERMINANT returns the determinant of the CHOW matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the ALPHA value.  A typical value is 1.0.
%
%    Input, real BETA, the BETA value.  A typical value is 0.0.
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  k = n - floor ( n / 2 );

  for i = 1 : k
    angle = i * pi / ( n + 2 );
    determ = determ * ( beta + 4.0 * alpha * ( cos ( angle ) )^2 );
  end

  determ  = determ * beta^( n - k );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fullpath = export_figure(handle, filename, format, varargin)
% export_figure Export a figure to various formats
%
% Utility function that exports a figure to various formats taking care
% of caching as well as some other options.
%
% Input:
% - handle (handle): Handle of a figure.
% - filename (string): Filename of the target file without extension.
% - format (string): Target format (fig, eps or png).
% - varargin[Cache] (boolean): Use caching. Disabled by default.
%
% Output:
% - fullpath (string): Path to the resulting file.
%

cache = false;

for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'cache'
            cache = varargin{i+1};
        otherwise
            error(['Unknown switch ', varargin{i},'!']) ;
    end
end

fullpath = [filename, '.', format];

if isa(handle, 'function_handle')
    autoclose = true;
    cache = true;
else
    autoclose = false;
end

if cache && exist(fullpath, 'file')
    return;
end;

if isa(handle, 'function_handle')
    handle = handle();
end;

switch lower(format)
    case 'fig'
        saveas(handle, filename, 'fig');
    case 'eps'
        allca = num2cell(findall(handle, 'type', 'axes'));
        allbackup = cellfun(@update_axis, allca, 'UniformOutput', false);
        print( handle, '-depsc', [filename, '.eps']);
        cellfun(@restore_axis, allca, allbackup, 'UniformOutput', false);
    case 'png'
        ah = get(handle, 'CurrentAxes');
        title = get(ah, 'Title');
        set(title, 'Visible', 'off');
        print( handle, '-dpng', '-r130', [filename, '.png']);
        set(title, 'Visible', 'on');
    otherwise
        error('Unknown format');
end;

if autoclose
    close(handle);
end

end

function [backup] = update_axis(ha)
    backup.fontSize = get(ha, 'FontSize');
    backup.fontWeight = get(ha, 'FontWeight');
    backup.lineWidth = get(ha, 'LineWidth');
    set(ha, 'FontSize', 10, 'FontWeight', 'bold', 'LineWidth', 2);
end

function [backup] = restore_axis(ha, backup)
    set(ha, 'FontSize', backup.fontSize, 'FontWeight', backup.fontWeight, 'LineWidth', backup.lineWidth);
end

author: votchallenge

repo: toolkit-legacy


text:
function [rRelmTest] = relm_LTest(vRatesH, vRatesN, nNumberSimulation, fMagThreshold, bOptimized, bDrawFigure)
% function [rRelmTest] = relm_LTest(vRatesH, vRatesN, nNumberSimulation, fMagThreshold, bOptimized, bDrawFigure)
% --------------------------------------------------------------------------------------------------------------
% Computation of the L-test for the RELM framework
%
% Input parameters:
%   vRatesH                       Matrix with rates of the test hypothesis
%   vRatesN                       Matrix with rates of the null hypothesis
%   nNumberSimulation             Number of random simulations
%   fMagThreshold                 Magnitude threshold (Use only bins with magnitude >= threshold
%   bOptimized                    0 (default): use a for loop, 1: matrix-wise calculation (needs a lot of memory)
%   bDrawFigure                   Draw the cumulative density plot after testing (default: off)
%
% Output paramters:
%   rRelmTest.fAlpha              Alpha-value of the cumulative density
%   rRelmTest.fBeta               Beta-value of the cumulative density
%   rRelmTest.vSimValues_H        Vector containing the sorted simulated numbers of events for the test hypothesis
%   rRelmTest.vSimValues_N        Vector containing the sorted simulated numbers of events for the null hypothesis
%   rRelmTest.nNumberSimulation   Number of random simulations
%   rRelmTest.fObservedData       Observed total number of events
%
% Copyright (C) 2002-2006 by Danijel Schorlemmer
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the
% Free Software Foundation, Inc.,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

% Exit on empty rate matrices
if isempty(vRatesH)  ||  isempty(vRatesN)
  rRelmTest.fAlpha = nan;
  rRelmTest.fBeta = nan;
  rRelmTest.vSimValues_H = nan;
  rRelmTest.vSimValues_N = nan;
  rRelmTest.nNumberSimulation = nan;
  rRelmTest.fObservedData = nan;
  return;
end

if ~exist('bDrawFigure')
  bDrawFigure = 0;
end

if ~exist('bOptimized')
  bOptimized = 0;
end

% Get the necessary data from the rate matrices and weight them properly
[vLambdaH, vLambdaN, vNumberQuake] = relm_PrepareData(vRatesH, vRatesN, fMagThreshold);
nNumberQuake = sum(vNumberQuake);

fLikelihood_H = sum(calc_logpoisspdf(vNumberQuake, vLambdaH));
fLikelihood_N = sum(calc_logpoisspdf(vNumberQuake, vLambdaN));

% Get the number of bins (rows)
[nRow, nColumn] = size(vLambdaH);

if bOptimized
  % Create the random numbers for the simulation
  vRandom = rand(nRow, nNumberSimulation);

  % Replicate the rate vectors
  vLambdaH = repmat(vLambdaH, 1, nNumberSimulation);
  vLambdaN = repmat(vLambdaN, 1, nNumberSimulation);

  % Compute the simulated number of events and sum them up
  vNum_H = poissinv(vRandom, vLambdaH);
  vNum_N = poissinv(vRandom, vLambdaN);
  vSimNum_H = sum(calc_logpoisspdf(vNum_H, vLambdaH)) - fLikelihood_H;
  vSimNum_N = sum(calc_logpoisspdf(vNum_N, vLambdaN)) - fLikelihood_N;
else
  % Create empty vectors for the total number of events
  vSimNum_H = [];
  vSimNum_N = [];

  % Loop over the simulations
  for nCnt = 1:nNumberSimulation
    % Create the random numbers for the simulation
    vRandom = rand(nRow, 1);

    % Compute the simulated number of events and sum them up
    vNum_H = poissinv(vRandom, vLambdaH);
    vNum_N = poissinv(vRandom, vLambdaN);
    vSimNum_H = [vSimNum_H; (sum(calc_logpoisspdf(vNum_H, vLambdaH)) - fLikelihood_H)];
    vSimNum_N = [vSimNum_N; (sum(calc_logpoisspdf(vNum_N, vLambdaN)) - fLikelihood_N)];
  end
end

% Sort them for the cumulative density plot
rRelmTest.vSimValues_H = sort(vSimNum_H);
rRelmTest.vSimValues_N = sort(vSimNum_N);

% Compute Alpha and Beta and store the important parameters
rRelmTest.fAlpha = sum(rRelmTest.vSimValues_N > 0)/nNumberSimulation;
rRelmTest.fBeta = sum(rRelmTest.vSimValues_H < 0)/nNumberSimulation;
rRelmTest.nNumberSimulation = nNumberSimulation;
rRelmTest.fObservedData = 0;

if bDrawFigure
  relm_PaintCumPlot(rRelmTest, 'Log-likelihood');
end


text:
function [jointTorques_star, LDot_tilde, L_tilde, intL_tilde, xiDes, f_des, xiDot_star] = momentumJerkControlYoga(intLTilde_angular, nu, J_CoM, JL, JR, xCoM, w_H_l_sole, w_H_r_sole, L, M, h, JDotL_nu, JDotR_nu, x_dx_ddx_dddx_CoM_des, ...
                                                                                gainsPCOM, gainsDCOM, LEFT_RIGHT_FOOT_IN_CONTACT, qj, qj_dqj_ddqj_des, impedances, dampings, ...
                                                                                contactForces_est, xi_star, xi, contactForces_star, Config, Gain, Reg)

    % MOMENTUMJERKCONTROL implements a momentum-based jerk controller.
    %
    % REFERENCES: IEEE-RAL 2018, "Momentum Control of an Underactuated Flying Humanoid Robot";
    %

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

    % parameters
    ndof            = size(M(7:end,7:end),1);

    % ROBUSTNESS w.r.t. modeling errors (inertial parameters)
    if Config.testRobustness

        M           = M * Config.robustnessFactor;
        h           = h * Config.robustnessFactor;
    end

    % get the robot total mass, and gravity forces, and feet positions
    m               = M(1,1);
    f_grav          = m*Config.GRAVITY_ACC*[0; 0; 1; zeros(3,1)];
    posLeftFoot     = w_H_l_sole(1:3,4);
    posRightFoot    = w_H_r_sole(1:3,4);

    % compute momentum references
    [LDDot_des, LDot_des, L_des, intL_des] = computeMomentumReferences(x_dx_ddx_dddx_CoM_des, m);

    %% %%%%%%%%%%%%%%% COMPUTE THE MOMENTUM ACCELERATION %%%%%%%%%%%%%%% %%
    %
    % WHILE BALANCING: the momentum acceleration is defined as follows:
    %
    %   LDDot = Ac * Beta* xiDot + AcDot * contactForces. (1)
    %
    % The contact forces derivative have been parametrized as follows:
    %
    %   contactForcesDot = Beta * xiDot
    %
    % xiDot is then used as control input. Beta is an invertible matrix.
    % It is therefore necessary to compute AcDot and Ac.

    % compute matrix Ac for both feet
    r_left          = posLeftFoot  - xCoM;
    r_right         = posRightFoot - xCoM;

    Ac_leftFoot     = [eye(3),           zeros(3);
                       wbc.skew(r_left), eye(3)];

    Ac_rightFoot    = [eye(3),            zeros(3);
                       wbc.skew(r_right), eye(3)];

    % compute matrix Beta_left and Beta_right
    Beta_left       = computeParametrizationGradient(xi(1:6),   Config);
    Beta_right      = computeParametrizationGradient(xi(7:end), Config);

    % multiplier of contact forces derivatives
    Ac              = [Ac_leftFoot .* LEFT_RIGHT_FOOT_IN_CONTACT(1), ...
                       Ac_rightFoot .* LEFT_RIGHT_FOOT_IN_CONTACT(2)];
    Ac_Beta         = [Ac_leftFoot * Beta_left .* LEFT_RIGHT_FOOT_IN_CONTACT(1), ...
                       Ac_rightFoot * Beta_right .* LEFT_RIGHT_FOOT_IN_CONTACT(2)];

    % compute the time derivative of r_left and r_right
    nu_leftFoot     = JL * nu;
    nu_rightFoot    = JR * nu;
    v_CoM           = J_CoM(1:3,:) * nu;
    rDot_left       = nu_leftFoot(1:3)  - v_CoM;
    rDot_right      = nu_rightFoot(1:3) - v_CoM;

    % compute AcDot
    AcDot_leftFoot  = [zeros(3),             zeros(3);
                       wbc.skew(rDot_left),  zeros(3)];

    AcDot_rightFoot = [zeros(3),             zeros(3);
                       wbc.skew(rDot_right), zeros(3)];

    AcDot           = [AcDot_leftFoot .* LEFT_RIGHT_FOOT_IN_CONTACT(1), ...
                       AcDot_rightFoot .* LEFT_RIGHT_FOOT_IN_CONTACT(2)];

    %% %%%%%%%%%%%%%%% MOMENTUM CONTROLLER DEFINITION  %%%%%%%%%%%%%%%%% %%
    %
    % The control input of Eq. (1) is:
    %
    %  u1 = xiDot
    %
    % For the purpose of this documentation, it suffices to know that the
    % robot momentum converges to the desired values if the following
    % equation is verified:
    %
    %  ATilde * u1 + deltaTilde = 0 (2)
    %
    % where: ATilde and deltaTilde must be properly calculated.
    % For details on the controller implementation, see the RAL 2018 paper
    % at: https://ieeexplore.ieee.org/document/7997895.

    % Compute the momentum error derivative/integral
    LDot_estimated  = Ac * [contactForces_est(1:6) .* LEFT_RIGHT_FOOT_IN_CONTACT(1); ...
                            contactForces_est(7:end) .* LEFT_RIGHT_FOOT_IN_CONTACT(2)] - f_grav;
    LDot_tilde      = LDot_estimated - LDot_des;
    L_tilde         = L - L_des;
    intL_tilde      = [(m * xCoM - intL_des(1:3)); intLTilde_angular];

    % Get the gains matrices
    KP_momentum     = blkdiag(diag(gainsPCOM), Gain.KP_AngularMomentum .*eye(3));
    KD_momentum     = blkdiag(diag(gainsDCOM), Gain.KD_AngularMomentum .*eye(3));

    % CONTROL LAW: momentum-based control with Lyapunov stability (IEEE-RAL 2018)
    KTilde          = KP_momentum + inv(Gain.KO_momentum) + KD_momentum;
    deltaTilde      = KTilde * L_tilde + (KD_momentum + eye(6)) * LDot_tilde + ...
                      KP_momentum * intL_tilde - LDDot_des + AcDot * contactForces_est;

    % Primary task control input
    pinvAc_Beta       = [wbc.pinvDamped(Ac_leftFoot * Beta_left,1e-5) .* (1-LEFT_RIGHT_FOOT_IN_CONTACT(2)); zeros(6)] + ...
                        [zeros(6); wbc.pinvDamped(Ac_rightFoot * Beta_right,1e-5) .* (1-LEFT_RIGHT_FOOT_IN_CONTACT(1))] + ...
                         wbc.pinvDamped(Ac_Beta,1e-5).* LEFT_RIGHT_FOOT_IN_CONTACT(1).* LEFT_RIGHT_FOOT_IN_CONTACT(2);

    xiDot_primaryTask = -pinvAc_Beta * deltaTilde;

    % Null space projector (forces)
    Null_Ac_Beta      = (eye(12) - pinvAc_Beta*Ac_Beta)*LEFT_RIGHT_FOOT_IN_CONTACT(1)*LEFT_RIGHT_FOOT_IN_CONTACT(2);

    %% %%%%%%%%%%%%%%%%% TORQUE CONTROLLER DEFINITION  %%%%%%%%%%%%%%%%% %%

    % Get the required matrices
    Jc                = [JL .*LEFT_RIGHT_FOOT_IN_CONTACT(1); JR .*LEFT_RIGHT_FOOT_IN_CONTACT(2)];
    JcDot_nu          = [JDotL_nu .*LEFT_RIGHT_FOOT_IN_CONTACT(1); JDotR_nu .*LEFT_RIGHT_FOOT_IN_CONTACT(2)];
    B                 = [zeros(6,ndof); eye(ndof)];
    Lambda            = Jc/M*B;
    pinvLambda        = wbc.pinvDamped(Lambda, Reg.pinvDamp);
    NLambda           = eye(ndof) - pinvLambda*Lambda;

    % terms required for calculating u_0
    qjDot             = nu(7:end);
    KP_torqueControl  = diag(impedances);
    KD_torqueControl  = diag(dampings);
    jointPos_err      = qj    - qj_dqj_ddqj_des(:,1);
    jointVel_err      = qjDot - qj_dqj_ddqj_des(:,2);

    % consitency with the current gain tuning
    USE_MBAR          = false;
    USE_ACC           = 0;

    if USE_MBAR

        MBar          = (M(7:end,7:end) - M(7:end,1:6)/M(1:6,1:6)*M(1:6,7:end)); %#ok<UNRCH>
    else
        MBar          = eye(ndof);
    end

    u_0               = MBar*(USE_ACC*qj_dqj_ddqj_des(:,3) - KD_torqueControl * jointVel_err - KP_torqueControl * jointPos_err);

    % terms required for calculating tau_0
    tau_0             = h(7:end) -M(7:end,1:6)/M(1:6,1:6)*h(1:6) ...
                       -transpose(Jc(:,7:end))*contactForces_est + ...
                        M(7:end,1:6)/M(1:6,1:6)*transpose(Jc(:,1:6))*contactForces_est + u_0;

    jointTorques_star = pinvLambda*(Jc/M*(h-Jc'*contactForces_star) -JcDot_nu) + NLambda*tau_0;

    %% %%%%%%%%%%%%%%%%%%%%% NULL SPACE OF MOMENTUM %%%%%%%%%%%%%%%%%%%% %%

    % Now, rewrite the joint torques equations as:
    %
    %    tau = Sigma * f_star + gamma_full
    %
    % where we isolated the terms multiplied by the forces in the
    % joint torques equations. Finally we compute the time derivative of tau:
    %
    %    tauDot = SigmaDot * f_star + Sigma * fDot_star + gamma_fullDot
    %    tauDot ~ Sigma * fDot_star (approximated) + gamma_redDot
    %
    % Then, consider the following Lyapunov function candidate:
    %
    %    V    = 1/2*|tau|^2
    %    VDot = tau^T*tauDot
    %
    % It is clear that choosing tauDot = -K_tau*tau will ensures the
    % convergence of the joint torques to zero.
    %

    % theoretically, the forces inside tau_0 are the measured ones, not the
    % desired. Therefore only the derivative of the desired forces should
    % be used to minimize the joint torques, however I am not really sure
    % about this.
    DEACTIVATE_NULL_SPACE_FORCES = 0;

    Sigma               = -pinvLambda*(Jc/M*Jc') -DEACTIVATE_NULL_SPACE_FORCES.*NLambda*(transpose(Jc(:,7:end)) - M(7:end,1:6)/M(1:6,1:6)*transpose(Jc(:,1:6)));
    Sigma_Beta          =  Sigma * Null_Ac_Beta;
    gamma_redDot        =  Sigma * xiDot_primaryTask -Gain.useVelFeedbackForMinTorques*NLambda*MBar*KP_torqueControl*jointVel_err;
    K_tau               =  Gain.K_tau;

    % Secondary task
    pinvSigma_Beta      =  wbc.pinvDamped(Sigma_Beta, Reg.pinvDamp);
    xiDot_secondaryTask = -pinvSigma_Beta*(K_tau*jointTorques_star + gamma_redDot);

    % FINAL control input (Parametrized)
    pinvAc              = [wbc.pinvDamped(Ac_leftFoot,1e-5) .* (1-LEFT_RIGHT_FOOT_IN_CONTACT(2)); zeros(6)] + ...
                          [zeros(6); wbc.pinvDamped(Ac_rightFoot,1e-5) .* (1-LEFT_RIGHT_FOOT_IN_CONTACT(1))] + ...
                           wbc.pinvDamped(Ac,1e-5).* LEFT_RIGHT_FOOT_IN_CONTACT(1).* LEFT_RIGHT_FOOT_IN_CONTACT(2);

    f_des               = pinvAc*(LDot_des + f_grav);

    if LEFT_RIGHT_FOOT_IN_CONTACT(1) > 0.5

        xi_des_left = fromForcesToParametrization(f_des(1:6),Config);
    else
        xi_des_left = zeros(6,1);
    end
    if LEFT_RIGHT_FOOT_IN_CONTACT(2) > 0.5

        xi_des_right = fromForcesToParametrization(f_des(7:end),Config);
    else
        xi_des_right = zeros(6,1);
    end

    xiDes               = [xi_des_left;xi_des_right];
    Ke                  = Gain.Ke;
    xiDot_star          = xiDot_primaryTask + Null_Ac_Beta*xiDot_secondaryTask -Ke*(xi_star-[xiDes(1:6)*LEFT_RIGHT_FOOT_IN_CONTACT(1);xiDes(7:end)*LEFT_RIGHT_FOOT_IN_CONTACT(2)]);
end

author: robotology

repo: whole-body-controllers


text:
function C=searchintial(X,method,varargin)
switch lower(method(1))
    case 's'
        K=varargin{1};
        C=X(randsample(size(X,1),K),:);
    case 'u'
        Xmins=min(X,[],1);
        Xmaxs=max(X,[],1);
        K=varargin{1};
        C=unifrnd(Xmins(ones(K,1),:), Xmaxs(ones(K,1),:));
end

author: decouples

path: github-repos/MATLAB/decouples-Matlab_deep_learning/Matlab_deep_learning-1b823b82686080e32b03e1f1a4648896bd6e3c44/第 22 章 基于 K-means 聚类算法的图像区域分割/searchintial.m


text:
%% FTAlignStack

function S=AlignStack(stack)

%% 081010 Tobias Henn

S=stack;
stackcontainer=stack.spectr;

clear S.spectr

[ymax,xmax,emax]=size(stackcontainer);

xresolution=S.Xvalue/xmax;
yresolution=S.Yvalue/ymax;

center=ceil(emax/4*3);

spectr=zeros(ymax,xmax,emax);

shifts=zeros(emax,4);

%calculate image shifts for each energy,perform shift with FT method

for k=1:emax

    shifts(k,:)=dftregistration(fft2(stackcontainer(:,:,center)),fft2(stackcontainer(:,:,k)),10);
    spectr(:,:,k)=FTMatrixShift(stackcontainer(:,:,k),-shifts(k,3),-shifts(k,4));

end

%Reduce image size

shiftymax=ceil(max(shifts(:,3)));
shiftxmax=ceil(max(shifts(:,4)));
shiftymin=ceil(abs(min(shifts(:,3))));
shiftxmin=ceil(abs(min(shifts(:,4))));

shiftmatrix=zeros(ymax-shiftymin-shiftymax,xmax-shiftxmax-shiftxmin,emax);

shiftmatrix(:,:,:)=spectr((1+shiftymax):(ymax-shiftymin),(1+shiftxmax):(xmax-shiftxmin),:);

S.spectr=abs(shiftmatrix);

S.Xvalue=size(S.spectr,2)*xresolution;
S.Yvalue=size(S.spectr,1)*yresolution;

return


text:
% Copyright (C) 2016, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function speed_demo_scara
global robot
% joint position
q = [pi/4 pi/4 0.2 0]';
% joint speed
qd = [1 1 1 1]';
%robot = load_robot

J = manipulator_jacobian(robot, q);

Ve=J*qd;

%plot speed
T = directkinematic(robot, q);
p0 = T(1:3,4);
drawrobot3d(robot, q)
draw_vector(Ve(1:3), p0, 'linear speed V', 2)
draw_vector(Ve(4:6), p0, 'angular speed W', 1)

repo: ARTE


text:
%% Sequential update for Gaussian
close all; clear;
d = 2;
n = 100;
X = randn(d,n);
x = randn(d,1);

mu = mean(X,2);
Xo = bsxfun(@minus,X,mu);
Sigma = Xo*Xo'/n;
p1 = logGauss(x,mu,Sigma);

gauss = Gauss(X(:,3:end)).addSample(X(:,1)).addSample(X(:,2)).addSample(X(:,3)).delSample(X(:,3));
p2 = gauss.logPdf(x);
maxdiff(p1,p2)
%% Sequential update for Gaussian-Wishart
close all; clear;
d = 2;
n = 100;
X = randn(d,n);
x = randn(d,1);

kappa0 = 1;
m0 = zeros(d,1);
nu0 = d;
S0 = eye(d);

xbar = mean(X,2);
kappa = kappa0+n;
nu = nu0+n;
m = (n*xbar+kappa0*m0)/kappa;
Xo = bsxfun(@minus,X,m);
X0 = m0-m;
S = S0+Xo*Xo'+kappa0*(X0*X0');

v = (nu-d+1);
r = (1+1/kappa)/v;
p1 = logSt(x,m,r*S,v);

gw0 = GaussWishart(kappa0,m0,nu0,S0);
gw0 = gw0.addData(X);
p0 = gw0.logPredPdf(x);

gw = GaussWishart(kappa0,m0,nu0,S0);
for i=1:n
    gw = gw.addSample(X(:,i));
end
p2 = gw.logPredPdf(x);
maxdiff(p1,p2)
%

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [ out ] = VideoToHistogramList( vidPath, s, w )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

video = VideoReader(vidPath);
NumberOfFramesToRead = ceil(video.NumberOfFrames/10);
HistsPerFrame = video.Height*video.Width/(s^2);
out = zeros(256/w,3,NumberOfFramesToRead*HistsPerFrame);

for i = 1:10:video.NumberOfFrames
   frame = read(video, i);
   out(:,:,(i-1)/10*HistsPerFrame+1:((i-1)/10+1)*HistsPerFrame) = Histograms1D(frame, s, w);
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function gr=bvalue(catalogObject, mcType, manual_on)
    %BVALUE evaluate b-value, a-value and magnitude of completeness
    % of an earthquake catalog stored in a Catalog object.
    %
    % gr = BVALUE(catalogObject, MCTYPE) produces a Gutenberg-Richter type plot
    %    with the best fit line and display of b-,a-values and Mc
    %    for catalogObject. These values are also returned in a structure.
    %    MCTYPE is a number from 1-5
    %    to select the algorithm used for calculation of the
    %    magnitude of completeness. Options are:
    %
    %    1: Maximum curvature
    %    2: Fixed Mc = minimum magnitude (Mmin)
    %    3: Mc90 (90% probability)
    %    4: Mc95 (95% probability)
    %    5: Best combination (Mc95 - Mc90 - maximum curvature)
    %
    % * Note: it seems only 1 only really works, and 5 is same as 1 *'
    %
    % gr = BVALUE(catalogObject, MCTYPE, manual_on) where the value of manual_on
    %    computes as true will give the user the ability to manually pick a
    %    linear segment on the graph too (which is then plotted with a
    %    green line). The manual Mc and bvalue are returned in the gr
    %    structure as gr.Mc_manual and gr.bvalue_manual

    % Liberally adapted from original code in ZMAP.
    % Author: Silvio De Angelis, 27/07/2012 00:00:00
    % Modified and included in Catalog by Glenn Thompson,
    % 14/06/2014

    % This program is free software; you can redistribute it and/or modify
    % it under the terms of the GNU General Public License as published by
    % the Free Software Foundation; either version 2 of the License, or
    % (at your option) any later version.
    %
    % This program is distributed in the hope that it will be useful,
    % but WITHOUT ANY WARRANTY; without even the implied warranty of
    % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    % GNU General Public License for more details.
    %
    % You should have received a copy of the GNU General Public License
    % along with this program; if not, write to the
    % Free Software Foundation, Inc.,
    % 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

    if nargin < 2
        disp('* Note: it seems only 1 only really works, and 5 is same as 1 *')
        mcType = menu('mcType can be:','Maximum curvature','Fixed Mc = minimum magnitude (Mmin)', ...
            'Mc90 (90% probability)', 'Mc95 (95% probability)', ...
            'Best combination (Mc95 - Mc90 - maximum curvature)')
    end
    if ~exist('manual_on','var')
        manual_on = false;
    end

    % form magnitude vector - removing any NaN values with find
    good_magnitude_indices = find(catalogObject.mag > -3.0);
    mag = catalogObject.mag(good_magnitude_indices);
    %MIN AND MAX MAGNITUDE IN catalogObject
    minimum_mag = min(mag);
    maximum_mag = max(mag);

    %COUNT EVENTS IN EACH MAGNITUDE BIN
    [bval, xt2] = hist(mag, (minimum_mag:0.1:maximum_mag));

    %CUMULATIVE NUMBER OF EVENTS IN EACH MAGNITUDE BIN
    bvalsum = cumsum(bval);

    %NUMBER OF EVENTS IN EACH BIN IN REVERSE ORDER
    bval2 = bval(length(bval):-1:1);

    %NUMBER OF EVENTS IN EACH MAGNITUDE BIN IN REVERSE ORDER
    bvalsum3 = cumsum(bval(length(bval):-1:1));

    %BINS IN REVERSE ORDER
    xt3 = (maximum_mag:-0.1:minimum_mag);
    backg_ab = log10(bvalsum3);

    %CREATE FIGURE WINDOW AND MAKE FREQUENCY-MAGNITUDE PLOT
    figure('Color','w','Position',[0 0 600 600])

    pl = semilogy(xt3,bvalsum3,'sb'); % semilogy is same as plot, except a log (base10) scale is used for Y-axis
    set(pl, 'LineWidth', [1.0],'MarkerSize', [10],'MarkerFaceColor','r','MarkerEdgeColor','k');
    axis square
    hold on

    pl1 = semilogy(xt3,bval2,'^b');
    set(pl1, 'LineWidth',[1.0],'MarkerSize',[10],'MarkerFaceColor','w','MarkerEdgeColor','k');
    xlabel('Magnitude','Fontsize', 12)
    ylabel('Cumulative Number','Fontsize',12)
    set(gca,'visible','on','FontSize',12,'FontWeight','normal',...
        'FontWeight','bold','LineWidth',[1.0],'TickDir','in','Ticklength',[0.01 0.01],...
        'Box','on','Tag','cufi','color','w')

    %ESTIMATE B-VALUE (MAX LIKELIHOOD ESTIMATE)
    Nmin = 10;
    fMccorr = 0;
    fBinning = 0.1;

    if length(mag) >= Nmin

        %GOODNESS-OF-FIT TO POWER LAW
        %%%%%%%%%%%%%%%%%% mcperc_ca3.m start %%%%%%%%%%%%%%%%%%%%
        % This is a completeness determination test

        [bval,xt2] = hist(mag,-2:0.1:6);
        l = max(find(bval == max(bval)));
        magco0 =  xt2(l);

        dat = [];

        %for i = magco0-0.6:0.1:magco0+0.2
        for i = magco0-0.5:0.1:magco0+0.7
            l = mag >= i - 0.0499;
            nu = length(mag(l));
            if length(mag(l)) >= 25;
                %[bv magco stan av] =  bvalca3(catZmap(l,:),2,2);
                [mw bv2 stan2 av] =  Catalog.bvalue_lib.bmemag(mag(l));
                Catalog.bvalue_lib.synthb_aut;
                dat = [ dat ; i res2];
            else
                dat = [ dat ; i nan];
            end

        end

        j =  min(find(dat(:,2) < 10 ));
        if isempty(j) == 1; Mc90 = nan ;
        else;
            Mc90 = dat(j,1);
        end

        j =  min(find(dat(:,2) < 5 ));
        if isempty(j) == 1; Mc95 = nan ;
        else;
            Mc95 = dat(j,1);
        end

        j =  min(find(dat(:,2) < 10 ));
        if isempty(j) == 1; j =  min(find(dat(:,2) < 15 )); end
        if isempty(j) == 1; j =  min(find(dat(:,2) < 20 )); end
        if isempty(j) == 1; j =  min(find(dat(:,2) < 25 )); end
        j2 =  min(find(dat(:,2) == min(dat(:,2)) ));
        %j = min([j j2]);

        Mc = dat(j,1);
        magco = Mc;
        prf = 100 - dat(j2,2);
        if isempty(magco) == 1; magco = nan; prf = 100 -min(dat(:,2)); end
        %display(['Completeness Mc: ' num2str(Mc) ]);
        %%%%%%%%%%%%%%%%%% mcperc_ca3.m end %%%%%%%%%%%%%%%%%%%%%%

        %CALCULATE MC
        [fMc] = Catalog.bvalue_lib.calc_Mc(mag, mcType, fBinning, fMccorr);
        l = mag >= fMc-(fBinning/2);
        if length(mag(l)) >= Nmin
            [fMeanMag, fBValue, fStd_B, fAValue] =  Catalog.bvalue_lib.calc_bmemag(mag(l), fBinning);
        else
            [fMc, fBValue, fStd_B, fAValue] = deal(NaN);
        end

        %STANDARD DEV OF a-value SET TO NAN;
        [fStd_A, fStd_Mc] = deal(NaN);

    else
        [fMc, fStd_Mc, fBValue, fStd_B, fAValue, fStd_A, ...
            fStdDevB, fStdDevMc] = deal(NaN);
    end

    magco = fMc;
    index_low=find(xt3 < magco+.05 & xt3 > magco-.05);
try
    mag_hi = xt3(1);
    index_hi = 1;
    mz = xt3 <= mag_hi & xt3 >= magco-.0001;
    mag_zone=xt3(mz);
    y = backg_ab(mz);

    %PLOT MC IN FIGURE
    Mc = semilogy(xt3(index_low),bvalsum3(index_low)*1.5,'vk');
    set(Mc,'LineWidth',[1.0],'MarkerSize',7)
    Mc = text(xt3(index_low)+0.2,bvalsum3(index_low)*1.5,'Mc');
    set(Mc,'FontWeight','normal','FontSize',12,'Color','k')

    %CREATE AND PLOT FIT LINE
    sol_type = 'Maximum Likelihood Solution';
    bw=fBValue;
    aw=fAValue;
    ew=fStd_B;
    p = [ -1*bw aw];
    f = polyval(p,mag_zone);
    f = 10.^f;
    hold on
    ttm= semilogy(mag_zone,f,'k');
    set(ttm,'LineWidth',[2.0])
    std_backg = ew;

    %ERROR CALCULATIONS
    %b = mag;
    bv = [];
    si = [];

    set(gca,'XLim',[min(mag)-0.5  max(mag+0.5)])
    %set(gca,'YLim',[0.9 length(mag+30)*2.5]);

    p=-p(1,1);
    p=fix(100*p)/100;
    tt1=num2str(bw,3);
    tt2=num2str(std_backg,1);
    tt4=num2str(bv,3);
    tt5=num2str(si,2);
    tmc=num2str(magco,2);
    rect=[0 0 1 1];
    h2=axes('position',rect);
    set(h2,'visible','off');
    t=catalogObject.gettimerange();
    a0 = aw-log10((t(2)-t(1))/365);

    text(.53,.88, ['b-value = ',tt1,' +/- ',tt2,',  a value = ',num2str(aw,3)],'FontSize',12);
    text(.53,.85,sol_type,'FontSize',12 );
    text(.53,.82,['Magnitude of Completeness = ',tmc],'FontSize',12);

    %% Added by Glenn 2018-05-01 to return a structure
    gr.bvalue = str2num(tt1);
    gr.bvalue_error = str2num(tt2);
    gr.avalue = aw;
    gr.Mc = str2num(tmc);

    %% Manually fit a line (added by Glenn 2018-05-01)
    if manual_on
        [xmag, yN] = ginput(2);
        slope = (log10(yN(1)) - log10(yN(2) )) / (xmag(2) - xmag(1) )
        plot(xmag, yN, 'g');
        gr.Mc_manual = xmag(1);
        gr.bvalue_manual = slope;
        text(xmag(1) + (xmag(2)-xmag(1))*0.5, yN(1) + (yN(2)-yN(1))*0.5, sprintf('manual b=%.2f',slope),'Color','g');
    end
catch
    gr.bvalue = NaN;
    gr.Mc = NaN;
    gr.avalue = NaN;
    gr.bvalue_error = NaN;
end

end

author: geoscience-community-codes

repo: GISMO


text:
function [OffDec,OffMask] = Operator(Problem,ParentDec,ParentMask,rbm,dae,Site,allZero,allOne)
% The operator of MOEA/PSL
%
%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Parameter setting
    Parent1Mask = ParentMask(1:end/2,:);
    Parent2Mask = ParentMask(end/2+1:end,:);
    Parent1Dec  = ParentDec(1:end/2,:);
    Parent2Dec  = ParentDec(end/2+1:end,:);

    %% Binary variation
    if any(Site)
        other   = ~allZero & ~allOne;
        OffTemp = BinaryCrossover(rbm.reduce(Parent1Mask(Site,other)),rbm.reduce(Parent2Mask(Site,other)));
        OffTemp = rbm.recover(OffTemp);
        OffMask = false(size(OffTemp,1),size(Parent1Mask,2));
        OffMask(:,other)  = OffTemp;
        OffMask(:,allOne) = true;
    else
        OffMask = [];
    end
    OffMask = [OffMask;BinaryCrossover(Parent1Mask(~Site,:),Parent2Mask(~Site,:))];
    OffMask = BinaryMutation(OffMask);

    %% Real variation
    if any(Problem.encoding~=4)
        if any(Site)
            OffDec = RealCrossover(dae.reduce(Parent1Dec(Site,:)),dae.reduce(Parent2Dec(Site,:)));
            OffDec = dae.recover(OffDec);
        else
            OffDec = [];
        end
        OffDec = [OffDec;RealCrossover(Parent1Dec(~Site,:),Parent2Dec(~Site,:))];
        OffDec = RealMutation(OffDec,Problem.lower,Problem.upper);
        OffDec(:,Problem.encoding==4) = 1;
    else
        OffDec = ones(size(OffMask));
    end
end

function Offspring = BinaryCrossover(Parent1,Parent2)
% Uniform crossover

    k = rand(size(Parent1)) < 0.5;
    Offspring1    = Parent1;
    Offspring2    = Parent2;
    Offspring1(k) = Parent2(k);
    Offspring2(k) = Parent1(k);
    Offspring     = [Offspring1;Offspring2];
end

function Offspring = BinaryMutation(Offspring)
% Bitwise mutation

    Site = rand(size(Offspring)) < 1/size(Offspring,2);
    Offspring(Site) = ~Offspring(Site);
end

function Offspring = RealCrossover(Parent1,Parent2)
% Simulated binary crossover

    disC  = 20;
    [N,D] = size(Parent1);
    beta  = zeros(N,D);
    mu    = rand(N,D);
    beta(mu<=0.5) = (2*mu(mu<=0.5)).^(1/(disC+1));
    beta(mu>0.5)  = (2-2*mu(mu>0.5)).^(-1/(disC+1));
    beta = beta.*(-1).^randi([0,1],N,D);
    Offspring = [(Parent1+Parent2)/2+beta.*(Parent1-Parent2)/2
                 (Parent1+Parent2)/2-beta.*(Parent1-Parent2)/2];
end

function Offspring = RealMutation(Offspring,Lower,Upper)
% Polynomial mutation

    disM  = 20;
    [N,D] = size(Offspring);
    Lower = repmat(Lower,N,1);
    Upper = repmat(Upper,N,1);
    Site  = rand(N,D) < 1/D;
    mu    = rand(N,D);
    temp  = Site & mu<=0.5;
    Offspring       = min(max(Offspring,Lower),Upper);
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                      (1-(Offspring(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp = Site & mu>0.5;
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                      (1-(Upper(temp)-Offspring(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
    Offspring       = min(max(Offspring,Lower),Upper);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function output = bicubic_interpolation_at(input,uu,vv,nx,ny,border_out,BOUNDARY_CONDITION)
output = 0.0;
sx = 0;
sy = 0;
if(uu < 0)
    sx = -1;
else
    sx = 1;
end
if(vv < 0)
    sy = -1;
else
    sy = 1;
end
out = 0;

switch(BOUNDARY_CONDITION)
    case 0
          [out,x]   = neumann_bc(uu, nx);
         [out,y]   = neumann_bc(vv, ny);
         [out,mx]  = neumann_bc(uu - sx, nx);
         [out,my]  = neumann_bc(vv - sx, ny);
         [out,dx]  = neumann_bc( uu + sx, nx);
         [out,dy]  = neumann_bc( vv + sy, ny);
         [out,ddx] = neumann_bc( uu + 2*sx, nx);
         [out,ddy] = neumann_bc( vv + 2*sy, ny);
    case 1
        [out,x]   =periodic_bc(uu, nx);
         [out,y]   = periodic_bc(vv, ny);
         [out,mx]  = periodic_bc(uu - sx, nx);
         [out,my]  = periodic_bc(vv - sx, ny);
         [out,dx]  = periodic_bc( uu + sx, nx);
         [out,dy]  = periodic_bc( vv + sy, ny);
         [out,ddx] = periodic_bc( uu + 2*sx, nx);
         [out,ddy] = periodic_bc( vv + 2*sy, ny);
    case 2
         [out,x]  = symmetric_bc(uu, nx);
         [out,y]   = symmetric_bc(vv, ny);
         [out,mx]  = symmetric_bc(uu - sx, nx);
         [out,my]  = symmetric_bc(vv - sx, ny);
         [out,dx]  = symmetric_bc( uu + sx, nx);
         [out,dy]  = symmetric_bc( vv + sy, ny);
         [out,ddx] = symmetric_bc( uu + 2*sx, nx);
         [out,ddy] = symmetric_bc( vv + 2*sy, ny);
    otherwise
        [out,x]   = neumann_bc(uu, nx);
         [out,y]   = neumann_bc(vv, ny);
         [out,mx]  = neumann_bc(uu - sx, nx);
         [out,my]  = neumann_bc(vv - sx, ny);
         [out,dx]  = neumann_bc( uu + sx, nx);
         [out,dy]  = neumann_bc( vv + sy, ny);
         [out,ddx] = neumann_bc( uu + 2*sx, nx);
         [out,ddy] = neumann_bc( vv + 2*sy, ny);
         if((out == 1) && (border_out == 0))

             output = 0.0;
         else
            p11 = input(mx  + nx * my);
         p12 = input(x   + nx * my);
          p13 = input(dx  + nx * my);
          p14 = input(ddx + nx * my);

          p21 = input(mx  + nx * y);
          p22 = input(x   + nx * y);
          p23 = input(dx  + nx * y);
          p24 = input(ddx + nx * y);

          p31 = input(mx  + nx * dy);
          p32 = input(x   + nx * dy);
          p33 = input(dx  + nx * dy);
          p34 = input(ddx + nx * dy);

          p41 = input(mx  + nx * ddy);
          p42 = input(x   + nx * ddy);
          p43 = input(dx  + nx * ddy);
          p44 = input(ddx + nx * ddy);
             pol = [p11, p21, p31, p41;
            p12, p22, p32, p42;
            p13, p23, p33, p43;
            p14, p24, p34, p44];
      output =   bicubic_interpolation_cell(pol, uu-x, vv-y);
         end

end


text:
clear
clc

% mysql default configure
mdl_data = get_default_mysql_data('iris');
class_index = 1;
node_sizes = get_node_sizes(mdl_data);
% make a dag
[mdl_data_row,mdl_data_col] = size(mdl_data);

% make a naive struct
dag = mk_naive_struct(length(node_sizes),class_index);

%
% order = catty_order_algorithm(mdl_data,'class_node',1,'subtract_type','sum');
% dag = learn_struct_K2(mdl_data,node_size,order);
% draw_graph(dag);
%
% % make a bnet
bnet = get_bnet(mdl_data,dag,'tabular');
% % set classified times
times = 10;
class_test_algorithm = 'CV-5';
[correct_rate, confusion_matrix] = times_classified_test(bnet,mdl_data,...
    class_index,class_test_algorithm,times);
%
[avg_correct_rate,avg_confusion_matrix] ...
    = computer_avg_classified(correct_rate,confusion_matrix);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function runtest(opt)
% Run examples in the reference:
%   B.N. Li, C.K. Chui, S. Chang, S.H. Ong (2011) Integrating spatial fuzzy
%   clustering with level set methods for automated medical image
%   segmentation. Computers in Biology and Medicine 41(1) 1-10.
%--------------------------------------------------------------------------

if opt==1
    img=imread('ctlivertumor.bmp');
    ncluster=3;
elseif opt==2
    img=imread('mrihead.bmp');
    ncluster=4;
else
    error('Invalid opt: 1 or 2 only!')
end

MF = SFCM2D(img,ncluster);

figure
subplot(231); imshow(img,[])
for i=1:ncluster
    imgfi=reshape(MF(i,:,:),size(img,1),size(img,2));
    subplot(2,3,i+1); imshow(imgfi,[])
    title(['Index No: ' int2str(i)])
end

temp=1;
while temp
    nopt = input('Input the Index No that you are interested\n');
    if ~isempty(nopt), temp=0; end
end

close(gcf);

imgfcm=reshape(MF(nopt,:,:),size(img,1),size(img,2));

fuzzyLSM(img,imgfcm,.5);

author: Sable

repo: mcbench-benchmarks


text:
function result=run(obj,result,test_partition_index,test_partition_count)
% Run the test suite
%
% result=run(obj[,result])
%
% Inputs:
%   obj                     MoxUnitTestSuite instance with tests to be run.
%   result                  MoxUnitTestResult instance to which test
%                           results are to be reported (default: empty
%                           MoxUnitTestResult instance).
%   partition_index         } Either both or neither of these arguments
%   partition_count         } must be passed to this function. It causes
%                             a subset of all tests to be run, namely the
%                             ones indexed by:
%                                partition_index+partition_count*K
%                             for all values of K. Default values are
%                             partition_index=1 and partition_count=1,
%                             meaning that all tests are run. A use case is
%                             parallelization of test cases over multiple
%                             processes.
%
% Output:
%   result          MoxUnitTestResult containing tests results
%                   after running the tests in obj.
%
% See also: MoxUnitTestResult
%
% NNO 2015

    if nargin<2
        result=MOxUnitTestReport();
    end

    if nargin<3
        test_partition_index=1;
        test_partition_count=1;
    elseif nargin<4
            error(['With third argument provided, the fourth argument '...
                        'is required']);
    end

    ensure_integer(test_partition_index,'test_partition_index');
    ensure_integer(test_partition_count,'test_partition_count');

    if test_partition_index>test_partition_count
        error('Test index %d greater than number of partitions %d',...
                test_partition_index,test_partition_count);
    end

    n_nodes=countTestNodes(obj);
    for j=test_partition_index:test_partition_count:n_nodes
        test_case=getTestNode(obj,j);
        result=run(test_case, result);
    end

function ensure_integer(value, label)
    msg='';

    if ~isnumeric(value) || ~isscalar(value)
        msg='input must be scalar numeric';
    elseif round(value)~=value
        msg='input must be integer';
    elseif ~(value>0)
        msg='input must be positive';
    end

    if ~isempty(msg)
        error('''%s'' argument: %s',label,msg);
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


function box_overlaps=get_gt_overlap_box(reg2sp, sp, instimg)
if(all(instimg==0)) box_overlaps=zeros(0,size(reg2sp,2)); return; end
boxes=get_region_boxes(sp, reg2sp);
instimg=double(instimg);
insts=unique(instimg(instimg~=0));
gt_boxes=zeros(numel(insts),4);
for i=1:numel(insts)
    [I,J]=find(instimg==insts(i));
    I=I(:); J=J(:);
    gt_boxes(i,:)=[min(J) min(I) max(J) max(I)];
end
boxes(:,3:4)=boxes(:,3:4)-boxes(:,1:2)+1;
gt_boxes(:,3:4)=gt_boxes(:,3:4)-gt_boxes(:,1:2)+1;
int=rectint(boxes, gt_boxes);
uni=bsxfun(@plus, prod(boxes(:,3:4),2), prod(gt_boxes(:,3:4),2)')-int;
box_overlaps=int./(uni+double(uni~=0));
box_overlaps=box_overlaps';


text:
% test single frame
caffe('set_batch_size', 8);
for i = 1 : 8 : 128
    box = zeros(1, 1, 4, 8, 'single');
    image_batch = zeros(100, 100, 3, 8, 'single');
    for j = i : i + 8
        im = imread(['/home/winsty/caffe/dump/' num2str(j) '.jpg']);
        im = single(im);
        im = im(:, :, [3 2 1]) - 120;
        image_batch(:, :, :, j - i + 1) = permute(im, [2, 1, 3]);
        % dummy target
        box(:, :, :, j - i + 1) = [0, 0, -1, -1];
    end
    input = {image_batch; box};
    caffe('forward', input);
    fea = caffe('extract_feature', 'fc11');
    fea = fea{1};
    fea = reshape(fea, [50, 50, 8]);
    fea = 1 ./ (1 + exp(-fea));
    j = i + 7;
    im = single(imread(['/home/winsty/caffe/dump/' num2str(j) '.jpg'])) / 255;
    mask = imresize(fea(:, :, j - i + 1)', [100, 100]);
    for k = 1 : 3
        masked(:, :, k) = min(im(:, :, k), mask);
    end
    h = figure; imshow([im, masked]);
    saveas(h, [num2str(j) '.jpg']);
end

author: flyers

repo: drone-tracking


text:
function y = atan_pos(x,rnd)
% Rigorous atan(x) for nonnegative double vector x, rounding corresponding
%   to rnd
% For internal use in atan, asin, acos, asec, acsc

% written  12/30/98     S.M. Rump
% modified 08/31/99     S.M. Rump  extreme input, improved speed and accuracy
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  INTLAB_STDFCTS_PI = getappdata(0,'INTLAB_STDFCTS_PI');

  y = x;

  % transformation of large arguments
  index = ( x>=4 );
  if any(index(:))
    setround(-rnd)
    xx = x(index);
    xx = (-2) ./ ( 1./xx - xx );       % with correct rounding -rnd
    yy = atan_pos_small(xx,-rnd);
    setround(rnd)
    if rnd==-1
      y(index) = INTLAB_STDFCTS_PI.PI2INF - yy/2;
    else
      y(index) = INTLAB_STDFCTS_PI.PI2SUP - yy/2;
    end
  end

  index = ~index;
  if any(index(:))
    y(index) = atan_pos_small(x(index),rnd);
  end

function y = atan_pos_small(x,rnd)
% Rigorous atan(x) for nonnegative double vector x with  0 <= x < 4
% rounding corresponding to rnd

  INTLAB_STDFCTS_ATAN = getappdata(0,'INTLAB_STDFCTS_ATAN');

  setround(0)                               % maximum first 15 bits, no bit
  xs = pow2( floor(pow2(x,13)) , -13 );     %   below 2^-13
  d = x - xs;                               % 0 <= d < 2^-15*x
  atanxs = atan(xs);

  % atan(xs+d) = atan(xs) + atan(d/(1+x*xs))
  setround(-rnd)
  E = 1 + x.*xs;
  setround(rnd)
  E = d./E;                                 % 0 <= E < d/(1+x*xs) < 2^-13

  if rnd==-1                         % 0 <= err <= E^5/5 <= d^4/5*d < 5e-17*d
    % atanE <= atan(E)
    atanE = E + (((-E).*E).*E)/3;
    y = atanxs + ( atanE + (-INTLAB_STDFCTS_ATAN.EPS)*atanxs );
 else
    % atanE >= atan(E)
    atanE = ((( E.*E/5 + (-1)/3 ).*E).*E).*E + E;
    y = atanxs + ( atanE + INTLAB_STDFCTS_ATAN.EPS*atanxs );
  end

repo: OpenMAS


text:
function outstr = GE_dumpExamHeader(exam)
%
% outstr = GE_dumpExamHeader(exam)
%
% Writes the exam header to the string outstr
%
% Souheil J. Inati
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%

outstr = sprintf('ex_suid: %s\n'     , char(exam.ex_suid));
outstr = strcat(outstr, sprintf('ex_uniq: %d\n'     , exam.ex_uniq));
outstr = strcat(outstr, sprintf('ex_diskid: %s\n'   , char(exam.ex_diskid)));
outstr = strcat(outstr, sprintf('ex_no: %d\n'       , exam.ex_no));
outstr = strcat(outstr, sprintf('hospname: %s\n'    , cstring(exam.hospname)));
outstr = strcat(outstr, sprintf('detect: %d\n'      , exam.detect));
outstr = strcat(outstr, sprintf('numcells: %d\n'    , exam.numcells));
outstr = strcat(outstr, sprintf('zerocell: %f\n'    , exam.zerocell));
outstr = strcat(outstr, sprintf('cellspace: %f\n'   , exam.cellspace));
outstr = strcat(outstr, sprintf('srctodet: %f\n'    , exam.srctodet));
outstr = strcat(outstr, sprintf('srctoiso: %f\n'    , exam.srctoiso));
outstr = strcat(outstr, sprintf('tubetyp: %d\n'     , exam.tubetyp));
outstr = strcat(outstr, sprintf('dastyp: %d\n'      , exam.dastyp));
outstr = strcat(outstr, sprintf('num_dcnk: %d\n'    , exam.num_dcnk));
outstr = strcat(outstr, sprintf('dcn_len: %d\n'     , exam.dcn_len));
outstr = strcat(outstr, sprintf('dcn_density: %d\n' , exam.dcn_density));
outstr = strcat(outstr, sprintf('dcn_stepsize: %d\n', exam.dcn_stepsize));
outstr = strcat(outstr, sprintf('dcn_shiftcnt: %d\n', exam.dcn_shiftcnt));
outstr = strcat(outstr, sprintf('magstrength: %d\n' , exam.magstrength));
outstr = strcat(outstr, sprintf('patid: %s\n'       , cstring(exam.patid)));
outstr = strcat(outstr, sprintf('patname: %s\n'     , cstring(exam.patname)));
outstr = strcat(outstr, sprintf('patage: %d\n'      , exam.patage));
outstr = strcat(outstr, sprintf('patian: %d\n'      , exam.patian));
outstr = strcat(outstr, sprintf('patsex: %d\n'      , exam.patsex));
outstr = strcat(outstr, sprintf('patweight: %d\n'   , exam.patweight));
outstr = strcat(outstr, sprintf('trauma: %d\n'      , exam.trauma));
outstr = strcat(outstr, sprintf('hist: %s\n'        , cstring(exam.hist)));
outstr = strcat(outstr, sprintf('reqnum: %s\n'      , cstring(exam.reqnum)));
outstr = strcat(outstr, sprintf('ex_datetime: %d\n' , exam.ex_datetime));
outstr = strcat(outstr, sprintf('refphy: %s\n'      , cstring(exam.refphy)));
outstr = strcat(outstr, sprintf('diagrad: %s\n'     , cstring(exam.diagrad)));
outstr = strcat(outstr, sprintf('op: %s\n'          , cstring(exam.op)));
outstr = strcat(outstr, sprintf('ex_desc: %s\n'     , cstring(exam.ex_desc)));
outstr = strcat(outstr, sprintf('ex_typ: %s\n'      , cstring(exam.ex_typ)));
outstr = strcat(outstr, sprintf('ex_format: %d\n'   , exam.ex_format));
outstr = strcat(outstr, sprintf('firstaxtime: %f\n' , exam.firstaxtime));
outstr = strcat(outstr, sprintf('ex_sysid: %s\n'    , sprintf('%s',exam.ex_sysid)));
outstr = strcat(outstr, sprintf('ex_lastmod: %d\n'  , exam.ex_lastmod));
outstr = strcat(outstr, sprintf('protocolflag: %d\n', exam.protocolflag));
outstr = strcat(outstr, sprintf('ex_alloc_key: %s\n', sprintf('%s',exam.ex_alloc_key)));
outstr = strcat(outstr, sprintf('ex_delta_cnt: %d\n', exam.ex_delta_cnt));
outstr = strcat(outstr, sprintf('ex_verscre: %s\n'  , char(exam.ex_verscre)));
outstr = strcat(outstr, sprintf('ex_verscur: %s\n'  , char(exam.ex_verscur)));
outstr = strcat(outstr, sprintf('ex_checksum: %d\n' , exam.ex_checksum));
outstr = strcat(outstr, sprintf('ex_complete: %d\n' , exam.ex_complete));
outstr = strcat(outstr, sprintf('ex_seriesct: %d\n' , exam.ex_seriesct));
outstr = strcat(outstr, sprintf('ex_numarch: %d\n'  , exam.ex_numarch));
outstr = strcat(outstr, sprintf('ex_numseries: %d\n', exam.ex_numseries));
outstr = strcat(outstr, sprintf('ex_series.length: %d\n'   , exam.ex_series.length));
outstr = strcat(outstr, sprintf('ex_series.data: %d\n'   , exam.ex_series.data));
outstr = strcat(outstr, sprintf('ex_numunser: %d\n' , exam.ex_numunser));
outstr = strcat(outstr, sprintf('ex_unseries.length: %d\n' , exam.ex_unseries.length));
outstr = strcat(outstr, sprintf('ex_unseries.data: %d\n' , exam.ex_unseries.data));
outstr = strcat(outstr, sprintf('ex_toarchcnt: %d\n', exam.ex_toarchcnt));
outstr = strcat(outstr, sprintf('ex_toarchive.length: %d\n', exam.ex_toarchive.length));
outstr = strcat(outstr, sprintf('ex_toarchive.data: %d\n', exam.ex_toarchive.data));
outstr = strcat(outstr, sprintf('ex_prospcnt: %d\n' , exam.ex_prospcnt));
outstr = strcat(outstr, sprintf('ex_prosp.length: %d\n'    , exam.ex_prosp.length));
outstr = strcat(outstr, sprintf('ex_prosp.data: %d\n'    , exam.ex_prosp.data));
outstr = strcat(outstr, sprintf('ex_modelnum: %d\n' , exam.ex_modelnum));
outstr = strcat(outstr, sprintf('ex_modelcnt: %d\n' , exam.ex_modelcnt));
outstr = strcat(outstr, sprintf('ex_models.length: %d\n'   , exam.ex_models.length));
outstr = strcat(outstr, sprintf('ex_models.data: %d\n'   , exam.ex_models.data));
outstr = strcat(outstr, sprintf('ex_stat: %d\n'     , exam.ex_stat));
outstr = strcat(outstr, sprintf('uniq_sys_id: %s\n' , sprintf('%s',exam.uniq_sys_id)));
outstr = strcat(outstr, sprintf('service_id: %s\n'  , sprintf('%s',exam.service_id)));
outstr = strcat(outstr, sprintf('mobile_loc: %s\n'  , sprintf('%s',exam.mobile_loc)));
outstr = strcat(outstr, sprintf('study_uid: %s\n'   , sprintf('%s',exam.study_uid)));
outstr = strcat(outstr, sprintf('study_status: %d\n', exam.study_status));
outstr = strcat(outstr, sprintf('ex_padding: %s\n' , cstring(exam.ex_padding)));

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% SELECT.M          (universal SELECTion)
%
% This function performs universal selection. The function handles
% multiple populations and calls the low level selection function
% for the actual selection process.

%
% Syntax:  SelCh = select(SEL_F, Chrom, FitnV, GGAP, SUBPOP)
%
% Input parameters:
%    SEL_F     - Name of the selection function
%    Chrom     - Matrix containing the individuals (parents) of the current
%                population. Each row corresponds to one individual.
%    FitnV     - Column vector containing the fitness values of the
%                individuals in the population.
%    GGAP      - (optional) Rate of individuals to be selected
%                if omitted 1.0 is assumed
%    SUBPOP    - (optional) Number of subpopulations
%                if omitted 1 subpopulation is assumed
%
% Output parameters:
%    SelCh     - Matrix containing the selected individuals.

% Author:     Hartmut Pohlheim
% History:    10.03.94     file created

function SelCh = select(SEL_F, Chrom, FitnV, GGAP, SUBPOP);

% Check parameter consistency
   if nargin < 3, error('Not enough input parameter'); end

   % Identify the population size (Nind)
   [NindCh,Nvar] = size(Chrom);
   [NindF,VarF] = size(FitnV);
   if NindCh ~= NindF, error('Chrom and FitnV disagree'); end
   if VarF ~= 1, error('FitnV must be a column vector'); end

   if nargin < 5, SUBPOP = 1; end
   if nargin > 4,
      if isempty(SUBPOP), SUBPOP = 1;
      elseif isnan(SUBPOP), SUBPOP = 1;
      elseif length(SUBPOP) ~= 1, error('SUBPOP must be a scalar'); end
   end

   if (NindCh/SUBPOP) ~= fix(NindCh/SUBPOP), error('Chrom and SUBPOP disagree'); end
   Nind = NindCh/SUBPOP;  % Compute number of individuals per subpopulation

   if nargin < 4, GGAP = 1; end
   if nargin > 3,
      if isempty(GGAP), GGAP = 1;
      elseif isnan(GGAP), GGAP = 1;
      elseif length(GGAP) ~= 1, error('GGAP must be a scalar');
      elseif (GGAP < 0), error('GGAP must be a scalar bigger than 0'); end
   end

% Compute number of new individuals (to select)
   NSel=max(floor(Nind*GGAP+.5),2);

% Select individuals from population
   SelCh = [];
   for irun = 1:SUBPOP,
      FitnVSub = FitnV((irun-1)*Nind+1:irun*Nind);
      ChrIx=feval(SEL_F, FitnVSub, NSel)+(irun-1)*Nind;
      SelCh=[SelCh; Chrom(ChrIx,:)];
   end

% End of function

author: vonsylvia


text:
function [a] = save_warper_3(string, proposal,mask,class,proposal_iou,proposal_pole_mask,response_index,original_index,mask_remained,class_remained,original_index_remained,all_class_flag)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

 a=0;
 save(string,'proposal','mask','class','proposal_iou','proposal_pole_mask','response_index','original_index','mask_remained','class_remained','original_index_remained','all_class_flag');
end

author: xjqicuhk

repo: SIMS

save_path: github-repos/MATLAB/xjqicuhk-SIMS


text:
function featureS = calcRadiomicsForImgType(volOrig3M,maskBoundingBox3M,paramS,gridS)
%calcRadiomicsForImgType.m
%Derive user-defined image type and extract radiomics features.
%
%AI 3/28/19
%AI 5/1/19    Turned off flag for diffAvg since this is equivalent to dissimilarity.

% Voxel volume for Total Energy calculation
xValsV = gridS.xValsV;
yValsV = gridS.yValsV;
zValsV = gridS.zValsV;
PixelSpacingX = gridS.PixelSpacingV(1);
PixelSpacingY = gridS.PixelSpacingV(2);
PixelSpacingZ = gridS.PixelSpacingV(3);
VoxelVol = PixelSpacingX*PixelSpacingY*PixelSpacingZ*1000; % convert cm to mm

whichFeatS = paramS.whichFeatS;
featureS = struct;

% Get image types with various parameters
fieldNamC = fieldnames(paramS.imageType);
imageTypeC = {};
for iImg = 1:length(fieldNamC)
    for iFilt = 1:length(paramS.imageType.(fieldNamC{iImg}))
        filtParamS = struct();
        filtParamS.imageType = fieldNamC{iImg};
        filtParamS.paramS = paramS.imageType.(fieldNamC{iImg})(iFilt);
        imageTypeC{end+1} = filtParamS;
    end
end

% ---- Calc. shape features (same across img. types) ----
tic
if whichFeatS.shape.flag
    rcsV = [];
    if isfield(paramS,'shapeParamS') && isfield(paramS.shapeParamS,'rcs')
        rcsV = paramS.shapeParamS.rcs.';
    end
    featureS.shapeS = getShapeParams(maskBoundingBox3M, ...
        {xValsV, yValsV, zValsV},rcsV);
end
toc

%% Loop over image types
maskOrig3M = maskBoundingBox3M;
for k = 1:length(imageTypeC)

    %Generate volume based on original/derived imageType
    if strcmpi(imageTypeC{k}.imageType,'original')
        if isfield(paramS,'toQuantizeFlag')
            quantizeFlag = paramS.toQuantizeFlag;
        else
            quantizeFlag = 0;
        end
        minClipIntensity = [];
        maxClipIntensity = [];
        if isfield(paramS,'textureParamS') && ...
                isfield(paramS.textureParamS,'minClipIntensity')
            minClipIntensity = paramS.textureParamS.minClipIntensity;
        end
        if isfield(paramS,'textureParamS') && ...
                isfield(paramS.textureParamS,'maxClipIntensity')
            maxClipIntensity = paramS.textureParamS.maxClipIntensity;
        end
        volToEval = volOrig3M;
        [minr, maxr, minc, maxc, mins, maxs] = compute_boundingbox(maskBoundingBox3M);
        maskBoundingBox3M = maskOrig3M(minr:maxr, minc:maxc, mins:maxs);
        volToEval = volToEval(minr:maxr, minc:maxc, mins:maxs);
    else
        %Add voxel size in mm to paramS
        voxSizV = [PixelSpacingX, PixelSpacingY, PixelSpacingZ]*10; %convert cm to mm
        imageTypeC{k}.paramS.VoxelSize_mm.val = voxSizV;
        if paramS.whichFeatS.padding.flag && ...
                ~strcmpi(whichFeatS.padding.method,'none')
            paddingSizeV = paramS.whichFeatS.padding.size;
            if length(paddingSizeV)==2
                paddingSizeV(end+1) = 0;
            end
            imageTypeC{k}.paramS.padding.size = paddingSizeV;
        else
            imageTypeC{k}.paramS.padding.size = [0,0,0];
        end

        outS = processImage(imageTypeC{k}.imageType,volOrig3M,maskOrig3M,...
            imageTypeC{k}.paramS);
        [minr, maxr, minc, maxc, mins, maxs] = compute_boundingbox(maskOrig3M);
        maskBoundingBox3M = maskOrig3M(minr:maxr, minc:maxc, mins:maxs);
        derivedImgName = fieldnames(outS);
        volToEval = outS.(derivedImgName{1});
        volToEval = volToEval(minr:maxr, minc:maxc, mins:maxs);
        if whichFeatS.texture.flag
            % always quantize the derived image for texture calc.
            quantizeFlag = true;
        else
            if isfield(paramS,'toQuantizeFlag')
                quantizeFlag = paramS.toQuantizeFlag;
            else
                quantizeFlag = false;
            end
        end
        minClipIntensity = []; % no clipping imposed for derived images
        maxClipIntensity = [];
        if isfield(paramS,'textureParamS') && ...
                isfield(paramS.textureParamS,'minClipIntensity')
            minClipIntensity = paramS.textureParamS.minClipIntensity;
        end
        if isfield(paramS,'textureParamS') && ...
                isfield(paramS.textureParamS,'maxClipIntensity')
            maxClipIntensity = paramS.textureParamS.maxClipIntensity;
        end
    end

    % Volume without NaNs for Peak/Valley computation
    volToEvalOrigM = volToEval;

    % Quantize the volume of interest
    if quantizeFlag
        numGrLevels = [];
        binwidth = [];

        if isfield(paramS,'textureParamS')
            if isfield(paramS.textureParamS,'numGrLevels')
                numGrLevels = paramS.textureParamS.numGrLevels;
            end
            if isfield(imageTypeC{k}.paramS,'textureParamS') && ...
                    isfield(imageTypeC{k}.paramS.textureParamS,'numGrLevels')
                numGrLevels = imageTypeC{k}.paramS.textureParamS.numGrLevels;
            end
            if isfield(paramS.textureParamS,'binwidth')
                binwidth = paramS.textureParamS.binwidth;
            end
            if isfield(imageTypeC{k}.paramS,'textureParamS') &&...
                    isfield(imageTypeC{k}.paramS.textureParamS,'binwidth')
                binwidth = imageTypeC{k}.paramS.textureParamS.binwidth;
            end
        end

        % Don't use intensities outside the ROI in discretization
        volToEval(~maskBoundingBox3M) = NaN;
        quantizedM = imquantize_cerr(volToEval,numGrLevels,...
            minClipIntensity,maxClipIntensity,binwidth);
        % Reassign the number of gray levels in case they were computed for the
        % passed binwidth
        numGrLevels = max(quantizedM(:));
        %paramS.textureParamS.numGrLevels = numGrLevels;

    else
        quantizedM = volToEval;
    end

    quantizedM(~maskBoundingBox3M) = NaN;
    numVoxels = sum(~isnan(quantizedM(:)));

    %Feature calculation
    outFieldName = createFieldNameFromParameters...
        (imageTypeC{k}.imageType,imageTypeC{k}.paramS);

    % --- 1. First-order features ---
    if whichFeatS.firstOrder.flag
        binWidthEntropy = [];
        binNumEntropy = [];
        offsetForEnergy = 0;
        if isfield(paramS.firstOrderParamS,'offsetForEnergy')
            offsetForEnergy = paramS.firstOrderParamS.offsetForEnergy;
        end
        if isfield(paramS.firstOrderParamS,'binWidthEntropy') && ...
            ~isempty(paramS.firstOrderParamS.offsetForEnergy)
            binWidthEntropy = paramS.firstOrderParamS.binWidthEntropy;
        end
        if isfield(paramS.firstOrderParamS,'binNumEntropy') && ...
                ~isempty(paramS.firstOrderParamS.binNumEntropy)
            binNumEntropy = paramS.firstOrderParamS.binNumEntropy;
        end
        if isfield(imageTypeC{k}.paramS,'firstOrderParamS') && ...
                isfield(imageTypeC{k}.paramS.firstOrderParamS,'offsetForEnergy')
            offsetForEnergy = imageTypeC{k}.paramS.firstOrderParamS.offsetForEnergy;
        end
        if isfield(imageTypeC{k}.paramS,'binWidthEntropy') && ...
                isfield(imageTypeC{k}.paramS.firstOrderParamS,'binWidthEntropy')
            binWidthEntropy = imageTypeC{k}.paramS.firstOrderParamS.binWidthEntropy;
        end
        if isfield(imageTypeC{k}.paramS,'binNumEntropy') && ...
                isfield(imageTypeC{k}.paramS.firstOrderParamS,'binNumEntropy')
            binNumEntropy = imageTypeC{k}.paramS.firstOrderParamS.binNumEntropy;
        end
        volV = volToEval(logical(maskBoundingBox3M));
        volV = volV(:);
        featureS.(outFieldName).firstOrderS = radiomics_first_order_stats...
            (volV, VoxelVol,offsetForEnergy,binWidthEntropy,binNumEntropy);
    end

    %---2. Higher-order (texture) features ----

    %Get directionality and avg type
    if any([whichFeatS.glcm.flag,whichFeatS.glrlm.flag,whichFeatS.gtdm.flag,...
            whichFeatS.gldm.flag,whichFeatS.glszm.flag])

        directionality = paramS.textureParamS.directionality;
        avgType = paramS.textureParamS.avgType;
        switch lower(directionality)
            case '2d'
                dirctn = 2;
            case '3d'
                dirctn = 1;
            otherwise
                error('Invalid input. Directionality must be "2D" or "3D"');
        end

        switch lower(avgType)
            case 'texturematrix'
                %Haralick features with combined cooccurrence matrix
                cooccurType = 1;
            case 'feature'
                %'Haralick features from separate cooccurrence matrix per direction, averaged'
                cooccurType = 2;
            otherwise
                error('Invalid input. Directionality must be "2D" or "3D"');
        end

        %numGrLevels = paramS.textureParamS.numGrLevels;
        voxelOffset = paramS.textureParamS.voxelOffset;

        % a. GLCM
        if whichFeatS.glcm.flag

            featC = whichFeatS.glcm.featureList;
            glcmFlagS = getHaralickFlags(featC);
            featureS.(outFieldName).glcmFeatS = get_haralick(dirctn, voxelOffset, cooccurType, quantizedM, ...
                numGrLevels, glcmFlagS);

        end

        % b. GLRLM
        if whichFeatS.glrlm.flag
            featC = whichFeatS.glrlm.featureList;
            rlmFlagS = getRunLengthFlags(featC);
            rlmType = cooccurType;
            featureS.(outFieldName).rlmFeatS = get_rlm(dirctn, rlmType, quantizedM, ...
                numGrLevels, numVoxels, rlmFlagS);
        end

        %c. GTDM
        if whichFeatS.gldm.flag
            patchRadiusV = paramS.textureParamS.patchRadiusVox;
            [s,p] = calcNGTDM(quantizedM, patchRadiusV, ...
                numGrLevels);
            featureS.(outFieldName).ngtdmFeatS = ngtdmToScalarFeatures(s,p,numVoxels);
        end

        %d. GLDM
        if whichFeatS.gldm.flag
            patchRadiusV = paramS.textureParamS.patchRadiusVox;
            imgDiffThresh = paramS.textureParamS.imgDiffThresh;
            ngldM = calcNGLDM(quantizedM, patchRadiusV,numGrLevels,imgDiffThresh);
            featureS.(outFieldName).ngldmFeatS = ngldmToScalarFeatures(ngldM,numVoxels);
        end

        %e. GLSZM
        if whichFeatS.glszm.flag
            featC = whichFeatS.glszm.featureList;
            szmFlagS = getSizeZoneFlags(featC);
            szmType = dirctn; % 1: 3d, 2: 2d
            szmM = calcSZM(quantizedM, numGrLevels, szmType);
            numVoxels = sum(~isnan(quantizedM(:)));
            featureS.(outFieldName).szmFeatS = szmToScalarFeatures(szmM,numVoxels, szmFlagS);
        end

        %f. Peak-valley
        if whichFeatS.peakValley.flag
            radiusV = paramS.peakValleyParamS.peakRadius;
            units = paramS.peakValleyParamS.units; %'cm' or 'vox'
            featureS.(outFieldName).peakValleyFeatureS = getImPeakValley(maskBoundingBox3M,...
                volToEvalOrigM, radiusV, units);
        end

        %g. IVH
        if whichFeatS.ivh.flag
            IVHBinWidth = paramS.ivhParamS.binwidth; %IVH binwidth
            xForIxV = paramS.ivhParamS.xForIxPct; % percentage volume
            xAbsForIxV = paramS.ivhParamS.xForIxCc; % absolute volume [cc]
            xForVxV = paramS.ivhParamS.xForVxPct; % percent intensity cutoff
            xAbsForVxV = paramS.ivhParamS.xForVxAbs; % absolute intensity cutoff [HU]
            scanV = double(volToEval(maskBoundingBox3M));
            volV = repmat(VoxelVol,numel(scanV),1);
            featureS.(outFieldName).ivhFeaturesS = getIvhParams(scanV, volV, IVHBinWidth,...
                xForIxV, xAbsForIxV, xForVxV, xAbsForVxV);

        end
    end

end

%% -------------- Sub-functions ----------------------------
    function glcmFlagS = getHaralickFlags(varargin)

        %Default: all features
        glcmFlagS.energy = 1;
        glcmFlagS.jointEntropy = 1;
        glcmFlagS.jointMax = 1;
        glcmFlagS.jointAvg = 1;
        glcmFlagS.jointVar = 1;
        glcmFlagS.contrast = 1;
        glcmFlagS.invDiffMoment = 1;
        glcmFlagS.sumAvg = 1;
        glcmFlagS.corr = 1;
        glcmFlagS.clustShade = 1;
        glcmFlagS.clustProm = 1;
        glcmFlagS.haralickCorr = 1;
        glcmFlagS.invDiffMomNorm = 1;
        glcmFlagS.invDiff = 1;
        glcmFlagS.invDiffNorm = 1;
        glcmFlagS.invVar = 1;
        glcmFlagS.dissimilarity = 1;
        glcmFlagS.diffEntropy = 1;
        glcmFlagS.diffVar = 1;
        glcmFlagS.diffAvg = 0;  %Equivalent to dissimilarity
        glcmFlagS.sumVar = 1;
        glcmFlagS.sumEntropy = 1;
        glcmFlagS.clustTendency = 1;
        glcmFlagS.autoCorr = 1;
        glcmFlagS.invDiffMomNorm = 1;
        glcmFlagS.firstInfCorr = 1;
        glcmFlagS.secondInfCorr = 1;

        if nargin==1 && ~strcmpi(varargin{1},'all')
            featureC = varargin{1};
            allFeatC = fieldnames(glcmFlagS);
            idxV = find(~ismember(allFeatC,featureC));
            for n = 1:length(idxV)
                glcmFlagS.(allFeatC{idxV(n)})  = 0;
            end
        end

    end

    function rlmFlagS = getRunLengthFlags(varargin)

        rlmFlagS.shortRunEmphasis = 1;
        rlmFlagS.longRunEmphasis = 1;
        rlmFlagS.grayLevelNonUniformity = 1;
        rlmFlagS.grayLevelNonUniformityNorm = 1;
        rlmFlagS.runLengthNonUniformity = 1;
        rlmFlagS.runLengthNonUniformityNorm = 1;
        rlmFlagS.runPercentage = 1;
        rlmFlagS.lowGrayLevelRunEmphasis = 1;
        rlmFlagS.highGrayLevelRunEmphasis = 1;
        rlmFlagS.shortRunLowGrayLevelEmphasis = 1;
        rlmFlagS.shortRunHighGrayLevelEmphasis = 1;
        rlmFlagS.longRunLowGrayLevelEmphasis = 1;
        rlmFlagS.longRunHighGrayLevelEmphasis = 1;
        rlmFlagS.grayLevelVariance = 1;
        rlmFlagS.runLengthVariance = 1;
        rlmFlagS.runEntropy = 1;

        if nargin==1 && ~strcmpi(varargin{1},'all')
            featureC = lower(varargin{1});
            allFeatC = lower(fieldnames(rlmFlagS));
            idxV = find(~ismember(allFeatC,featureC));
            for n = 1:length(idxV)
                rlmFlagS.(allFeatC{idxV(n)})  = 0;
            end
        end

    end

function szmFlagS = getSizeZoneFlags(varargin)

        szmFlagS.grayLevelNonUniformity = 1;
        szmFlagS.grayLevelNonUniformityNorm = 1;
        szmFlagS.grayLevelVariance = 1;
        szmFlagS.highGrayLevelZoneEmphasis = 1;
        szmFlagS.lowGrayLevelZoneEmphasis = 1;
        szmFlagS.largeAreaEmphasis = 1;
        szmFlagS.largeAreaHighGrayLevelEmphasis = 1;
        szmFlagS.largeAreaLowGrayLevelEmphasis = 1;
        szmFlagS.sizeZoneNonUniformity = 1;
        szmFlagS.sizeZoneNonUniformityNorm = 1;
        szmFlagS.sizeZoneVariance = 1;
        szmFlagS.zonePercentage = 1;
        szmFlagS.smallAreaEmphasis = 1;
        szmFlagS.smallAreaLowGrayLevelEmphasis = 1;
        szmFlagS.smallAreaHighGrayLevelEmphasis = 1;
        szmFlagS.zoneEntropy = 1;

        if nargin==1 && ~strcmpi(varargin{1},'all')
            featureC = lower(varargin{1});
            allFeatC = fieldnames(szmFlagS);
            idxV = find(~ismember(allFeatC,featureC));
            for n = 1:length(idxV)
                szmFlagS.(allFeatC{idxV(n)})  = 0;
            end
        end

    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR