text:
function [ fighandle, surfh ] = DrawImageSlice3D( dicom_filename,  fighandle,  transparency )
%============================
% [ fighandle, surfh ] = DrawImageSlice3D( dicom_filename,  fighandle,  transparency )
% This is a function for drawing slices in 3D space based on the dicom header
% information.
% Input:
%   dicom_filename = full path to a DICOM image file. Default = '*'
%                   empty string ('') or and string with a wildcard (*) will prompt the file browser
%   figh = figure handle
%   transparency = alpha for the image slice.
%                   Use values < 1 to get partial transparency when
%                   plotting multiple slices in one figure
% Example:
% First pick and image (assuming the DICOM file extension is IMA) and draw the slice semi-transparent
%  [ fighandle, surfh(1) ] = DrawImageSlice3D( '*.IMA',  [],  0.5 );
% Next, call up a second image and draw a slice in the same figure
%  [ fighandle, surfh(2) ] = DrawImageSlice3D( 'mysecond_dicomim.IMA', fighandle,  0.5 );
%
% The axes are rotated at the end of this function.
% A simple edit setting 'show_rotate' to 0 will suppress this.
%
% NOTE: In the DICOM file there are several series UID and reference frame UID to
% ensure that the slices are from the same study and position.
% No check for this is included in this function. One could compose figures
% with slices of different subjects if so inclined.
% The DICOM fields used should be universal, the reference frame should be
% patient orientation dependend X = RL, Y =AP, Z = FH
%
%======================================================================

%======================================================================
% Ronald Ouwerkerk , NIH/NIDDK, June 2010
% Tested on Mac OS 10.5 with Martlab  2008b on Siemens DICOM images (VB15 export to offline)
% Images were all head first supine. Let me know if other orientations
% create problems.
%======================================================================

%% Define constants
topl = 1;
topr = 2;
botl = 3;
botr = 4;

X = 1;
Y = 2;
Z = 3;

%% Set defaults for input arguments
if nargin <1
    dicom_filename = '*';
end

if nargin <2
    fighandle = [];
end

if ~ishandle( fighandle )
    msgstr = sprintf('Warning in %s: Second argument is not a valid figure handle. Creating a new one', mfilename);
    disp( msgstr)
     fighandle = [];
end

if nargin <3
   transparency = 1;
end

if isempty( dicom_filename )
    dicom_filename = '*';
end

%save the current path
orig_path = pwd;

% Call the filebrowswer UI if the dicom filename contains a wildcard
if ~isempty( findstr( dicom_filename, '*') )
    [ filename, fpath ] = uigetfile( dicom_filename , 'Select DICOM file');
    cd( fpath );
    dicom_filename =filename;
end

%% Read DICOM header and image data
dinfo= dicominfo( dicom_filename );
imdata = dicomread(  dicom_filename );
cd( orig_path )

%% Calculate slice corner positions from the DICOM header info
% Get the top left corner position in XYZ coordinates
pos    = dinfo.ImagePositionPatient;

nc = double(dinfo.Columns);
nr = double(dinfo.Rows);

% Get the row and column direction vercors in XYZ coordinates
row_dircos(X) = dinfo.ImageOrientationPatient(1);
row_dircos(Y) = dinfo.ImageOrientationPatient(2);
row_dircos(Z) = dinfo.ImageOrientationPatient(3);
col_dircos(X) = dinfo.ImageOrientationPatient(4);
col_dircos(Y) = dinfo.ImageOrientationPatient(5);
col_dircos(Z) = dinfo.ImageOrientationPatient(6);

% % Check normality and orthogonality of the row and col vectors
% Crownorm = dot(row_dircos, row_dircos);
% Ccolnorm = dot(col_dircos, col_dircos);
% Cdotprod = dot(row_dircos, col_dircos);
%
% if abs(Cdotprod) > 1e-5
%     warnstr = sprintf('Possible dicominfo error: the dotproduct of the row and col vectors is %f should be 0',Cdotprod );
%     disp(warnstr)
% end

% Calculate image dimensions
row_length = dinfo.PixelSpacing(1) * nr;
col_length = dinfo.PixelSpacing(2) * nc;

%% Set up the corner positions matrix in XYZ coordinates
% Top Left Hand Corner
corners( topl, X) = pos(X);
corners( topl, Y) = pos(Y);
corners( topl, Z) = pos(Z);

% Top Right Hand Corner
corners( topr, X) = pos(X) + row_dircos(X) * row_length;
corners( topr, Y) = pos(Y) + row_dircos(Y) * row_length;
corners( topr, Z) = pos(Z) + row_dircos(Z) * row_length;

% Bottom Left Hand Corner
corners( botl, X) = pos(X) + col_dircos(X) * col_length;
corners( botl, Y) = pos(Y) + col_dircos(Y) * col_length;
corners( botl, Z) = pos(Z) + col_dircos(Z) * col_length;

% Bottom Right Hand Corner
corners( botr, X) = pos(X) + row_dircos(X) * row_length + col_dircos(X) * col_length;
corners( botr, Y) = pos(Y) + row_dircos(Y) * row_length + col_dircos(Y) * col_length;
corners( botr, Z) = pos(Z) + row_dircos(Z) * row_length + col_dircos(Z) * col_length;

%% Prepare the figure
% Select active figure, set hold on to alow multiple slices in one figure
if isempty(  fighandle )
     fighandle = figure;
     colormap( gray );
     newfig = 1;
else
    newfig = 0;
end

figure( fighandle );
hold on;

%Tidy up the figure
% set aspect ratio
daspect( [1,1,1]);
set( gca, 'color', 'none')

%% Display slice
%  normalize image data
imdata = double( imdata );
imdata = imdata / max( imdata(:));
% scale the image
I = imdata*255;
% create an alternative matrix for corner points
A( 1,1 , 1:3 ) = corners( topl, : );
A( 1,2 , 1:3 ) = corners( topr, : );
A( 2,1 , 1:3 ) = corners( botl, : );
A( 2,2 , 1:3 ) = corners( botr, : );
% extract the coordinates for the surfaces
x = A( :,:,X );
y = A( :,:,Y );
z = A( :,:,Z );

% plot surface
surfh = surface('XData',x,'YData',y,'ZData',z,...
'CData', I,...
'FaceColor','texturemap',...
'EdgeColor','none',...
'LineStyle','none',...
'Marker','none',...
'MarkerFaceColor','none',...
'MarkerEdgeColor','none',...
'CDataMapping','direct');
%set transparency level
set( surfh, 'FaceAlpha', transparency );
% label axes and optimize figure
xlabel('RL');
ylabel('AP');
zlabel('FH');

% if only one slice is in the figure this may flatten the 3rd axis
if ~newfig
    axis tight
end

%% Optional: rotate to show all
% edit this to create a movie
do_movie = 0;
show_rotate = 1

if do_movie
    % open avifile
    aviobj = avifile('slice3Drotate.avi');
    show_rotate = 1;
end

if show_rotate
    % tilt and rotate
    el = 22;
    for azplus = 0:10:360
        az = mod( 45+azplus, 360);
        view( [az, el] );
        if do_movie
            % add farem to the movie
            frame = getframe(fh);
            aviobj = addframe(aviobj,frame);
        else
            % needed to see the intemediate steps
           drawnow;
        end
    end
end

if do_movie
    % close avifile
    close( aviobj );
end

return


text:
function [modelIDsMissingInTable,incorrectIDsInTable] = testInputFiles
% Tests if the names of all reconstructed microbes are present in the input
% files for the AGORA2 pipeline.
%
% OUTPUT
% modelIDsMissingInTable       AGORA2 reconstruction IDs missing in the table
% incorrectIDsInTable          Reconstruction IDs in the table not
%                              corresponding to any reconstruction
% Almut Heinken, 09/2019

[~, infoFile, ~] = xlsread('AGORA2_infoFile.xlsx');
models=infoFile(2:end,1);

inputTables = {
    'AromaticAATable'
    'BileAcidTable'
    'CarbonSourcesTable'
    'FermentationTable'
    'GenomeAnnotation'
    'GrowthRequirementsTable'
    'PutrefactionTable'
    'uptakeTable'
    'secretionProductTable'
    };

for i = 1:length(inputTables)
    table = readtable([inputTables{i} '.txt'], 'ReadVariableNames', false, 'Delimiter', 'tab','TreatAsEmpty',['UND. -60001','UND. -2011','UND. -62011']);
    table = table2cell(table);
    notintable = setdiff(models,table(2:end,1));
    modelIDsMissingInTable.(inputTables{i})= notintable;
    incorrectIDs=setdiff(table(2:end,1),models);
    incorrectIDsInTable.(inputTables{i})= incorrectIDs;
    if ~isempty(notintable)
        warning(['Some reconstruction names are not present in table ' inputTables{i}])
    end
    if ~isempty(incorrectIDs)
        warning(['Some reconstruction IDs are incorrect in table ' inputTables{i}])
    end
end

end

author: opencobra

repo: cobratoolbox


text:
function [ output_args ] = test_gmm_ex( input_args )
%TEST_GMM_HU_EX Summary of this function goes here
%   Detailed explanation goes here
close all;

dataset = '7';
[endmembers,I,Y,R_gt,A_gt,names,wl] = prepare_supervised_unmixing(dataset);
[rows,cols,B] = size(I);

% noise covariance matrix. it can be zero.
D = 0.001^2 * eye(B);

% smoothness and sparsity constraints on the abundances
options.beta1 = 0;
options.beta2 = 0;

% show intermediate results (1) or not (0)
options.show_fig = 1;

options.names = names;
options.D = D;

% project_mode can be
%   'image' - apply PCA on the pixels of the image
%   'endmembers' - concatenate the spectra in the library and use PCA on them
options.project_mode = 'image';

% threshold of convergence
options.convergence_thresh = 0.0001;

% calculate abundance maps
[A,R,w_jk,mu_jk,sigma_jk,extra] = gmm_hu_ex(I, endmembers, options);

% calculate endmembers per pixel
E = gmm_hu_endmember(I,A,D,w_jk,mu_jk,sigma_jk);

save('result_gmm.mat','A','R','E','w_jk','mu_jk','sigma_jk');

mdiff(A,A_gt);
show_abundances(A,rows,cols);
replay_scatter_abund(extra.frames_scatter, extra.frames_abund);

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [r1,r2] = quad_roots(a1, a2, a3)

t1 = -a2/2./a1;
t2 = sqrt(a2.^2 - 4*a1.*a3)/2./a1;
r1 = t1 + t2;
r2 = t1 - t2;

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


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

%%
clear; close all; clc;

%% Syntax
% |[Fs,Vs,Cs,Cs_rgb,Cs_d]=polyTube(Vg,optStruct);|

%% Description
% The polyTube function can be used to
%% Examples

%%
% Plot settings
markerSize=15;
lineWidth=2;
fontSize=10;

%% Example 1: Creating a swept tube

%%
% Creating example curve

t=linspace(0,4*pi,50);
Vg=[ t(:) sin(t(:)) cos(t(:))];
Vg=evenlySampleCurve(Vg);

%%
% Create tube
clear optStruct
optStruct.r=1.2; %Radius
optStruct.nr=25; %Number of points allong circumference
optStruct.patchType='quad'; %Output mesh type
optStruct.fixOpt=0; %Option to fix self intersection

[Fs,Vs,Cs]=polyTube(Vg,optStruct);

%%
% Visualization

cFigure;
subplot(1,2,1); hold on;
hp1=plotV(Vg,'k.-','MarkerSize',markerSize,'LineWidth',lineWidth);
axisGeom;

subplot(1,2,2); hold on;
hp2=gpatch(Fs,Vs,Cs,'k',1,1);
axisGeom; camlight headlight;

legend([hp1 hp2],{'Guide curve','Tube'});
gdrawnow;

%% Example 2: Variable radius

%%
% Creating example curve

t=linspace(0,1.5*pi,50);
r=linspace(0.25,1,50)'; %Defining radius allong curve
Vg=[t(:) zeros(size(t(:))) sin(t(:))];
Vg=evenlySampleCurve(Vg);

%%
% Create tube
clear optStruct
optStruct.r=r; %Radius
optStruct.nr=25; %Number of points allong circumference
optStruct.patchType='quad'; %Output mesh type
optStruct.fixOpt=0; %Option to fix self intersection

[Fs,Vs,Cs]=polyTube(Vg,optStruct);

%%
% Visualization

cFigure; hold on;
title('Variable radius');
gpatch(Fs,Vs,Cs,'k',1,1);
axisGeom; camlight headlight;
gdrawnow;

%% Example 3: Triangulated surface output and closed ends

%%
% Creating example curve

r=linspace(0.25,1,30)'; %Defining radius allong curve
t=linspace(0,pi,30);
Vg=[t(:) zeros(size(t(:))) -2*sin(t(:))];
Vg=evenlySampleCurve(Vg);

%%
% Create tube with open ends
clear optStruct
optStruct.r=1;
optStruct.nr=25;
optStruct.patchType='tri';
optStruct.closeOpt=0;
[Fs1,Vs1,Cs1]=polyTube(Vg,optStruct);

%%
% Create tube with closed ends
clear optStruct
optStruct.r=1;
optStruct.nr=25;
optStruct.patchType='tri';
optStruct.closeOpt=1;
[Fs2,Vs2,Cs2]=polyTube(Vg,optStruct);

%%
% Visualization

cFigure;
subplot(1,2,1); hold on;
title('Open tube')
gpatch(Fs1,Vs1,Cs1,'k',1,1);
axisGeom; camlight headlight;

subplot(1,2,2); hold on;
title('Closed tube')
gpatch(Fs2,Vs2,Cs2,'k',1,1);
axisGeom; camlight headlight;

gdrawnow;

%% Example 4: Using the fix option to attempt to fix self intersections

%%
% Creating example curve
t=linspace(0,pi,30);
Vg=[t(:) 2*sin(t(:)) zeros(size(t(:))) ];
Vg=evenlySampleCurve(Vg);

%%
% Create tube without self-intersection fix
clear optStruct
optStruct.r=1.2;
optStruct.nr=25;
optStruct.patchType='quad';
optStruct.fixOpt=0;
[Fs1,Vs1,Cs1]=polyTube(Vg,optStruct);

%%
% Create tube with self-intersection fix
clear optStruct
optStruct.r=1.2;
optStruct.nr=25;
optStruct.patchType='quad';
optStruct.fixOpt=1;
[Fs2,Vs2,Cs2]=polyTube(Vg,optStruct);

%%
% Visualization

cFigure;
subplot(1,2,1); hold on;
title('Self intersecting tube')
gpatch(Fs1,Vs1,'rw','k',0.5,2);
axisGeom; camlight headlight;
view(2);

subplot(1,2,2); hold on;
title('Fixed tube')
gpatch(Fs2,Vs2,'gw','k',1,2);
axisGeom; camlight headlight;
view(2);

gdrawnow;

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function kern = matern52KernExpandParam(kern, params)

% MATERN52KERNEXPANDPARAM Create kernel structure from MATERN52 kernel's parameters.
% FORMAT
% DESC returns a matern kernel with nu=5/2 kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : matern52KernParamInit, matern52KernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

kern.lengthScale = params(1);
kern.variance = params(2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function calcFeatureImportanceRF_Aerts(pathRF,nPerms,seed)

startpath = pwd;
cd(pathRF), load('testingVariableImportance') % 'variableImportance' gets out of there

[percentAUCdecrease,varNames] = calcFeatureImportanceRF_HN(pathRF,nPerms,'DeathSign','CT',seed);
variableImportance.DeathSign.CT.percentAUCdecrease = percentAUCdecrease;
variableImportance.DeathSign.CT.varNames = varNames;
cd(pathRF), save('testingVariableImportance','variableImportance')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%SerialLink.INERTIA Manipulator inertia matrix
%
% I = R.inertia(Q) is the symmetric joint inertia matrix (NxN) which relates
% joint torque to joint acceleration for the robot at joint configuration Q.
%
% If Q is a matrix (KxN), each row is interpretted as a joint state
% vector, and the result is a 3d-matrix (NxNxK) where each plane corresponds
% to the inertia for the corresponding row of Q.
%
% Notes::
% - The diagonal elements I(J,J) are the inertia seen by joint actuator J.
% - The off-diagonal elements I(J,K) are coupling inertias that relate
%   acceleration on joint J to force/torque on joint K.
% - The diagonal terms include the motor inertia reflected through the gear
%   ratio.
%
% See also SerialLink.RNE, SerialLink.CINERTIA, SerialLink.ITORQUE.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function M = inertia(robot, q)
    if numcols(q) ~= robot.n
        error('q must have %d columns', robot.n);
    end

    if numrows(q) > 1
        M = [];
        for i=1:numrows(q)
            M = cat(3, M, robot.inertia(q(i,:)));
        end
        return
    end

	n = robot.n;

	if numel(q) == robot.n
		q = q(:).';
	end

	M = zeros(n,n,0);
	for Q = q.'
		m = rne(robot, ones(n,1)*Q.', zeros(n,n), eye(n), 'gravity', [0 0 0]);
		M = cat(3, M, m);
	end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [min_width, min_angle] = minimumCaliperDiameter(points)
%MINIMUMCALIPERDIAMETER Minimum caliper diameter of a set of points
%
%   WIDTH = minimumCaliperDiameter(POINTS)
%   Computes the minimum width of a set of points. As polygons and
%   polylines are represented as point lists, this function works also for
%   polygons and polylines.
%
%   [WIDTH THETA] = minimumCaliperDiameter(POINTS)
%   Also returns the direction of minimum width. The direction corresponds
%   to the horizontal angle of the edge that minimizes the width. THETA is
%   given in radians, between 0 and PI.
%
%
%   Example
%   % Compute minimal caliper diameter, and check coords of rotated points
%   % have expected extent
%     points = randn(30, 2);
%     [width theta] = minimumCaliperDiameter(points);
%     points2 = transformPoint(points, createRotation(-theta));
%     diff = max(points2) - min(points2);
%     abs(width - diff(2)) < 1e-10
%     ans =
%         1
%
%   References
%   Algorithms use rotating caliper. Implementation was based on that of
%   Wikipedia:
%   http://en.wikipedia.org/wiki/Rotating_calipers
%
%   See also
%   polygons2d, convexHull, orientedBox
%

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

% first, compute convex hull of the polygon
inds = convhull(points(:,1), points(:,2));
hull = points(inds, :);

% if first and last points are the same, remove the last one
if inds(1) == inds(end)
    hull = hull(1:end-1, :);
end

% number of hull vertices
nV = size(hull, 1);

% default values
rotated_angle = 0;
min_width = inf;
min_angle = 0;

% avoid degenerated cases
if nV < 3
    return;
end

[tmp, p_a] = min(hull(:, 2)); %#ok<ASGLU>
[tmp, p_b] = max(hull(:, 2)); %#ok<ASGLU>

caliper_a = [ 1 0];    % Caliper A points along the positive x-axis
caliper_b = [-1 0];    % Caliper B points along the negative x-axis

while rotated_angle < pi
    % compute the direction vectors corresponding to each edge
    ind_a2 = mod(p_a, nV) + 1;
    vector_a = hull(ind_a2, :) - hull(p_a, :);

    ind_b2 = mod(p_b, nV) + 1;
    vector_b = hull(ind_b2, :) - hull(p_b, :);

    % Determine the angle between each caliper and the next adjacent edge
    % in the polygon
    angle_a = vectorAngle(caliper_a, vector_a);
    angle_b = vectorAngle(caliper_b, vector_b);

    % Determine the smallest of these angles
    minAngle = min(angle_a, angle_b);

    % Rotate the calipers by the smallest angle
    caliper_a = rotateVector(caliper_a, minAngle);
    caliper_b = rotateVector(caliper_b, minAngle);

    rotated_angle = rotated_angle + minAngle;

    % compute current width, and update opposite vertex
    if angle_a < angle_b
        line = createLine(hull(p_a, :), hull(ind_a2, :));
        width = distancePointLine(hull(p_b, :), line);
        p_a = mod(p_a, nV) + 1;

    else
        line = createLine(hull(p_b, :), hull(ind_b2, :));
        width = distancePointLine(hull(p_a, :), line);
        p_b = mod(p_b, nV) + 1;

    end

    % update minimum width and corresponding angle if needed
    if width < min_width
        min_width = width;
        min_angle = rotated_angle;
    end

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


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

%%
clear; close all; clc;

%% Syntax
% |[F,V]=quadPlate(plateDim,plateEl)|

%% Description
% This function generates a quadrilateral mesh for a rectangular place

%% Examples

%%
% Plot settings
fontSize=20;

%% Example: Creating a meshed plate

%Define input parameters
plateDim = [2 2]; %Plate dimensions in X and Y direction
plateEl  = [6 7]; %Number of elements in X and Y direction

%Use |quadPlate| to create a quadrilateral mesh
[F,V]=quadPlate(plateDim,plateEl);

%%
% Plotting model
cFigure; hold on;
gpatch(F,V,'bw','k');
% patchNormPlot(F,V);
axisGeom(gca,fontSize);
drawnow;

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

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


text:
% NMFLAB for Signal Processing written by A. Cichocki and R. Zdunek
% in cooperation with other members of Laboratory for Advanced Brain Signal
% Processing, BSI, RIKEN, Saitama, JAPAN

function [A,X,Distance_output]=nmf_cascade(Y,r,Index_norm,A_true,S,mc_on, NoAlts,restart_mc_on,AL,Y_true,type_alg,max_restart,no_iter, alphaA1, alphaA2, alphaX, alphaSa, alphaS, alpha)
%
%
% Non-negative Matrix Factorization (NMF) with the cascade implementation
% of some NMF algorithms
%
% [A,X]=nmf_cascade(Y, r, Index_norm, A_true, S, mc_on, NoAlts, restart_mc_on, AL, Y_true, type_alg, max_restart, no_iter,
%       alphaA1, alphaA2, alphaX, alphaSa, alphaS, alpha)
%       returns mixing matrix A of dimension [m by r],
%       and source matrix X of dimension [r by T], for the linear mixing model: AX = Y,
%       where Y is an observation matrix [m by T].
% Note: > m: number of sensors,
%       > r: number of sources,
%       > T: number of samples,
%
% INPUTS:
%       > Index_norm:    vector of 13 binary entries indicating which number
%                        divergence measures are turned on in the View Options,
%
%       > A_true:        true mixing matrix (only for synthetic data),
%       > S:             true source matrix (only for synthetic data),
%       > mc_on:         1 - Monte Carlo analysis enabled, 0 - otherwise,
%       > NoAlts:        number of alternating steps (only for Monte Carlo
%                        analysis and the option "Fixed Alternatings" is selected)
%       > restart_mc_on: 1 - restarts in Monte Carlo analysis are enabled,
%                        0 - otherwise,
%       > AL:            mixing matrix estimaed from the preceeding layer,
%       > Y_true:        the first layer mixed signals (mixtures),
%       > type_alg:      indicates the selected algorithm,
%       > max_restart:   number of restarts,
%       > no_iter:       number of inner iterations,
%       > alphaA1:       regularization parameter for computation of
%                        the first mixing matrix,
%       > alphaA2:       regularization parameter for computation of
%                        the second mixing matrix,
%       > alphaX:        regularization parameter for computation of
%                        the sources,
%       > alphaSa:       parameter of non-linear projection in computation
%                        of the mixing matrix,
%       > alphaS:        parameter of non-linear projection in computation
%                        of the sources,
%       > alpha:         parameter "alpha" in the Amari alpha-divergence,
%
% OUTPUTS:
%       > A:               estimated mixing matrix,
%       > X:               estimated source matrix,
%       > Distance_output: structures of different divergences measured between "Y" and estimated "AX" versus iterations,
%
%
% #########################################################################
A = [];
X = [];
if (nargin < 19) | isempty(alpha) | max(size(alpha) > 1)
   disp('Incorrect parameter alpha in alpha-divergence');
   return
end
if (nargin < 18) | isempty(alphaS) | max(size(alphaS) > 1)
   disp('Incorrect regularization parameter alphaS');
   return
end
if (nargin < 17) | isempty(alphaSa) | max(size(alphaSa) > 1)
   disp('Incorrect regularization parameter alphaSa');
   return
end
if (nargin < 16) | isempty(alphaX) | max(size(alphaX) > 1)
   disp('Incorrect regularization parameter alphaX');
   return
end
if (nargin < 15) | isempty(alphaA2) | max(size(alphaA2) > 1)
   disp('Incorrect regularization parameter alphaA2');
   return
end
if (nargin < 14) | isempty(alphaA1) | max(size(alphaA1) > 1)
   disp('Incorrect regularization parameter alphaA1');
   return
end
if (nargin < 13) | isempty(no_iter) | (no_iter < 1) | max(size(no_iter) > 1)
   disp('Incorrect number of iterations in the EMML algorithm');
   return
end
if (nargin < 12) | isempty(max_restart)  | (max_restart < 0) | max(size(max_restart) > 1)
   disp('Number of restarts must be given correctly');
   return
end
if (nargin < 11) | isempty(type_alg) | (type_alg < 0) | max(size(type_alg) > 1)
   disp('Unknown algorithm');
   return
end
if (nargin < 10) | isempty(Y_true)
   disp('The first layer mixed signals are unknown');
   Y_true = zeros(size(Y_true));
end
if (nargin < 9) | isempty(AL)
   disp('Mixing matrix from the preceeding layer unknown');
   AL = eye(size(Y,1));
end
if (nargin < 8) | isempty(restart_mc_on) | max(size(restart_mc_on) > 1)
   disp('Index od restarts in MC analysis unknown');
   restart_mc_on = 0;
end
if (nargin < 7) | isempty(NoAlts) | max(size(NoAlts) > 1)
   disp('Adjustable number of alternatings');
   NoAlts = [];
end
if (nargin < 6) | isempty(mc_on) | max(size(mc_on) > 1)
   disp('No Monte Carlo Analysis');
   mc_on = 0;
end
if (nargin < 5) | isempty(S)
   disp('X_true not given');
end
if (nargin < 4) | isempty(A_true)
   disp('A_true not given');
   index_fixed_A = 1;
else
   index_fixed_A = 0;
end
if (nargin < 3) | isempty(Index_norm)
   '"Index_norm" must be specified'
   return
end
if (nargin < 2) | isempty(r)
   'Rank of factorization must be given'
   return
end
if isempty(Y) | isnan(Y)
   error('No data');
   return
end
% test for negative values in Y
if min(min(Y)) < 0
    disp('Some matrix entries are changed from negative to small positive');
    Y(Y< 0) = eps;
end
if min(sum(Y,2)) == 0
    disp('Not all entries in a row can be zero');
    return
end

Y = Y + eps;
%Y = Y*Y';

if (alpha == 0) & (type_alg == 3)
    type_alg = 4;
end

[m,T]=size(Y);
niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)
niter_sample = 30; % maximum number of iterations for each random sample
epsil_normA = 1E-12; % tolerance for alternating

% Monte Carlo and alternatings adjustment
if mc_on & ~restart_mc_on
    max_restart = 0;
end
if ~isempty(NoAlts)
    niter_selected = NoAlts;
end

% Declaration for A and X
A=zeros(m,r);
Ap = zeros(m);
X=zeros(r,T);
Ainit = A;
Xinit = X;
Z = zeros(m,T);
KL_outer_temp = 0;
Z_outer = 0;
nr = 0; restart_on = 0; norm_A = 10; nr_best = -1;
m_sx = 1:m; r_sx = 1:r; T_sx = 1:T; s_dist = 0;

while (nr <= max_restart)

   % Initialize random A and X
       if ~nr & (~mc_on | restart_mc_on)
          A1_init(m_sx',m_sx) = abs(repmat(.1*sin(2*pi*.1*m_sx'),1,m) + repmat(.1*cos(2*pi*.1*m_sx),m,1) + repmat(cos(2*pi*.471*m_sx'),1,m) + repmat(sin(2*pi*.471*m_sx),m,1));
          A1_init = A1_init/max(max(A1_init));

          A2_init(m_sx',r_sx) = abs(repmat(.1*sin(2*pi*.1*m_sx'),1,r) + repmat(.1*cos(2*pi*.1*r_sx),m,1) + repmat(cos(2*pi*.471*m_sx'),1,r) + repmat(sin(2*pi*.471*r_sx),m,1));
          A2_init = A2_init/max(max(A2_init));

          Xinit(r_sx',T_sx) = abs(repmat(.1*sin(2*pi*.1*r_sx'),1,T) + repmat(.1*cos(2*pi*.1*T_sx),r,1) + repmat(cos(2*pi*.471*r_sx'),1,T) + repmat(sin(2*pi*.471*T_sx),r,1));
          Xinit = Xinit/max(max(Xinit));
       else
          A1_init=rand(m,m);
          A2_init=rand(m,r);
          Xinit=rand(r,T);
       end

        % Normalization of initial guess
       A1_init = A1_init*diag(1./sum(A1_init,1));
       A2_init = A2_init*diag(1./sum(A2_init,1));

        if (nr == max_restart)&(max_restart > 0)
           A1 = A1_best;
           A2 = A2_best;
           X = X_best;
        else
           A1 = A1_init;
           A2 = A2_init;
           X = Xinit;
        end % initial guess assignment

    Yx = zeros(m,T);
    n = 0; k = 0;

while ((k <= niter_sample)&(nr < max_restart)) | ((k <= niter_selected)&(nr == max_restart)&(norm_A > epsil_normA)& isempty(NoAlts)) | ((k <= niter_selected)&(nr == max_restart)& (NoAlts > 0))

k = k + 1;

% generalized divergence-reducing NMF iterations (main algorithm)
        if no_iter == 1

            if type_alg == 1 % KL

                Ap = A1;

                Nom = (Y./(A1*A2*X + eps))*(A2*X)';
                Nom(Nom <= 0) = eps;
                A1 = (A1.*(Nom./( repmat((A2*X*ones(T,1))',m,1) + eps))).^(1 + alphaSa);
                A1 = A1*diag(1./sum(A1,1));

                Nom = A1'*(Y./(A1*A2*X + eps))*X';
                Nom(Nom <= 0) = eps;
                A2 = (A2.*(Nom./(A1'*ones(m,1)*ones(1,T)*X' + eps))).^(1 + alphaSa);
                A2 = A2*diag(1./sum(A2,1));

                Nom = (A1*A2)'*(Y./(A1*A2*X + eps));
                Nom(Nom <= 0) = eps;
                X = (X.*(Nom./(repmat((A1*A2)'*ones(m,1),1,T) + eps))).^(1 + alphaS);

            elseif type_alg == 2 % Frobenius

                Ap = A1;
                Nom = Y*(A2*X)' - alphaA1;
                Nom(Nom <= 0) = eps;
                A1 = A1.*(Nom./(A1*A2*X*(A2*X)' + eps));
                A1 = A1*diag(1./sum(A1,1));

                Nom = A1'*Y*X' - alphaA2;
                Nom(Nom <= 0) = eps;
                A2 = A2.*(Nom./(A1'*A1*A2*X*X' + eps));
                A2 = A2*diag(1./sum(A2,1));

                Nom = (A1*A2)'*Y - alphaX;
                Nom(Nom <= 0) = eps;
                X = X.*(Nom./((A1*A2)'*A1*A2*X + eps));

            elseif type_alg == 3 % Amari alpha-divergence

                Ap = A1;
                Nom = ((Y./(A1*A2*X + eps)).^alpha )*(A2*X)';
                Nom(Nom <= 0) = eps;
                A1 = (A1.*(Nom./( repmat((A2*X*ones(T,1))',m,1) + eps)).^(1/alpha) ).^(1 + alphaSa);
                A1 = A1*diag(1./sum(A1,1));

                Nom = A1'*((Y./(A1*A2*X + eps)).^alpha )*X';
                Nom(Nom <= 0) = eps;
                A2 = (A2.*(Nom./(A1'*ones(m,1)*ones(1,T)*X' + eps)).^(1/alpha) ).^(1 + alphaSa);
                A2 = A2*diag(1./sum(A2,1));

                Nom = (A1*A2)'*(Y./(A1*A2*X + eps)).^alpha;
                Nom(Nom <= 0) = eps;
                X = (X.*(Nom./(repmat((A1*A2)'*ones(m,1),1,T) + eps)).^(1/alpha) ).^(1 + alphaS);

            elseif type_alg == 4 % SMART

                Ap = A1;
                Nom = ( log(Y./(A1*A2*X + eps) +eps) )*(A2*X)';
                A1 = (A1.*exp( Nom* diag(1./(sum( A2*X,2) +eps) ) ) ).^(1 + alphaSa);
                A1 = A1*diag(1./(sum(A1,1) + eps));

                Nom = A1'*( log(Y./(A1*A2*X + eps) + eps) )*X';
                A2 = (A2.*exp( diag(1./(sum( A1,1) + eps) )*Nom*diag(1./(sum( X,2) + eps) ) ) ).^(1 + alphaSa);
                A2 = A2*diag(1./(sum(A2,1) + eps));

                Nom = (A1*A2)'*log(Y./(A1*A2*X + eps) + eps);
                X = (X.*exp(diag(1./(sum( A1*A2,1) + eps) )*Nom ) ).^(1 + alphaS);

            end % type_alg

        else

          if type_alg == 1 % KL
                   Ap = A1;
               for t = 1:no_iter
                   Nom = (Y./(A1*A2*X + eps))*(A2*X)' - alphaA1;
                   Nom(Nom <= 0) = eps;
                   A1 = (A1.*(Nom./( repmat((A2*X*ones(T,1))',m,1) + eps))).^(1 + alphaSa);
                   A1 = A1*diag(1./sum(A1,1));

               end
               for t = 1:no_iter
                    Nom = A1'*(Y./(A1*A2*X + eps))*X' - alphaA2;
                    Nom(Nom <= 0) = eps;
                    A2 = (A2.*(Nom./(A1'*ones(m,1)*ones(1,T)*X' + eps))).^(1 + alphaSa);
                    A2 = A2*diag(1./sum(A2,1));
               end
               for t = 1:no_iter
                   Nom = (A1*A2)'*(Y./(A1*A2*X + eps)) - alphaX;
                   Nom(Nom <= 0) = eps;
                   X = (X.*(Nom./(repmat((A1*A2)'*ones(m,1),1,T) + eps))).^(1 + alphaS);
               end

          elseif type_alg == 2 % Frobenius
               Ap = A1;
               for t = 1:no_iter
                   Nom = Y*(A2*X)' - alphaA1;
                   Nom(Nom <= 0) = eps;
                   A1 = (A1.*(Nom./(A1*A2*X*(A2*X)' + eps))).^(1 + alphaSa);
                   A1 = A1*diag(1./sum(A1,1));
               end
               for t = 1:no_iter
                     Nom = A1'*Y*X' - alphaA2;
                     Nom(Nom <= 0) = eps;
                     A2 = (A2.*(Nom./(A1'*A1*A2*X*X' + eps))).^(1 + alphaSa);
                     A2 = A2*diag(1./sum(A2,1));
               end
               for t = 1:no_iter
                    Nom = (A1*A2)'*Y - alphaX;
                    Nom(Nom <= 0) = eps;
                    X = (X.*(Nom./((A1*A2)'*A1*A2*X + eps))).^(1 + alphaS);
               end

          elseif type_alg == 3 % Amari alpha-divergence
               Ap = A1;
               for t = 1:no_iter
                    Nom = ((Y./(A1*A2*X + eps)).^alpha )*(A2*X)';
                    Nom(Nom <= 0) = eps;
                    A1 = (A1.*(Nom./( repmat((A2*X*ones(T,1))',m,1) + eps)).^(1/alpha) ).^(1 + alphaSa);
                    A1 = A1*diag(1./sum(A1,1));
               end
               for t = 1:no_iter
                     Nom = A1'*((Y./(A1*A2*X + eps)).^alpha )*X';
                     Nom(Nom <= 0) = eps;
                     A2 = (A2.*(Nom./(A1'*ones(m,1)*ones(1,T)*X' + eps)).^(1/alpha) ).^(1 + alphaSa);
                     A2 = A2*diag(1./sum(A2,1));
               end
               for t = 1:no_iter
                     Nom = (A1*A2)'*(Y./(A1*A2*X + eps)).^alpha;
                     Nom(Nom <= 0) = eps;
                     X = (X.*(Nom./(repmat((A1*A2)'*ones(m,1),1,T) + eps)).^(1/alpha) ).^(1 + alphaS);
               end

            elseif type_alg == 4 % SMART
               Ap = A1;
               for t = 1:no_iter
                    Nom = ( log(Y./(A1*A2*X + eps) +eps) )*(A2*X)';
                    Nom(Nom <= 0) = eps;
                    A1 = (A1.*exp( Nom* diag(1./(sum( A2*X,2) +eps) ) ) ).^(1 + alphaSa);
                    A1 = A1*diag(1./(sum(A1,1) + eps));
               end
               for t = 1:no_iter
                    Nom = A1'*( log(Y./(A1*A2*X + eps) + eps) )*X';
                    Nom(Nom <= 0) = eps;
                    A2 = (A2.*exp( diag(1./(sum( A1,1) + eps) )*Nom*diag(1./(sum( X,2) + eps) ) ) ).^(1 + alphaSa);
                    A2 = A2*diag(1./(sum(A2,1) + eps));
               end
               for t = 1:no_iter
                    Nom = (A1*A2)'*log(Y./(A1*A2*X + eps) + eps);
                    Nom(Nom <= 0) = eps;
                    X = (X.*exp(diag(1./(sum( A1*A2,1) + eps) )*Nom ) ).^(1 + alphaS);
               end

          end % type_alg

        end % if no_iter

                if (nr == max_restart)&(mod(k,50)==0)& (~mc_on | restart_mc_on)
                    norm_A = norm(abs(A1 - Ap),'fro');
                    fprintf(1, 'Restart %d,  %d-th alternating step\n',nr_best+1,k);
                end

                if sum(Index_norm)
                   if (nr == max_restart) & (((k < 50) & (mod(k,5)==0)) | ((k>49) & ((mod(k,50)==0))))

                       s_dist = s_dist + 1;
                       k_select(s_dist) = k;
                       Z = A1*A2*X + eps;
                       Z = diag(1./(sqrt(var(Z')) + eps))*Z;

                        dist_Fro(s_dist) = norm(Y - Z,'fro');
                        dist_KL(s_dist) = sum(sum(Y.*log(Y./Z + eps) - Y + Z));
                        dist_KL2(s_dist) = sum(sum(Z.*log(Z./Y + eps) + Y - Z));
                        dist_Pearson(s_dist) = sum(sum( ((Y - Z).^2)./Z ));
                        dist_Hellinger(s_dist) = sum(sum( (sqrt(Z) - sqrt(Y)).^2 ));
                        dist_JS_rel(s_dist) = sum(sum(2*Y.*log(2*Y./(Y + Z) + eps) + Z - Y));
                        dist_JS_rel2(s_dist) = sum(sum(2*Z.*log(2*Z./(Y + Z) + eps) - Z + Y));
                        Zy = Y + Z;
                        dist_JS(s_dist) = sum(sum(Y.*log(2*Y./Zy + eps) + Z.*log(2*Z./Zy + eps) ));
                        dist_AG_rel(s_dist) = sum(sum(Zy.*log(.5*Zy./Y + eps) + Y - Z));
                        dist_AG(s_dist) = sum(sum(.5*Zy.*log(.5*Zy./sqrt(Y.*Z) + eps)));
                        dist_J(s_dist) = sum(sum( .5*(Y - Z).*log(Y./Z + eps) ));
                        dist_Chi(s_dist) = sum(sum( ((Y + Z).*(Y - Z).^2)./(Y.*Z) ));
                        dist_Tria(s_dist) = sum(sum( ((Y - Z).^2)./(Y + Z) ));
                    end % if multiple
                end % if sum

end % while (k)

% Outer KL divergence
Z = AL*A1*A2*X;
Z_outer = norm(Z,'fro') + eps;
KL_outer = sum(sum(Y_true.*log((Y_true + eps)./(Z + eps)) - Y_true + Z))/Z_outer;

          if (nr == 0) | (KL_outer < KL_outer_temp)
              A1_best = A1; A2_best = A2; X_best = X; KL_outer_temp = KL_outer; nr_best = nr;
           end % multi-conditions

   nr = nr + 1;

   if nr <=max_restart
      fprintf(1, 'Restart %d, Kullback-Leibler divergence = %e\n',	nr, KL_outer);
   end

end % while (restarts)

% One-Variance scaling
X(X <= 0) = eps;
A = A1*A2;

Distance_output = cell(length(s_dist),1);
Distance_output(1) = {[]};
Distance_output(2) = {[]};
Distance_output(3) = {[]};
Distance_output(4) = {[]};
Distance_output(5) = {[]};
Distance_output(6) = {[]};
Distance_output(7) = {[]};
Distance_output(8) = {[]};
Distance_output(9) = {[]};
Distance_output(10) = {[]};
Distance_output(11) = {[]};
Distance_output(12) = {[]};
Distance_output(13) = {[]};
Distance_output(14) = {[]};

if sum(Index_norm)
   Distance_output(1) = {k_select};
   Distance_output(2) = {dist_Fro};
   Distance_output(3) = {dist_KL};
   Distance_output(4) = {dist_KL2};
   Distance_output(5) = {dist_Pearson};
   Distance_output(6) = {dist_Hellinger};
   Distance_output(7) = {dist_JS_rel};
   Distance_output(8) = {dist_JS_rel2};
   Distance_output(9) = {dist_JS};
   Distance_output(10) = {dist_AG_rel};
   Distance_output(11) = {dist_AG};
   Distance_output(12) = {dist_J};
   Distance_output(13) = {dist_Chi};
   Distance_output(14) = {dist_Tria};
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%io_loadlcmdetail.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% [metabs,corrMatrix]=io_loadlcmdetail(filename);
%
% DESCRIPTION:
% This function loads in the "detailed output" of LCModel and returns the
% matrix of metabolite correlation coefficients.
%
% INPUTS:
% filename   = Filename of the lcmodel detailed output file.
%
% OUTPUTS:
% metabs     = A listing of the metabolites included in the correlation
%               coefficients table.
% corrMatrix = A matrix of correlation coefficients between metabolites,
%               with indices specified by the 'metabs' variable.

function [metabs,corrMatrix]=io_loadlcmdetail(filename);

%First open the detailed output file and search for the beginning of the
%Correlation Coefficients table
fid=fopen(filename);
line=fgets(fid);
search_index=findstr(line,'Correlation coefficients');

while isempty(search_index)
    line=fgets(fid);
    search_index=findstr(line,'Correlation coefficients');
end

%Now that we've found the top of the correlation coefficients table, we
%need to record the names of all of the metabolites that are listed in the
%table.  Start by going down two lines to where the metabolite names are
%listed:
line=fgets(fid);
line=fgets(fid);

%Read the first line of metabolite names and then go down one line.  The
%table of correlation coefficients can only fit 17 metabolites per row.
%Therefore if there are more than 17 metabolites to read, then we will have
%to go down another row and read again.  A simple trick to know whether to
%go down another row is to check and see if the first token in the next row
%is equal to the second token from the first row.  If so, then we're done
%reading metabolite names.  If not, then we need to keep reading.
metabstr=line;
[temp,remain]=strtok(line);
secondMetabName=strtok(remain);
line=fgets(fid);

while ~strcmp(strtok(line),secondMetabName)
    metabstr=[metabstr line];
    line=fgets(fid);
end

%We are now finished reading all of the metabolite names.  Now make a cell
%array containing the names of all of the meatbolites.  Note that at this
%point in the code, the final element in metabs will be a blank space,
%since the last metabolite name is missing from the top row of the correlation
%coefficient table in the LCModel detailed output.
remain=metabstr;
n=1;
while ~isempty(remain)
    [metabs{n},remain]=strtok(remain);
    n=n+1;
end

%Now calculate the total number of metabolites, and initialize the
%correlation coefficient matrix:
numMetabs=numel(metabs);
corrMatrix=zeros(numMetabs);
fillLine=1;

%Now start populating the correlation coefficient matrix.  The first token
%of each line is a metabolite name, so we have to remove it.  The rest of
%the line is a row in the correlation matrix.  The first 17 lines are
%simple to read becuase they all have 17 elements or less, so they don't
%wrap around the page (remember there are only 17 columns in the page).  NOTE:
%The correlation matrix is symmetric about the diagonal, but we start by only
%constructing the lower half of the diagonal (The upper half will be all
%zeros).  The diagonal elements themselves will be set to 0.5.  Then later,
%the full matrix, will be completed by adding the matrix to it's transpose,
%thus completing the symmetry, and making all of the diagonals ones.  We
%will also rename "metabs" as we go, becuase the full names of each
%metabolite are given with greater precision on the left column of the
%table in the lcmodel detailed output (compared to the top row of the table).
corrMatrix(1,1)=0.5;
if numMetabs>17
    upto=18;
else
    upto=numMetabs;
end
while fillLine<upto
    [name,corrs]=strtok(line);
    metabs{fillLine+1}=name;
    corrs=[corrs(1:end-1) ' 0.5'];
    corrMatrix(fillLine+1,1:length(str2num(corrs)))=str2num(corrs);
    line=fgets(fid);
    fillLine=fillLine+1;
end

line=fgets(fid);

%For the 18th to 34th lines, there is one line of 'wrap around'.
if numMetabs>34
    upto=35;
else
    upto=numMetabs;
end
while fillLine<upto
    [name,corrs]=strtok(line);
    metabs{fillLine+1}=name;
    line=fgets(fid);
    corrs=[corrs(1:end-1) ' ' line(1:end-1) ' 0.5'];
    corrMatrix(fillLine+1,1:length(str2num(corrs)))=str2num(corrs);
    line=fgets(fid);
    fillLine=fillLine+1;
end

%For the 35th to 52nd lines, there are two lines of 'wrap around'.  Assume
%that there will be no more than 52 elemens in the correlation matrix.
while fillLine<numMetabs
    [name,corrs]=strtok(line);
    metabs{fillLine+1}=name;
    line1=fgets(fid);
    line2=fgets(fid);
    corrs=[corrs(1:end-1) ' ' line1(1:end-1) ' ' line2(1:end-1) ' 0.5'];
    size(str2num(corrs));
    corrMatrix(fillLine+1,1:length(str2num(corrs)))=str2num(corrs);
    line=fgets(fid);
    fillLine=fillLine+1;
end
metabs{end}=name;

%Now add the correlation matrix (the lower have of which is currently only
%populated) to its transpose to complete the symmetry.
corrMatrix=corrMatrix+corrMatrix';

fclose(fid);

% imagesc(corrMatrix);
% impixelinfo;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function c = convv(a,b);
%function c = convv(a,b);
%
%Convolves 3-D matrices a and b

% S. de Waele, March 2003.

sa = kingsize(a);
sb = kingsize(b);
s = sb;
s(3) = sa(3)+sb(3)-1;
b_z = zeros(s);

b_z(:,:,1:sb(3)) = b;
c = armafilterv(b_z,1,a);

author: Sable

repo: mcbench-benchmarks


text:
% imPadded = padImage(im, padsize, padval)
%       padsize = [top, bottom, left, right]
%       padval = valid arguments for padval to padarray. e.g. 'replicate', or 0
%
%       for negative padsize, undoes the padding

%	Author:		Oliver Whyte <oliver.whyte@ens.fr>
%	Date:		November 2011
%	Copyright:	2011, Oliver Whyte
%	Reference:  O. Whyte, J. Sivic and A. Zisserman. "Deblurring Shaken and Partially Saturated Images". In Proc. CPCV Workshop at ICCV, 2011.
%	URL:		http://www.di.ens.fr/willow/research/saturation/

function imPadded = padImage(im, padsize, padval)
if nargin < 3, padval = 0; end

if any(padsize < 0)
    padsize = -padsize;
    imPadded = im(padsize(1)+1:end-padsize(2), padsize(3)+1:end-padsize(4), :);
else
    imPadded =  padarray( ...
                padarray(im, [padsize(1) padsize(3)],padval,'pre'), ...
                             [padsize(2) padsize(4)],padval,'post');
end

author: IVRL

repo: Kernel-Modeling-Super-Resolution


text:
function [MBox] = MBoxtest(X,alpha)
% Multivariate Statistical Testing for the Homogeneity of Covariance Matrices by the Box's M.
%
%   Syntax: function [MBox] = MBoxtest(X,alpha)
%
%     Inputs:
%          X - data matrix (Size of matrix must be n-by-(1+p); sample=column 1, variables=column 2:p).
%      alpha - significance level (default = 0.05).
%     Output:
%          MBox - the Box's M statistic.
%          Chi-sqr. or F - the approximation statistic test.
%          df's - degrees' of freedom of the approximation statistic test.
%          P - observed significance level.
%
%    If the groups sample-size is at least 20 (sufficiently large), Box's M test
%    takes a Chi-square approximation; otherwise it takes an F approximation.
%
%    Example: For a two groups (g = 2) with three independent variables (p = 3), we
%             are interested to test the homogeneity of covariances matrices with a
%             significance level = 0.05. The two groups have the same sample-size
%             n1 = n2 = 5.
%                                       Group
%                      ---------------------------------------
%                            1                        2
%                      ---------------------------------------
%                       x1   x2   x3             x1   x2   x3
%                      ---------------------------------------
%                       23   45   15             277  230   63
%                       40   85   18             153   80   29
%                      215  307   60             306  440  105
%                      110  110   50             252  350  175
%                       65  105   24             143  205   42
%                      ---------------------------------------
%
%             Total data matrix must be:
%          X=[1 23 45 15;1 40 85 18;1 215 307 60;1 110 110 50;1 65 105 24;
%          2 277 230 63;2 153 80 29;2 306 440 105;2 252 350 175;2 143 205 42];
%
%             Calling on Matlab the function:
%                MBoxtest(X,0.05)
%
%             Answer is:
%
%  ------------------------------------------------------------
%       MBox         F           df1          df2          P
%  ------------------------------------------------------------
%     27.1622     2.6293          6           463       0.0162
%  ------------------------------------------------------------
%  Covariance matrices are significantly different.
%

%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%             And the special collaboration of the post-graduate students of the 2002:2
%             Multivariate Statistics Course: Karel Castro-Morales, Alejandro Espinoza-Tenorio,
%             Andrea Guia-Ramirez, Raquel Muniz-Salazar, Jose Luis Sanchez-Osorio and
%             Roberto Carmona-Pina.
%  November 2002.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A., R. Hernandez-Walls, K. Castro-Morales, A. Espinoza-Tenorio, A. Guia-Ramirez
%    and R. Carmona-Pina. (2002). MBoxtest: Multivariate Statistical Testing for the Homogeneity of
%    Covariance Matrices by the Box's M. A MATLAB file. [WWW document]. URL http://www.mathworks.com/
%    matlabcentral/fileexchange/loadFile.do?objectId=2733&objectType=FILE
%
%  References:
%
%  Stevens, J. (1992), Applied Multivariate Statistics for Social Sciences. 2nd. ed.
%              New-Jersey:Lawrance Erlbaum Associates Publishers. pp. 260-269.

if nargin < 1,
    error('Requires at least one input arguments.');
end;

if nargin < 2,
    alpha = 0.05; %(default)
end;

if (alpha <= 0 | alpha >= 1)
   fprintf('Warning: significance level must be between 0 and 1\n');
   return;
end;

g = max(X(:,1)); %Number of groups.

n = []; %Vector of groups-size.
indice = X(:,1);
for i = 1:g
   Xe = find(indice==i);
   eval(['X' num2str(i) '= X(Xe,2);']);
   eval(['n' num2str(i) '= length(X' num2str(i) ') ;'])
   eval(['xn= n' num2str(i) ';'])
   n = [n,xn];
end;

[f,c] = size(X);
X = X(:,2:c);

[N,p]=size(X);
r=1;
r1=n(1);
bandera=2;
for k=1:g
   if n(k)>=20;
      bandera=1;
   end
end
%Partition of the group covariance matrices.
for k=1:g
   eval(['S' num2str(k) '=cov(X(r:r1,:));';]);
   if k<g
    r=r+n(k);
    r1=r1+n(k+1);
   end
end

deno=sum(n)-g;
suma=zeros(size(S1));

for k=1:g
   eval(['suma =suma + (n(k)-1)*S' num2str(k) ';']);
end

Sp=suma/deno;  %Pooled covariance matrix.
Falta=0;

for k=1:g
   eval(['Falta =Falta + ((n(k)-1)*log(det(S' num2str(k) ')));']);
end

MB=(sum(n)-g)*log(det(Sp))-Falta;  %Box's M statistic.
suma1=sum(1./(n(1:g)-1));
suma2=sum(1./((n(1:g)-1).^2));
C=(((2*p^2)+(3*p)-1)/(6*(p+1)*(g-1)))*(suma1-(1/deno));  %Computing of correction factor.
if bandera==1
X2=MB*(1-C);  %Chi-square approximation.
v=(p*(p+1)*(g-1))/2;  %Degrees of freedom.
P=1-chi2cdf(X2,v);  %Significance value associated to the observed Chi-square statistic.
disp(' ')
;
fprintf('------------------------------------------------\n');
disp('     MBox     Chi-sqr.         df          P')
fprintf('------------------------------------------------\n');
fprintf('%10.4f%11.4f%12.i%13.4f\n',MB,X2,v,P);
fprintf('------------------------------------------------\n');
   if P >= alpha
      disp('Covariance matrices are not significantly different.');
   else
      disp('Covariance matrices are significantly different.');
   end
else
%To obtain the F approximation we first define Co, which combined to the before C value
%are used to estimate the denominator degrees of freedom (v2); resulting two possible cases.
Co=(((p-1)*(p+2))/(6*(g-1)))*(suma2-(1/(deno^2)));
  if Co-(C^2)>= 0;
	v1=(p*(p+1)*(g-1))/2;  %Numerator degrees of freedom.
   v21=fix((v1+2)/(Co-(C^2)));  %Denominator degrees of freedom.
   F1=MB*((1-C-(v1/v21))/v1);  %F approximation.
   P1=1-fcdf(F1,v1,v21);  %Significance value associated to the observed F statistic.
disp(' ')
  ;
fprintf('------------------------------------------------------------\n');
disp('     MBox         F           df1          df2          P')
fprintf('------------------------------------------------------------\n');
fprintf('%10.4f%11.4f%11.i%14.i%13.4f\n',MB,F1,v1,v21,P1);
fprintf('------------------------------------------------------------\n');
   if P1 >= alpha
      disp('Covariance matrices are not significantly different.');
   else
      disp('Covariance matrices are significantly different.');
   end

 else
   v1=(p*(p+1)*(g-1))/2;  %Numerator degrees of freedom.
   v22=fix((v1+2)/((C^2)-Co));  %Denominator degrees of freedom.
   b=v22/(1-C-(2/v22));
   F2=(v22*MB)/(v1*(b-MB));  %F approximation.
   P2=1-fcdf(F2,v1,v22);  %Significance value associated to the observed F statistic.
disp(' ')
  ;
fprintf('------------------------------------------------------------\n');
disp('     MBox         F           df1          df2          P')
fprintf('------------------------------------------------------------\n');
fprintf('%10.4f%11.4f%11.i%14.i%13.4f\n',MB,F2,v1,v22,P2);
fprintf('------------------------------------------------------------\n');

   if P2 >= alpha
      disp('Covariance matrices are not significantly different.');
   else
      disp('Covariance matrices are significantly different.');
   end
 end
end

repo: REST


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

prob = sigmoid(X * theta);
pos = find(prob >= 0.5);
p(pos,1) = 1;

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

end

author: scruel

repo: Notes-ML-AndrewNg


text:
function [y] = local2globalPosM(x, X)

% SYNTAX:
%   [y] = local2globalPos(x, X);
%
% INPUT:
%   x = local position vector(s)
%   X = origin vector(s)
%
% OUTPUT:
%   y = global position vector(s)
%
% DESCRIPTION:
%   Rototraslation from local-level reference frame to Earth-fixed reference frame

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%initialize new position vector
y = zeros(size(x));

    %geodetic coordinates
    [phi, lam] = cart2geod(X(1,i), X(2,i), X(3,i));

    %rotation matrix from global to local reference system
    R = [-sin(lam) cos(lam) 0;
         -sin(phi)*cos(lam) -sin(phi)*sin(lam) cos(phi);
         +cos(phi)*cos(lam) +cos(phi)*sin(lam) sin(phi)];

    %rototraslation
    %y(:,i) = R\x(:,i) + X(:,i);
    y = R * x' + X;

end

author: goGPS-Project


text:
function [pp, err] = MitsunagaNayarCRFFull(stack_samples, stack_exposure, N, maxIterations)
%
%       [pp, err] = MitsunagaNayarCRFFull(stack, stack_exposure, N, nSamples, sampling_strategy)
%
%       This function computes camera response function using Mitsunaga and
%       Nayar method.
%
%        Input:
%           -stack_samples: a stack of samples from LDR images
%           -nSamples: number of samples for computing the CRF
%           -N: polynomial degree of the inverse CRF
%           -maxIterations:
%
%        Output:
%           -pp: a polynomial encoding the inverse CRF
%           -err: the error
%
%     Copyright (C) 2016  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

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

threshold = 1e-4;

col = size(stack_samples, 3);

Q = length(stack_exposure);

Mmax = 1.0;

%recovering the CRF
function d = MN_d(c, q1, q2, n)

    M_q   = stack_samples(:, q1, c);
    M_q_p = stack_samples(:, q2, c);

    indx = find(M_q > 0.0 & M_q_p > 0.0);

    d = M_q(indx).^n - R(q1, q2) * (M_q_p(indx).^n);
end

pp = zeros(N + 1, col);
pp_prev = zeros(N + 1, col);

x = (0:255) / 255;

err = 0.0;

R0 = ones(Q - 1, Q - 1);
for q1=1:(Q - 1)
    for q2=1:(Q - 1)
        if(q1 ~= q2)
            R0(q1, q2) = stack_exposure(q1) / stack_exposure(q2);
        end
    end
end

for c=1:col

    R = R0;

    bLoop = 1;
    iter = 0;

    while(bLoop)
        A = zeros(N, N);
        b = zeros(N, 1);

        for i=1:N
            %init A
            for j=1:N
                for q1=1:(Q - 1)
                    for q2=1:(Q - 1)
                        if(q1 ~= q2)
                            delta  = MN_d(c, q1, q2, j - 1) - MN_d(c, q1, q2, N);
                            A(i,j) = A(i,j) + sum(MN_d(c, q1, q2, i - 1) .* delta);
                        end
                    end
                end
            end

            %init b
            for q1=1:(Q - 1)
                for q2=1:(Q - 1)
                    if(q1 ~= q2)
                        b(i) = b(i) - sum(Mmax * MN_d(c, q1, q2, i - 1) .* MN_d(c, q1, q2, N));
                    end
                end
            end
        end

        coeff = A \ b;
        coeff_n = Mmax - sum(coeff);

        pp(:,c) = flip([coeff; coeff_n]);

        f_1 = polyval(pp(:,c),      x);
        f_2 = polyval(pp_prev(:,c), x);
        bLoop = max(abs(f_1 - f_2) > threshold);

        bLoop = bLoop & (iter < maxIterations);

        if(bLoop)
            pp_prev = pp;

            %update R
            for q=1:(Q - 1)
                for q2=1:(Q - 1)
                    if(q1 ~= q2)
                        s1 = stack_samples(:, q1    , c);
                        s2 = stack_samples(:, q2, c);

                        indx = find(s1 > 0.0 & s2 > 0.0);

                        e1 = polyval(pp(:,c), s1(indx));
                        e2 = polyval(pp(:,c), s2(indx));
                        R(q1, q2) = sum(e1) / sum(e2);
                    end
                end
            end

            iter = iter + 1;
        end
    end

    %compute err
    for q1=1:(Q - 1)
        for q2=1:(Q - 1)
            if(q1 ~= q2)
                s1 = stack_samples(:, q1, c);
                s2 = stack_samples(:, q2, c);

                indx = find(s1 > 0.0 & s2 > 0.0);

                e1 = polyval(pp(:, c), s1(indx));
                e2 = polyval(pp(:, c), s2(indx));
                err = err + sum((e1 - R(q1, q2) * e2).^2);
            end
        end
    end
end

end

author: banterle


text:
classdef TestWarpAffine
    %TestWarpAffine

    properties (Constant)
        img = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_1
            im = imread(TestWarpAffine.img);
            M = eye(2,3);  % identity transform
            dst = cv.warpAffine(im, M);
            validateattributes(dst, {class(im)}, {'size',size(im)});
            assert(isequal(im, dst));
        end

        function test_2
            im = cv.imread(TestWarpAffine.img, 'Grayscale',true);
            M = eye(2,3);  % identity transform
            dst = cv.warpAffine(im, M, 'DSize',[256 256], ...
                'Interpolation','Linear', 'WarpInverse',false, ...
                'BorderType','Constant', 'BorderValue',0);
            validateattributes(dst, {class(im)}, {'size',[256 256]});
        end

        function test_3
            im = imread(TestWarpAffine.img);
            [h,w,~] = size(im);
            M = cv.getRotationMatrix2D([w h]./2, 30, 0.7);
            dst = cv.warpAffine(im, M, 'Dst',im, 'BorderType','Transparent');
            validateattributes(dst, {class(im)}, {'size',size(im)});
        end

        function test_error_argnum
            try
                cv.warpAffine();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [in3] = cc2in3(cc)
% Convert volume from cubic centimeters to cubic inches.
% Chad Greene 2012
in3 = cc*0.061023744095;


text:
function gpuNet = update_net(gpuNet, resBack, lr, numDatabase, batchSize)

weightDecay = 5*1e-4 ;
numLayers = numel(gpuNet.layers) ;
for ii = 1:numLayers
    if isfield(gpuNet.layers{ii},'weights')
        gpuNet.layers{ii}.weights{1} = gpuNet.layers{ii}.weights{1} - ...
            lr*(resBack(ii).dzdw{1}/(batchSize * numDatabase) + weightDecay*gpuNet.layers{ii}.weights{1});
        gpuNet.layers{ii}.weights{2} = gpuNet.layers{ii}.weights{2} - ...
            lr*(resBack(ii).dzdw{2}/(batchSize * numDatabase) + weightDecay*gpuNet.layers{ii}.weights{2});
    end
end
end

author: jiangqy


text:
%%% INPUT THE IMAGE FILE NAME:

graphout = 0;

if ~exist('fc')|~exist('cc')|~exist('kc')|~exist('alpha_c'),
   fprintf(1,'No intrinsic camera parameters available. Maybe, need to load Calib_Results.mat\n');
   return;
end;

KK = [fc(1) alpha_c*fc(1) cc(1);0 fc(2) cc(2) ; 0 0 1];

disp('Program that undistorts a whole sequence of images (works with bmp only so far... needs some debugging)');
disp('The intrinsic camera parameters are assumed to be known (previously computed)');
disp('After undistortion, the intrinsic parameters fc, cc, alpha_c remain unchanged. The distortion coefficient vector kc is zero');

dir;

fprintf(1,'\n');

seq_name = input('Basename of sequence images (without number nor suffix): ','s');

format_image_seq = '0';

while format_image_seq == '0',
    format_image_seq =  input('Image format: ([]=''r''=''ras'', ''b''=''bmp'', ''t''=''tif'', ''p''=''pgm'', ''j''=''jpg'', ''m''=''ppm'') ','s');
    if isempty(format_image_seq),
        format_image_seq = 'ras';
    else
        if lower(format_image_seq(1)) == 'm',
            format_image_seq = 'ppm';
        else
            if lower(format_image_seq(1)) == 'b',
                format_image_seq = 'bmp';
            else
                if lower(format_image_seq(1)) == 't',
                    format_image_seq = 'tif';
                else
                    if lower(format_image_seq(1)) == 'p',
                        format_image_seq = 'pgm';
                    else
                        if lower(format_image_seq(1)) == 'j',
                            format_image_seq = 'jpg';
                        else
                            if lower(format_image_seq(1)) == 'r',
                                format_image_seq = 'ras';
                            else
                                disp('Invalid image format');
                                format_image_seq = '0'; % Ask for format once again
                            end;
                        end;
                    end;
                end;
            end;
        end;
    end;
end;

ima_sequence = dir( [ seq_name '*.' format_image_seq]);

if isempty(ima_sequence),
    fprintf(1,'No image found\n');
    return;
end;

ima_name = ima_sequence(1).name;
if format_image_seq(1) == 'p',
    if format_image_seq(2) == 'p',
        I = double(loadppm(ima_name));
    else
        I = double(loadpgm(ima_name));
    end;
else
    if format_image_seq(1) == 'r',
        I = readras(ima_name);
    else
        I = double(imread(ima_name));
    end;
end;

[ny,nx,nc] = size(I);

% Pre-compute the necessary indices and blending coefficients to enable quick rectification:
[Irec_junk,ind_new,ind_1,ind_2,ind_3,ind_4,a1,a2,a3,a4] = rect_index(zeros(ny,nx),eye(3),fc,cc,kc,alpha_c,KK);

n_seq = length(ima_sequence);

for kk = 1:n_seq,

    ima_name = ima_sequence(kk).name;

            fprintf(1,'Loading original image %s...',ima_name);

    %%% READ IN IMAGE:

    if format_image_seq(1) == 'p',
        if format_image_seq(2) == 'p',
            I = double(loadppm(ima_name));
        else
            I = double(loadpgm(ima_name));
        end;
    else
        if format_image_seq(1) == 'r',
            I = readras(ima_name);
        else
            I = double(imread(ima_name));
        end;
    end;

    [ny,nx,nc] = size(I);

    if graphout,
        figure(2);
        image(uint8(I));
        drawnow;
    end;

    I2 = zeros(ny,nx,nc);

    for ii = 1:nc,

        Iii = I(:,:,ii);
        I2ii = zeros(ny,nx);

        I2ii(ind_new) = uint8(a1 .* Iii(ind_1) + a2 .* Iii(ind_2) + a3 .* Iii(ind_3) + a4 .* Iii(ind_4));

        I2(:,:,ii) = I2ii;

    end;

    I2 = uint8(I2);

    if graphout,
        figure(3);
        image(I2);
        drawnow;
    end;

    ima_name2 = ['undist_' ima_name];

    fprintf(1,'Saving undistorted image under %s...\n',ima_name2);

    if format_image_seq(1) == 'p',
        if format_image_seq(2) == 'p',
            saveppm(ima_name2,I2);
        else
            savepgm(ima_name2,I2);
        end;
    else
        if format_image_seq(1) == 'r',
            writeras(ima_name2,I2,gray(256));
        else
            imwrite(I2,ima_name2,format_image_seq);
        end;
    end;

end;

repo: ekfmonoslam


text:
%%
% function eval_y_and_dy: (Np,1) --> (Ny,1)x(Ny,Np), [y,J]=F(p), value y and jacobian J=dF(p)/dp of the function.
%
% p ... size ((11+3)*K+3*N,1), p = [iP(1); u0(1); kappa(1); .. iP(K); u0(K); kappa(K); iX(1); .. iX(N)],
%   where Np = (11+2+raddeg)*K + 3*N
%
% Optional parameters:
%   TP{1:K} ... size (12,11), Pk(:) = P0k + TP{k}*iPk
%   TX{1:N} ... size (4,3), Xn = X0n + TX{n}*iXn
%   qivis ... size (K,N), has 1 in entries corresponding to visible image points
% y ... size (2*nnz(qivis),1), visible image observations, stacked as [q(1,1); q(2,1); .. q(2*K,N)], unobserved q are omitted
%%
% We consider the imaging model as
% u = hom(x)
% x = P*X
% P = P0 + TP*iP
% X = X0 + TX*iX
%%
function [y,J] = eval_y_and_dy(p,P0,TP,X0,TX,y,qivis,RADIAL)
%
[K,N] = size(qivis);
NR = RADIAL*3; % number of radial distortion paramters
%
% rearrange p to P, X [,radial]
for k = 1:K
  P(k2i(k),:) = P0(k2i(k),:) + reshape(TP{k}*p([1:11]+(k-1)*(11+NR)),[3 4]);
  if RADIAL
    radial(k).u0 = p([12:13]+(k-1)*(11+NR));
    radial(k).kappa = p([14]+(k-1)*(11+NR));
  end
end
X = zeros(4,N);
for n = 1:N
  X(:,n) = TX{n}*p((11+NR)*K+[1:3]+(n-1)*3);
end
X = X0 + X;
%
% compute retina points q
for k = 1:K
  x(k2i(k),:) = P(k2i(k),:)*X;
  q(k2i(k,2),:) = nhom(x(k2i(k),:));
  if RADIAL
	q(k2i(k,2),:) = raddist(q(k2i(k,2),:),radial(k).u0,radial(k).kappa);
  end
end
qq = reshape(q,[2 K*N]);
qq = qq(:,qivis);
y = qq(:)-y; % function value
%
% compute Jacobian J = dF(p)/dp
if nargout<2
  return
end
[kvis,nvis] = find(qivis);
Ji = zeros(2*length(kvis)*(11+NR+3),1); Jj = zeros(size(Ji)); Jv = zeros(size(Ji));
cnt = 0;
for l = 1:length(kvis)  % loop for all VISIBLE points in all cameras
  k = kvis(l);
  n = nvis(l);
  xl = x(k2i(k),n);
  ul = nhom(xl);

  % Compute derivatives (Jacobians). Notation: E.g., dudx = du(x)/dx, etc.
  dxdP = kron(X(:,n)',eye(3))*TP{k}; % dx(iP,iX)/diP
  dxdX = P(k2i(k),:)*TX{n}; % dx(iP,iX)/diX
  dudx = [eye(2) -ul]/xl(3);
  if RADIAL
    [dqdu,dqdu0,dqdkappa] = raddist(ul,radial(k).u0,radial(k).kappa);
  else
    dqdu = eye(2); dqdu0 = []; dqdkappa = [];
  end
  dqdP = dqdu*dudx*dxdP;  % dq(iP,iX)/diP
  dqdX = dqdu*dudx*dxdX;         % dq(iP,iX)/dX
  c = cnt+[1:2*(11+NR+3)];
  [Ji(c),Jj(c),Jv(c)] = spidx([1:2]+(l-1)*2,[[1:(11+NR)]+(k-1)*(11+NR) (11+NR)*K+[1:3]+(n-1)*3],[dqdP dqdu0 dqdkappa dqdX]);
  cnt = cnt + 2*(11+NR+3);
end
J = sparse(Ji,Jj,Jv,length(y),length(p));
return

function [i,j,v] = spidx(I,J,V)
%
i = I'*ones(1,length(J));
j = ones(length(I),1)*J;
i = i(:);
j = j(:);
v = V(:);
return

author: strawlab

repo: MultiCamSelfCal


text:
function [l, L_rf, L_sf, L_sk, L_hm, L_beta] = ...
    retroProjPlkLinFromPinHoleOnRob(Rf, Sf, Sk, hm, beta)

% RETROPROJPLKLINFROMPINHOLEONROB Retro-project Plucker line from pinhole on robot.
%
%   L = RETROPROJPLKLINFROMPINHOLEONROB(RF, SF, SK, HML, beta) gives the
%   retroprojected Plucker line L in World Frame from an observed
%   homogeneous line HM. RF and SF are Robot and Sensor Frames, SK is the
%   camera calibration parameters vector, HML is the detected homogeneous 2D
%   line and BETA is the non-measurable prior. L is a 6-vector :
%     L = [nx ny nz vx vy vz]'
%   with the 6 Plucker coordinates of a line in 3D space.
%
%   [L, L_rf, L_sf, L_k, L_hm, L_n] = ... returns the
%   Jacobians wrt RF.x, SF.x, SK, SC, HML and N.
%
%   See also INVPINHOLEIDP, FROMFRAMEIDP.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

% Frame World -> Robot  :  Rf
% Frame Robot -> Sensor :  Sf

if nargout == 1

    % L in Sensor Frame
    ls = invPinHolePlucker(Sk,hm,beta) ;

    % L in world frame
    lr = fromFramePlucker(Sf, ls) ;
    l  = fromFramePlucker(Rf, lr) ;

else

    % L in Sensor Frame
    [ls, LS_sk, LS_hm, LS_beta] = invPinHolePlucker(Sk,hm,beta) ;

    % L in world frame
    [lr, LR_sf, LR_ls] = fromFramePlucker(Sf, ls) ;
    [l,  L_rf,  L_lr]  = fromFramePlucker(Rf, lr) ;

    % chain rule for Jacobians
    L_sf   = L_lr*LR_sf ;
    L_ls   = L_lr*LR_ls ;
    L_sk   = L_ls*LS_sk ;
    L_hm   = L_ls*LS_hm ;
    L_beta = L_ls*LS_beta ;

end

return

%% test
syms rx ry rz ra rb rc rd sx sy sz sa sb sc sd u0 v0 au av h1 h2 h3 b1 b2 real

Rf.x = [rx;ry;rz;ra;rb;rc;rd];
Sf.x = [sx;sy;sz;sa;sb;sc;sd] ;

Rf   = updateFrame(Rf);
Sf   = updateFrame(Sf);
Sk   = [u0;v0;au;av];
hm   = [h1 h2 h3]';
beta = [b1 b2]';

%% test jacobian
[l,L_r,L_s,L_k,L_hm,L_beta] = retroProjPlkLinFromPinHoleOnRob(Rf,Sf,Sk,hm,beta)

%% down here tha Jac test - WARNING! IT TAKES AGES TO COMPUTE !!
simplify(L_r - jacobian(l,Rf.x))
simplify(L_s - jacobian(l,Sf.x))
% simplify(L_k - jacobian(l,Sk))
% simplify(L_hm - jacobian(l,hm))

%% numerical test
Rf.x = epose2qpose([0 0 0 deg2rad([0 0 0])]');
Sf.x = epose2qpose([0 0 0 deg2rad([90 0 90])]');
Rf   = updateFrame(Rf);
Sf   = updateFrame(Sf);
Sk   = [100 100 100 100]';
hm   = [1 0 0]';
beta = [1 0]';

[l,L_r,L_s,L_k,L_hm,L_beta]  = retroProjPlkLinFromPinHoleOnRob(Rf,Sf,Sk,hm,beta)
[l2,L2_s,L2_k,L2_hm,L2_beta] = retroProjectPlucker(Sf,Sk,hm,beta)

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function i4vec_index_search_test ( )

%*****************************************************************************80
%
%% I4VEC_INDEX_SEARCH_TEST tests I4VEC_INDEX_SEARCH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n_max = 20;
  n = 0;
  x = [];
  indx = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_INDEX_SEARCH_TEST\n' );
  fprintf ( 1, '  I4VEC_INDEX_SEARCH searches for an entry with\n' );
  fprintf ( 1, '  a given value.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generate some random values:\n' );

  b = 0;
  c = n_max;
  seed = 123456789;

  for i = 1 : n_max
    [ xval, seed ] = i4_uniform_ab ( b, c, seed );
    [ n, x, indx ] = i4vec_index_insert_unique ( n, x, indx, xval );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', i, indx(i), x(i), x(indx(i)) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Results of search for given XVAL:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  XVAL  Less Equal More\n' );
  fprintf ( 1, '\n' );

  for xval = 0 : 20
    [ less, equal, more ] = i4vec_index_search ( n, x, indx, xval );
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', xval, less, equal, more );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function par_est2a
% parameter estimation with derivatives                      Holzbecher September 2005
% using differential equations
% Idea from FEMLAB - there R instead of Q
% see COMSOL News, Nr. 1, 2005, page 15
global xfit cfit Q

% specify fitting data
xfit = [0.05:0.1:0.95];
cfit = [0.9256859756097451       0.7884908536585051       0.6665396341462926...
        0.559832317073104        0.4683689024389414       0.39214939024380824...
        0.33117378048770196      0.28544207317062964      0.25495426829258294      0.23971036585356142];
Q = -2;
D = fzero(@myfun,1.8);
display (['Best fit for D = ' num2str(D)]);
x = [0:0.01:1];
plot (xfit,cfit,'o',x,-(Q/D/2)*x.*x + (Q/D)*x + 1,'-');
legend ('given','modelled');
xlabel ('x'); ylabel ('c');

function f = myfun(D);
global xfit cfit Q

options = bvpset;
% solve diffusion equation for c with c(0)=1 and dc/dx(1)=0
solinit = bvpinit([0 xfit 1],@guess);
c = bvp4c (@mat4ode,@mat4bc,solinit,options,Q/D,1);
%plot (c.x,c.y(1,:),'r',xfit,-(Q/D/D)*xfit.*xfit+(Q/D)*xfit+ones(1,size(xfit,2)));

% solve Poisson equation for dc/dD (cD) with boundary conditions
solinit = bvpinit([0 xfit 1],@guess1);
cD = bvp4c (@mat4ode,@mat4bc,solinit,options,Q/D/D,0);

% specify function f to vanish
f = 2*(c.y(1,2:size(c.y,2)-1)-cfit)*cD.y(1,2:size(c.y,2)-1)';

function dydx = mat4ode(x,y,Q,c0)
dydx = [y(2); -Q];
% ------------------------------------------------------------
function res = mat4bc(y0,y1,Q,c0)
res = [y0(1)-c0; y1(2)];
% ------------------------------------------------------------
function v = guess(x)
v = [x*(x-2)+1; 2*(x-1)];
% ------------------------------------------------------------
function v = guess1(x)
v = [x*(x-2); 2*(x-1)];


text:
function [ R t s ]=computeOrientation(x,xp,method,varargin)
% Compute Absolute or Exterior Orientation (Pose Estimation)
%
% Absolute is : given x and xp (2 arrays of 3D coordinates), find the best
% transformation such that x=s*R*xp+t
%
% Exterior is : 1 array of 3D position x is known and its 2D orthographic
% projection xp.
% Find the best transformation such that xp=projection*(s*R*x+t)
% (same as Pose Estimation, ePNP). Projection is ortho for now
%
% The routines below are only for the orthographic case for now
%
% USAGE
%  [R,t,s]=computeOrientation(x,xp,method)
%
% INPUTS
%  x,xp       - 3xN or 2xN array of points or right/left.
%               For 'exteriorSequence', it can be of size
%               2 x nPoint x nFrame and 3 x nPoint x nFrame
%  method     - 'absolute', 'absoluteHard' (try more if Horn's fails)
%               or 'exterior' and 'exteriorSequence'
%  varargin   - list of paramaters in quotes alternating with their values
%       - 'RIni' initial rotation from which the optimization will be
%       started (only for 'exterior' and 'exteriorSequence')
%
% OUTPUTS
%  R       - rotation matrix
%  t       - translation vector
%  s       - scale factor (if not requested, s=1)
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

nPoint=size(x,2);

switch method
  case { 'absolute', 'absoluteHard' }
    % Computes the absolute orientation between 2 sets of 3D points
    % Xright=s*R*Xleft + t   . We want to recover s,R,T
    % Ref: B.K.P. Horn, H.M. Hilden, and S. Negahdaripour, Closed-Form
    % Solution of Absolute Orientation Using Orthonormal Matrices
    rr=normalizePoint(x,4); rl=normalizePoint(xp,4);
    rrBar=mean(rr,2); rlBar=mean(rl,2);
    rrp=rr-rrBar(:,ones(1,nPoint)); rlp=rl-rlBar(:,ones(1,nPoint));

    M=zeros(3); for i=1:nPoint; M=M+rrp(:,i)*rlp(:,i)'; end
    [V,D]=eig(M'*M);
    V=real(V);

    R=1/sqrt(D(2,2))*V(:,2)*V(:,2)' + 1/sqrt(D(3,3))*V(:,3)*V(:,3)';
    temp=V(:,1)*V(:,1)';
    if D(1,1)>0
      R=R+1/sqrt(D(1,1))*temp;
    else
      if det(R+temp)>0; R=R+temp; else R=R-temp; end
    end
    R=M*R;

    R=rotationMatrix(R);

    if det(R)<0 && strcmp(method,'absoluteHard') % if Horn's method fails
      warning('Horn''s method failed. Trying GloptiPoly3 solution');
      RHorn=R;
      try
        % Sedumi and GloptiPoly3 must be installed and in the path !
        mpol R 3 3;

        g0=rlp-R*rrp;
        g0=g0(1,:)*g0(1,:)'+g0(2,:)*g0(2,:)'+g0(3,:)*g0(3,:)';

        % define the rotation constraints
        K = [ R(1,:)*R(1,:)' == 1, R(2,:)*R(2,:)'==1, R(3,:)*R(3,:)'==1,...
          R(1,:)*R(2,:)' == 0, R(1,:)*R(3,:)' == 0, R(2,:)*R(3,:)' == 0,...
          det(R)>=0 ];

        % define the problem and solve it
        P=msdp(min(g0),K);
        [ status obj ] = msol(P);

        R=rotationMatrix( double(R) );
      catch %#ok<CTCH>
        R = RHorn;
      end
    end

    % Figure out s and t
    if nargout==2
      s=1;
    else
      s=norm(rrp,'fro')/norm(rlp,'fro');
      if norm( rrp+s*R*rlp, 'fro' ) < norm( rrp-s*R*rlp, 'fro' ); s=-s; end
    end

    t=rrBar-s*R*rlBar;
  case 'exterior'
    if ~isempty(varargin); RIni = varargin{0};
    else RIni=[];
    end

    if isempty(RIni)
      try
        % Sedumi and GloptiPoly3 must be installed and in the path !
        mpol R 2 3;
        if nargout<2; t=0; else mpol t 2 1; end
        if nargout<3; s=1; else mpol s; end

        g0=0;
        for i=1:nPoint
          tmp=xp(:,i) - s*(R*x(:,i)+t); %#ok<NODEF>
          g0 = g0 + tmp'*tmp;
        end

        % define the rotation constraints
        K = [ R(1,:)*R(1,:)' == 1, R(2,:)*R(2,:)'==1, R(1,:)*R(2,:)' == 0];

        % define the problem and solve it
        mset('verbose',false);
        P=msdp(min(g0),K);

        [ status obj ] = msol(P);

        R=rotationMatrix( double(R) );

        if nargout>=2; t=[ double(t); 0]; end
        if nargout==3; s=double(s); end
      catch %#ok<CTCH>
        warning(['GloptiPoly3 not installed or Gloptypoly crashed,' ...
          'using ePnP']);
        try
          mset clear;
        catch %#ok<CTCH>
          [ R t ] = efficient_pnp( x', xp', eye(3,3) );
          t=t(1:2);
        end
      end
      if any(isnan(R))
        warning(['GloptiPoly failed, using ePnP']);
        [ R t ] = efficient_pnp( x', xp', eye(3,3) );
        t=t(1:2);
      end
    else
      R=RIni;
    end

    % perform gradient descent to optimize the rotation
    Q = refineExteriorOrientation(x,xp,quaternion(R));
    R=quaternion( Q );
  case 'exteriorSequence'
    RIni=getPrmDflt( varargin, {'RIni' [] }, 1 );
    if length(RIni)>1
      Q=quaternion(RIni);
    else
      Q=2*rand(4,size(xp,3))-1;
    end
    [ Q err ]=refineExteriorOrientation(x,xp,Q);

    for i = 1 : 10
      [ QNew errNew ]=refineExteriorOrientation(x,xp,2*rand(4,size(xp,3))-1);
      temp = errNew<err;
      Q(:,temp) = QNew(:,temp);
    end

    R=quaternion( Q );
end
end

author: vrabaud


text:
function sonorobot_keyboard(r,ns,ds,mapoff)
if nargin<4
    mapoff = false;
    if nargin<3
        ds = .1;
        if nargin<2
            ns = 24;
            if nargin<1;
                r = .37;
            end
        end
    end
end
f_h = figure('name','Sono Robot','numbertitle','off','menubar','none');
set(f_h,'keypressfcn',@robot_move_fcn)
uicontrol(f_h,'style','toggle','position',[505,5,50,25],...
    'string','Quit','callback','close');
msg1 = uicontrol(f_h,'style','text','position',[15,5,480,55]);
msg2 = uicontrol(f_h,'style','text','position',[15,65,480,25],...
    'string','Measurement results...','fontsize',11.5,...
    'horizontal','left');
hold on
axisx = 12.5;
axisy = 10;
axis([0 axisx 0 axisy]);
axis off
%Input walls & robot position
%Environment
nmax = 50;
stpmax = 100;
segset = zeros(nmax,4);%line segment set:[Ax,Ay,Bx,By]
areax = 1;
areay = .5;
dftarea = [axisx-areax,axisy-areay;axisx-areax,axisy;axisx,axisy;axisx,axisy-areay];
okarea = [axisx-areax,axisy-2*areay;axisx-areax,axisy-areay;axisx,axisy-areay;axisx,axisy-2*areay];
fill(okarea(:,1),okarea(:,2),'k')
fill(dftarea(:,1),dftarea(:,2),'w')
stt = true;%true: on, ready to insert; false: off, waiting for the 2nd point
polygon = false;%%
n = 0;%number of segments
%debug%
mdldisp = ['A_1';'A_2';'A_3';'B_1';'B_2';'B_3';'B_4';'B_5';'B_6';'C  ';...
    'D_1';'D_2';'D_3';'D_4';'E  ';'nil';'DK '];
%--%
while n <= nmax
    %display current operation
    if n == 0
        set(msg1,'string',['Insert the first edge. Or click the white '...
            'area and load default data.'],'fontsize',11.5,'horizontal','left');
    else
        if stt
            set(msg1,'string',['Insert a new line: The ',orderdisp(n+1),...
                ' edge. First end point. Or click the black area to finish.'],...
                'fontsize',11.5,'horizontal','left');
        else
            set(msg1,'string',['Insert a new line: The ',orderdisp(n),...
                ' edge. Second end point.'],'fontsize',11.5,...
                'horizontal','left');
        end
    end
    %ginput
    [x,y] = ginput(1);
    if (x-(axisx-areax))*(x-axisx)<0&&(y-(axisy-2*areay))*(y-(axisy-areay))<0 ||...
           ((x-(axisx-areax))*(x-axisx) < 0&&(y-(axisy-areay))*(y-axisy)<0&&n == 0)
        break
    end
    plot(x,y,'ro')
    if stt
        n = n+1;
        segset(n,1:2) = [x,y];
    else
        segset(n,3:4) = [x,y];
        plot(segset(n,[1,3]),segset(n,[2,4]),'b');
    end
    stt = ~stt;
end
if n ~= nmax
    segset = segset(1:n,:);
end

%default map
dftmap = [2,2,2,7;2,7,3,8;3,8,8,8;8,8,8,2;8,2,4,2;4,2,4,6;4,6,6,6;6,6,6,4];
if n == 0
    segset = dftmap;
    n = 8;
end

%Robot
% display current operation
set(msg1,'string','Set original position of the robot. Or click the white area and load default data.',...
    'fontsize',11.5,'horizontal','left');
[rx,ry] = ginput(1);
if (rx-(axisx-areax))*(rx-axisx)<0&&(ry-(axisy-areay))*(ry-axisy)<0
    rx = 5;
    ry = 5;
end
plot(rx,ry,'ro')
a = 2*pi/ns;

%Probing & Moving
set(msg1,'string','Looking for exit...','fontsize',11.5,'horizontal','left');
rtrack = zeros(stpmax,2);rtrack(1,:) = [rx,ry];
stp = 0;
    function robot_move_fcn( ~,evt)
        switch evt.Key
            case {'w','uparrow'}
                ry = ry+ds;
            case {'s','downarrow'}
                ry = ry-ds;
            case {'a','leftarrow'}
                rx = rx-ds;
            case {'d','rightarrow'}
                rx = rx+ds;
        end
        rxo = rx;
        ryo = ry;
        if approachwall(segset,[rx,ry],2*r)
            set(msg1,'string','Approaching wall...',...
                'fontsize',11.5,'horizontal','left');
        else
            set(msg1,'string','Searching exit...',...
                'fontsize',11.5,'horizontal','left');
        end
        [istouch,~] = approachwall(segset,[rx,ry],r);
        if istouch
            set(msg1,'string','Robot has touched wall... Move back.',...
                'fontsize',11.5,'horizontal','left');
            rx = rxo;
            ry = ryo;
        end
        if ~istouch
            stp = stp+1;
            if stp>stpmax
                rtrack = [rtrack;zeros(stpmax,2)];
                stpmax = 2*stpmax;
            end
            rtrack(stp,:) = [rx,ry];
            %reprobe
            rs = sonomeasure(segset,[rx,ry],polygon,ns);
            [modelpt,mptadd,finished,mdl] = sonomodel(rs,r);%!!
            %redraw
            cla
            % environment
            if ~mapoff
                for k = 1:n
                    plot(segset(k,[1,3]),segset(k,[2,4]),'--r','LineWidth',2)
                end
            end
            % probing data
            for k = 1:ns
                plot(rx+rs(k)*cos(((k-1):.1:k)*a),ry+rs(k)*sin(((k-1):.1:k)*a),'b:')
            end
            plot(rx,ry,'ro');%text(rx-.15,ry-.15,'Robot');
            plot(rx+r*cos((0:.1:2)*pi),ry+r*sin((0:.1:2)*pi),'b','LineWidth',1.45);
            % model report
            mdlrpt = [0 0];
            for k = 1:ns
                if finished(k)
                    plot(rx+modelpt(k,[1,3]),ry+modelpt(k,[2,4]),'k','LineWidth',1.5)
                    switch mdl(k)
                        case {1,2,3}
                            plot(rx+(rs(k)-.2)*cos((k-1:.1:k)*a),ry+(rs(k)-.2)*sin((k-1:.1:k)*a),'c','LineWidth',2.5);
                            mdlrpt(1) = 1;
                        case {4,5,6,7,8,9}
                            plot(rx+(rs(k)-.2)*cos((k-1:.1:k)*a),ry+(rs(k)-.2)*sin((k-1:.1:k)*a),'g','LineWidth',2.5);
                            mdlrpt(2) = 1;
                    end
                end
            end
            for k = 1:size(mptadd,1)
                plot(rx+mptadd(k,[1,3]),ry+mptadd(k,[2,4]),'k','LineWidth',1.5)
            end
            if mdlrpt(1)&&mdlrpt(1)
                set(msg2,'string','Faultages found...',...
                    'fontsize',11.5,'horizontal','left');
            elseif mdlrpt(1)
                set(msg2,'string','Faultages on two sides...',...
                    'fontsize',11.5,'horizontal','left');
            elseif mdlrpt(2)
                set(msg2,'string','Faultages on one side...',...
                    'fontsize',11.5,'horizontal','left');
            else
                set(msg2,'string','Measurement results...',...
                    'fontsize',11.5,'horizontal','left');
            end
            % robot track
            plot(rtrack(1:stp,1),rtrack(1:stp,2),'b:')
            %debug%
            % report model
            for k = 1:ns
                text(rx+rs(k)*cos((k-.5)*a),ry+rs(k)*sin((k-.5)*a),mdldisp(mdl(k),:))
            end
        end
    end

end


text:
function [result] = ft_test_report(varargin)

% FT_TEST_REPORT

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

narginchk(1, inf);
command = varargin{1};
assert(isequal(command, 'report'));
varargin = varargin(2:end);

optbeg = find(ismember(varargin, {'matlabversion', 'fieldtripversion', 'user', 'hostname', 'branch', 'arch', 'showdate', 'showid', 'timeout'}));
if ~isempty(optbeg)
  optarg   = varargin(optbeg:end);
  varargin = varargin(1:optbeg-1);
else
  optarg = {};
end

% varargin contains the file (or files) to test
% optarg contains the command-specific options
showdate = ft_getopt(optarg, 'showdate', false);  if ischar(showdate), showdate = istrue(showdate);    end
showid   = ft_getopt(optarg, 'showid', false);    if ischar(showid),   showid   = istrue(showid);      end
timeout  = ft_getopt(optarg, 'timeout', 30);      if ischar(timeout),  timeout  = str2double(timeout); end

% construct the query string that will be passed in the URL
query = '?';
queryparam = {'matlabversion', 'fieldtripversion', 'hostname', 'user', 'branch', 'arch'};
for i=1:numel(queryparam)
  val = ft_getopt(optarg, queryparam{i});
  if ~isempty(val)
    query = [query sprintf('%s=%s&', queryparam{i}, val)];
  end
end

options = weboptions('ContentType', 'json', 'Timeout', timeout); % this returns the result as MATLAB structure
url = 'http://dashboard.fieldtriptoolbox.org/api/';

if isempty(varargin)
  result = webread([url query], options);
  assert(~isempty(result), 'no results were returned');
  result = mergecellstruct(result);
else
  results = cell(size(varargin));
  for i=1:numel(varargin)
    result = webread([url query sprintf('&functionname=%s', varargin{i})], options);
    assert(~isempty(result), 'no results were returned for %s', varargin{i});
    results{i} = mergecellstruct(result);
  end
  % merge all results
  result = mergecellstruct(results);
end

% rename the automatically added fields
result = renamefields(result, 'x_id', 'id');
result = renamefields(result, 'createDate', 'date');

% remove some of the fields
if ~showid
  result = removefields(result, 'id');
end
if ~showdate
  result = removefields(result, 'date');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function fileList = resortFileOrder(fileList)
% this patches the file ordering for matlab, ie. to take into account the re-sorting of length wise as well
% fileList = char(fileList);

sortedT = sort(fileList);

for i = 1:length(sortedT)
   lengthOfEachString(i) = length(char(sortedT(i))) ;
end

sortByLength = sort(lengthOfEachString);

% remove repeated numbers
arrayOfIndices = sortByLength;
uniqueNumbers = unique(arrayOfIndices);
sortByLength = uniqueNumbers;

numberOfGroups = length(sortByLength);

for i = 1:numberOfGroups
   indexToCheckFor = find(sortByLength(i) == lengthOfEachString);
   sortedWithLengthTakenCare{i} = char(sortedT(indexToCheckFor));
end

sortedIndex = [];
% append the clustered list
for i = 1:numberOfGroups
   % fprintf('%s \n', sortedWithLengthTakenCare{i});

   numberOfItemPerCluster = size(sortedWithLengthTakenCare{i}, 1);

   clusterToBeAddressed = cellstr(sortedWithLengthTakenCare{i});

   for j=1:numberOfItemPerCluster
       sortedIndex = [sortedIndex, '', clusterToBeAddressed(j)];
   end
   % sortedIndex = [sortedIndex sortedWithLengthTakenCare{i}];
end

fileList = char(sortedIndex);

end

%{
clc;
clear all;

targetFolder = 'repositoryArchive';
imageFileType = '.jpg';

[imageFileNames, imageFileCount] = listFilesInDirectory( targetFolder, imageFileType );

for i = 1:imageFileCount
    imageChoiceIndex = i;
    imageFileAddressed = cell2str(imageFileNames(imageChoiceIndex));
    % remove '{' and '}'
    imageFileAddressed = regexprep(imageFileAddressed, '{', '');
    imageFileAddressed = regexprep(imageFileAddressed, '}', '');
    imageFileAddressed = regexprep(imageFileAddressed,'''', ''); % remove character '

    numberIndex(i) = str2num(extractNumberIndexFromFileName(imageFileAddressed));
end

sortedFileOrder = sort(numberIndex, 'ascend');

% depicts the file picked in order

for i = 1:imageFileCount
    addressedFile(i) = find (sortedFileOrder(i) == numberIndex);
end
%}


text:
%TRPRINT Compact display of homogeneous transformation
%
% TRPRINT(T, OPTIONS) displays the homogoneous transform in a compact
% single-line format.  If T is a homogeneous transform sequence then each
% element is printed on a separate line.
%
% S = TRPRINT(T, OPTIONS) as above but returns the string.
%
% TRPRINT T  is the command line form of above, and displays in RPY format.
%
% Options::
% 'rpy'        display with rotation in roll/pitch/yaw angles (default)
% 'euler'      display with rotation in ZYX Euler angles
% 'angvec'     display with rotation in angle/vector format
% 'radian'     display angle in radians (default is degrees)
% 'fmt', f     use format string f for all numbers, (default %g)
% 'label',l    display the text before the transform
%
% Examples::
%        >> trprint(T2)
%        t = (0,0,0), RPY = (-122.704,65.4084,-8.11266) deg
%
%        >> trprint(T1, 'label', 'A')
%               A:t = (0,0,0), RPY = (-0,0,-0) deg
%
% See also TR2EUL, TR2RPY, TR2ANGVEC.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function out = trprint(T, varargin)

    if ischar(T)
        % command form: trprint T
        trprint( evalin('base', T) );
        return;
    end

    opt.fmt = [];
    opt.mode = {'rpy', 'euler', 'angvec'};
    opt.radian = false;
    opt.label = '';

    opt = tb_optparse(opt, varargin);

    s = '';

    if size(T,3) == 1
        if isempty(opt.fmt)
            opt.fmt = '%g';
        end
        s = tr2s(T, opt);
    else
        if isempty(opt.fmt)
            opt.fmt = '%8.2g';
        end

        for i=1:size(T,3)
            % for each 4x4 transform in a possible 3D matrix
            s = char(s, tr2s(T(:,:,i), opt) );
        end
    end

    % if no output provided then display it
    if nargout == 0
        disp(s);
    else
        out = s;
    end
end

function s = tr2s(T, opt)
    % print the translational part if it exists
    if ~isempty(opt.label)
        s = sprintf('%8s: ', opt.label);
    else
        s = '';
    end
    if ~isrot(T)
        s = strcat(s, sprintf('t = (%s),', vec2s(opt.fmt, transl(T)')));
    end

    % print the angular part in various representations
    switch (opt.mode)
        case {'rpy', 'euler'}
            % angle as a 3-vector
            if strcmp(opt.mode, 'rpy')
                ang = tr2rpy(T);
                label = 'RPY';
            else
                ang = tr2eul(T);
                label = 'EUL';
            end
            if opt.radian
                s = strcat(s, ...
                    sprintf(' %s = (%s) rad', label, vec2s(opt.fmt, ang)) );
            else
                s = strcat(s, ...
                    sprintf(' %s = (%s) deg', label, vec2s(opt.fmt, ang*180.0/pi)) );
            end
        case 'angvec'
            % as a vector and angle
            [th,v] = tr2angvec(T);
            if isnan(v(1))
                s = strcat(s, sprintf(' R = nil') );
            elseif opt.radian
                s = strcat(s, sprintf(' R = (%sdeg | %s)', ...
                    sprintf(opt.fmt, th), vec2s(opt.fmt, v)) );
            else
                s = strcat(s, sprintf(' R = (%sdeg | %s)', ...
                    sprintf(opt.fmt, th*180.0/pi), vec2s(opt.fmt,v)) );
            end
    end
end

function s = vec2s(fmt, v)
    s = '';
    for i=1:length(v)
        s = strcat(s, sprintf(fmt, v(i)));
        if i ~= length(v)
            s = strcat(s, ', ');
        end
    end
end

author: Allopart

repo: rbpf-gmapping


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: midpoint quadrature rule for 1D spline
%
%==============================================================================

clear, close all, help(mfilename)

omega = [0,6];   I = 1;  psi = @(x) spline1D(3,x); h = []; Q = [];
for j=1:10,
  m    = 2^j;
  h(j) = diff(omega)/m;
  xc   = getCellCenteredGrid(omega,m);
  Q(j) = h(j)*sum(psi(xc));
end;
figure(1); clf; p1=semilogx(h/h(1),Q+eps,'kx',h/h(1),Q,'k-');
figure(2); clf; p2=loglog(h/h(1),abs(I-Q)+eps,'kx',h/h(1),abs(I-Q),'k-');
%==============================================================================