text:
function varargout = CannyEdgeDetector(varargin)
% CannyEdgeDetector Detects edges in the intensity image.
%     imo = CannyEdgeDetector(imi) takes the intensity image imi as input
%     and returns a binary image imo with the same size as imi, where 1's
%     are the edge pixels detected by Canny's method (default). The default
%     value for the thresholds are specified as
%     0.3*max(GradientMagnitudeImage), 0.1*max(GradientMagnitudeImage) and
%     1 for high and low thresholds and sigma for Gaussian kernel,
%     respectively. There are also some additional features listed as
%     below:
%
%     imo = CannyEdgeDetector(imi,method) takes an additional string
%     'method' as input which specifies the function to detect the edges
%     whether using Canny's method, or Laplacian of Gaussian. The possible
%     valid entries for method input are 'canny' and 'log'. The default
%     threshold value for log method is automatically chosen by 'edge'
%     function, which is present in "Image Processing Toolbox".
%
%     imo = CannyEdgeDetector(imi,sigma) returns the same output as
%     described at the top, only with the specified sigma value for
%     Gaussian kernel in stead of the default value. The default value for
%     method is 'canny'.
%
%     imo = CannyEdgeDetector(imi,sigma,method) returns the edge detected
%     binary image imo, with the specified sigma value and the specified
%     method, in stead of the defaults.
%
%     imo = CannyEdgeDetector(...,Th) returns the same output as described
%     above, with an additional input to specify the threshold for 'log'
%     method, or the high threshold for 'canny'. In the second case, the
%     low threshold is 1/3 of Th by default.
%
%     imo = CannyEdgeDetector(...,Th,Tl) takes an additional input Tl,
%     which specifies the low threshold in case of Canny's method. When the
%     method input is 'log', this usage is invalid and the function will
%     return an error message.
%
%     imo = CannyEdgeDetector(...,Th,Tl,feature) In the case of Canny's
%     method, there are two distinct ways to apply hysteresis thresholding
%     defined. The first one is a recursive search (namely, depth first
%     search), and the second one is to use some morphological operations.
%     Obviously, the recursive method takes longer to operate. The string
%     input 'feature' specifies which one to use. Valid options for this
%     input are 'dfs' and 'morph'. It is defined as 'dfs' by default. In
%     the case of log method, this usage is invalid and will return an
%     error message.
%
%     CannyEdgeDetector(...) usage (without the output arguments) plots the
%     binary image using imshow, in stead of making any assignments.
%
%     Author: Halim Can ALBASAN
%
%     Written for the term project of the course "Robot Vision" (EE701)
%     given in Middle East Technical University, Ankara, Turkey; by the
%     lecturer A. Aydin Alatan.
%
%     Date: 15/06/2008

[imi,sigma,method,Th,Tl,feature] = parse_inputs (varargin{:});

switch method
    case 'canny'

%% Detection of Gaussian kernel size ( 3 <= N <= 15 )

G = gaussmf (-10:10,[sigma 0]);
ind = find (G<=.1);
if length(ind)<=6
    N = 15;                             % Upper boundary
    warning('MATLAB:LargeSigma','Sigma too large!')

elseif length(ind)>=18
    N = 3;                              % Lower boundary
    warning('MATLAB:SmallSigma','Sigma too small!')
else
    G(ind) = [];
    N = length (G);
end

%% Gaussian kernel and its first derivatives ( horizontal (Gh) and vertical
%% (Gv) )

G = fspecial ('gaussian', [N N], sigma);
[Gh,Gv] = gradient(G);

%% Applying the operators to the image and obtaining the gradient
%% magnitudes and directions

imx = conv2(double(imi),Gh);
imx = imx ((N-1)/2+1:size(imx,1)-(N-1)/2,(N-1)/2+1:size(imx,2)-(N-1)/2);
imy = conv2(double(imi),Gv);
imy = imy ((N-1)/2+1:size(imy,1)-(N-1)/2,(N-1)/2+1:size(imy,2)-(N-1)/2);

gra_mag = sqrt(imx.^2+imy.^2);
gra_dir = atan2(imy,imx);

%% Gradient direction discretization

[r c] = size (gra_dir);

for i=1:r
    for j=1:c
        if ( gra_dir(i,j)<pi/8 && gra_dir(i,j)>=0 ) ...
                || ( gra_dir(i,j)<=2*pi && gra_dir(i,j)>=15*pi/8 ) ...
                || ( gra_dir(i,j)<9*pi/8 && gra_dir(i,j)>=7*pi/8 )
            gra_dir(i,j) = 0;
        elseif ( gra_dir(i,j)>=pi/8 && gra_dir(i,j)<3*pi/8 ) ...
                || ( gra_dir(i,j)>=9*pi/8 && gra_dir(i,j)<11*pi/8 )
            gra_dir(i,j) = 45;
        elseif ( gra_dir(i,j)>=3*pi/8 && gra_dir(i,j)<5*pi/8 ) ...
                || ( gra_dir(i,j)>=11*pi/8 && gra_dir(i,j)<13*pi/8)
            gra_dir(i,j) = 90;
        else
            gra_dir(i,j) = 135;
        end
    end
end

%% Nonmaxima suppression

gra_mag_s = gra_mag;
gra_mag_s(1,:) = zeros (1,c);
gra_mag_s(r,:) = zeros (1,c);
gra_mag_s(:,1) = zeros (r,1);
gra_mag_s(:,c) = zeros (r,1);             % Suppress boundaries first

for i=2:r-1
    for j=2:c-1
        if gra_dir(i,j)==0
            if ( gra_mag(i,j)<=gra_mag(i,j+1) ) ...
                    || ( gra_mag(i,j)<=gra_mag(i,j-1) )
                gra_mag_s(i,j) = 0;
            end
        elseif gra_dir(i,j)==45
            if ( gra_mag(i,j)<=gra_mag(i-1,j+1) ) ...
                    || ( gra_mag(i,j)<=gra_mag(i+1,j-1) )
                gra_mag_s(i,j) = 0;
            end
        elseif gra_dir(i,j)==90
            if ( gra_mag(i,j)<=gra_mag(i+1,j) ) ...
                    || ( gra_mag(i,j)<=gra_mag(i-1,j) )
                gra_mag_s(i,j) = 0;
            end
        else
            if ( gra_mag(i,j)<=gra_mag(i+1,j+1) ) ...
                    || ( gra_mag(i,j)<=gra_mag(i-1,j-1) )
                gra_mag_s(i,j) = 0;
            end
        end
    end
end

%% Hysteresis thresholding and forming the output image

if isempty(Th)
    Tl = .10*max(max(gra_mag_s));
    Th = .30*max(max(gra_mag_s));
end

h_thr_im = gra_mag_s;
h_thr_im(gra_mag_s<Th) = 0;

l_thr_im = gra_mag_s;
l_thr_im(gra_mag_s<Tl) = 0;

        switch feature
            case 'morph'

h_thr_im = logical(h_thr_im);
l_thr_im = logical(l_thr_im);

[ii jj] = find(h_thr_im);

imo = bwselect(l_thr_im , jj , ii , 8);
imo = bwmorph(imo , 'thin' , 1);

            case 'dfs'

for i=2:r-1
    for j=2:c-1
        if h_thr_im (i,j)>0
            h_thr_im = edge_follow(h_thr_im,l_thr_im,i,j);
        end
    end
end

imo = logical(h_thr_im);

        end

    case 'log'

imo = edge(imi,method,Th,sigma);

end

if nargout<1
    imshow(imo)
else
    varargout{1} = imo;
end

function thr_h = edge_follow(thr_h,thr_l,i,j)

    x = [-1  0  1 -1 1 -1 0 1];               % Relative coordinates of 8 neighbors
    y = [-1 -1 -1  0 0  1 1 1];

    for k=1:8
        if thr_h(i+x(k),j+y(k))==0 && thr_l(i+x(k),j+y(k))~=0
            thr_h(i+x(k),j+y(k)) = -thr_l(i+x(k),j+y(k));
            thr_h = edge_follow(thr_h,thr_l,i+x(k),j+y(k));
        end
    end

end

function [imi,sigma,method,Th,Tl,feature] = parse_inputs(varargin)

    imi = varargin{1};
    sigma = 1;                  % defaults
    method = 'canny';
    Th = [];                    % the defaults for the thresholds will be
    Tl = [];                    % specified after the gradient magnitude is
    feature = 'dfs';            % calculated.

    methods = {'canny','log'};
    features = {'dfs','morph'};
    sigma_default = false;

    if nargin<1
        error('Not enough input arguments.')
    end

    if nargin>=2

        if ischar(varargin{2})
            if length(strmatch(varargin{2},methods))~=1
                error('Wrong input for edge detection method')
            else
                method=varargin{2};
                sigma_default = true;
            end

        else
            sigma=varargin{2};
        end
    end

    if nargin>=3

        if sigma_default
            Th = varargin{3};
            Tl = Th/3;
        else
            if length(strmatch(varargin{3},methods))~=1
                error('Wrong input for edge detection method')
            else
                method=varargin{3};
            end
        end
    end

    if nargin>=4

        if sigma_default
            if strcmp(method,'canny')
                Tl = varargin{4};
            else
                error('Wrong sequence or number of inputs for "log" method')
            end
        else
            Th = varargin{4};
        end
    end

    if nargin>=5
        if strcmp(method,'log')
            error('Too many input arguments!')
        else
            if sigma_default
                if ischar(varargin{5})
                    if length(strmatch(varargin{5},features))~=1
                        error('Wrong input for hysteresis thresholding method')
                    else
                        feature = varargin{5};
                    end
                else
                    error('Wrong usage of input argument "feature"!')
                end
            else
                Tl = varargin{5};
            end
        end
    end

    if nargin>=6
        if sigma_default
            error('Too many input arguments!')
        else
            if ischar(varargin{6})
                if length(strmatch(varargin{6},features))~=1
                    error('Wrong input for hysteresis thresholding method')
                else
                    feature = varargin{5};
                end
            else
                error('Wrong usage of input argument "feature"!')
            end
        end
    end

    if nargin>6
        error('Too many input arguments')
    end
end

end


%% Example script to visualize the aircraft simulation data
% Add the path of the aircraft_3d_animation function
addpath('../src/');
% path of the *.mat file containing the 3d model information
model_info_file = '../3d_models/x15_3d_model.mat';
% Load the simulation data
% load('scissors_maneuver.mat')
% load('breakaway_maneuver.mat')
% load('split_s_maneuver.mat')
% load('departure.mat')
load('departure2.mat')
% define the reproduction speed factor
speedx = 1;
% Do you want to save the animation in a mp4 file? (0.No, 1.Yes)
isave_movie = 1;
% Movie file name
movie_file_name = 'departure2.mp4';

% -------------------------------------------------------------------------
% The frame sample time shall be higher than 0.02 seconds to be able to
% update the figure (CPU/GPU constraints)
frame_sample_time = max(0.02, tout(2)-tout(1));
% Resample the time vector to modify the reproduction speed
t_new   = tout(1):frame_sample_time*(speedx):tout(end);
% Resample the recorded data
act     = interp1(tout, act, t_new','linear');
stick   = interp1(tout, stick, t_new','linear');
y_new   = interp1(tout, yout, t_new','linear');
% We have to be careful with angles with ranges
y_new(:, 7)  = atan2(interp1(tout, sin(yout(:, 7)), t_new','linear'), interp1(tout, cos(yout(:, 7)), t_new','linear')) * 180 / pi;
y_new(:, 8)  = atan2(interp1(tout, sin(yout(:, 8)), t_new','linear'), interp1(tout, cos(yout(:, 8)), t_new','linear')) * 180 / pi;
y_new(:, 9)  = atan2(interp1(tout, sin(yout(:, 9)), t_new','linear'), interp1(tout, cos(yout(:, 9)), t_new','linear')) * 180 / pi;
% Assign the data
heading_deg           =  y_new(:, 7);
pitch_deg             =  y_new(:, 8);
bank_deg              =  y_new(:, 9);
roll_command          = -stick(:, 2);
pitch_command         = -stick(:, 1);
angle_of_attack_deg   =  y_new(:, 2) * 180 / pi;
angle_of_sideslip_deg =  y_new(:, 3) * 180 / pi;
fligh_path_angle_deg  =  y_new(:, 22) * 180 / pi;
mach                  =  y_new(:, 21);
altitude_ft           = -y_new(:, 12);
nz_g                  =  y_new(:, 19);
% Flight control surfaces
dr     = -act(:, 8);
df1    =  act(:, 6);
df2    =  act(:, 5);
df3    =  act(:, 4);
df4    =  act(:, 3);
% Control array assignation
% (modify the order according to your particular 3D model)
controls_deflection_deg = [ 0.25 * (df1(:) + df2(:) - df3(:) - df4(:)), ...
                           -0.25 * (df1(:) + df2(:) - df3(:) - df4(:)), ...
                            0.25 * (df1(:) + df2(:) + df3(:) + df4(:)), ...
                            0.25 * (df1(:) + df2(:) + df3(:) + df4(:)), ...
                            df1(:) * 0, ...
                            df1(:) * 0, ...
                            dr(:)];

%% Run aircraft_3d_animation function
% -------------------------------------------------------------------------
aircraft_3d_animation(model_info_file,...
    heading_deg, ...            Heading angle [deg]
    pitch_deg, ...              Pitch angle [deg]
    bank_deg, ...               Roll angle [deg]
    roll_command, ...           Roll  stick command [-1,+1] [-1 -> left,            +1 -> right]
    pitch_command, ...          Pitch stick command [-1,+1] [-1 -> full-back stick, +1 -> full-fwd stick]
    angle_of_attack_deg, ...    AoA [deg]
    angle_of_sideslip_deg, ...  AoS [deg]
    fligh_path_angle_deg, ...   Flight path angle [deg]
    mach, ...                   Mach number
    altitude_ft, ...            Altitude [ft]
    nz_g,  ...                  Vertical load factor [g]
    controls_deflection_deg, ...Flight control deflection (each column is a control surface)
    frame_sample_time, ...      Sample time [sec]
    speedx, ...                 Reproduction speed
    isave_movie, ...            Save the movie? 0-1
    movie_file_name);           % Movie file name


text:
function viewFourthOrderTensor(varargin)

% function viewFourthOrderTensor(C,numDigits,fontSizeIm,fontSize)

%% Parse input

numDigits=[];
fontSizeIm=[];
fontSize=[];
toleranceLevel=eps(1);
switch nargin
    case 1
        C=varargin{1};
        toleranceLevel=eps(1);
    case 2
        C=varargin{1};
        numDigits=varargin{2};
    case 3
        C=varargin{1};
        numDigits=varargin{2};
        fontSizeIm=varargin{3};
    case 4
        C=varargin{1};
        numDigits=varargin{2};
        fontSizeIm=varargin{3};
        fontSize=varargin{4};
    case 5
        C=varargin{1};
        numDigits=varargin{2};
        fontSizeIm=varargin{3};
        fontSize=varargin{4};
        toleranceLevel=varargin{5};
end

if isempty(C)
    I=eye(3,3); %The 2nd order identity tensor
    C=dyadicProduct(I,I,1); %A type of 4th order identity tensor
end

if isempty(numDigits)
    numDigits=3;
end

if isempty(fontSizeIm)
    fontSizeIm=15;
end

if isempty(fontSize)
    fontSize=15;
end

%%
% Plot settings
faceAlpha=0.5;

%%

ind_C_all=1:numel(C);
[I,J,K,L]=ind2sub(size(C),ind_C_all);

kronD_IJ=(I==J);
kronD_KL=(K==L);

switch class(C)
    case 'sym'
        CV=sym(zeros(6,6));
        logicSym=1;
    otherwise
        CV=zeros(6,6);
        logicSym=0;
end

%Create mapping indices
P=I.*kronD_IJ+(1-kronD_IJ).*(9-I-J);
Q=K.*kronD_KL+(1-kronD_KL).*(9-K-L);

%Convert to linear indices
[ind_CV]=sub2ind(size(CV),P,Q);

%Get the unique 36 entries
[~,indUni,~]=unique(ind_CV,'first');
ind_CV_uni=ind_CV(indUni);

ind_C_uni=ind_C_all(indUni);

CV(ind_CV_uni)=C(ind_C_uni);

%%

%Derive 9x9 matrix representation
[CM]=fourthOrderMat(C);

if logicSym
    symbolicVariableSet=symvar(CM);
    C_dummy=double(subs(CM,symbolicVariableSet,1:numel(symbolicVariableSet)));
    logicEntry=abs(C_dummy)>toleranceLevel;
    C_dummy_V=double(subs(CV,symbolicVariableSet,1:numel(symbolicVariableSet)));
    logicEntry_V=abs(C_dummy_V)>toleranceLevel;
else
    logicEntry=abs(CM)>toleranceLevel;
    logicEntry_V=abs(CV)>toleranceLevel;
end

%%

C_V=zeros(size(C));
C_V(ind_C_uni)=ind_CV_uni;
[CM_V]=fourthOrderMat(C_V);
logicEntry_CV=CM_V>eps(CM_V);

%%
hf1=cFigure;
title('9x9 array mapping of $4^{th}$ order tensor','fontSize',fontSize,'Interpreter','LATEX');
xlabel('$q=3(j-1)+l$','fontSize',fontSize,'Interpreter','LATEX');
ylabel('$p=3(i-1)+k$','fontSize',fontSize,'Interpreter','LATEX');

if logicSym
    [Fp,Vp,Cp]=ind2patch(logicEntry,C_dummy,'sk');
    gpatch(Fp,Vp,Cp,0.25.*ones(1,3),faceAlpha);
else
    [Fp,Vp,Cp]=ind2patch(logicEntry,CM,'sk');
    gpatch(Fp,Vp,Cp,0.25.*ones(1,3),faceAlpha);
end
colormap gjet;
if ~logicSym
    colorbar;
end

m=zeros(3,3);
[Fp,Vp,~]=ind2patch(true(size(m)),m,'sk');
Vp(:,[1 2])=3*Vp(:,[1 2])-1;
patch('Faces',Fp,'Vertices',Vp,'FaceColor','none','EdgeColor','k','lineWidth',5);

[Fp,Vp,~]=ind2patch(logicEntry_CV,logicEntry_CV,'sk');
patch('Faces',Fp,'Vertices',Vp,'FaceColor',0.25.*ones(1,3),'EdgeColor','k','lineWidth',1,'FaceAlpha',0.25);

image_numeric(CM,hf1,numDigits,fontSizeIm,'k','tex');

view(2); axis tight; axis ij; axis square;
set(gca,'FontSize',fontSize);
drawnow;

%%

hf2=cFigure;
title('Voigt array mapping of $4^{th}$ order tensor','fontSize',fontSize,'Interpreter','LATEX');
xlabel('$q=k\delta_{kl}+(1-\delta_{kl})(9-k-l)$','fontSize',fontSize,'Interpreter','LATEX');
ylabel('$p=i\delta_{ij}+(1-\delta_{ij})(9-i-j)$','fontSize',fontSize,'Interpreter','LATEX');

if logicSym
    [Fp,Vp,Cp]=ind2patch(logicEntry_V,C_dummy_V,'sk');
    gpatch(Fp,Vp,Cp,0.25.*ones(1,3),faceAlpha);
else
    [Fp,Vp,Cp]=ind2patch(logicEntry_V,CV,'sk');
    gpatch(Fp,Vp,Cp,0.25.*ones(1,3),faceAlpha);
end
colormap gjet;
if ~logicSym
    colorbar;
end

m=zeros(2,2);
[Fp,Vp,~]=ind2patch(true(size(m)),m,'sk');
Vp(:,[1 2])=3*Vp(:,[1 2])-1;
patch('Faces',Fp,'Vertices',Vp,'FaceColor','none','EdgeColor','k','lineWidth',5);

image_numeric(CV,hf2,numDigits,fontSizeIm,'k','tex');
view(2); axis tight; axis ij; axis square;
set(gca,'FontSize',fontSize);
drawnow;

%%
% _*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 [lo hi] = ir_dwt_filters(wname, varargin)
%|function [coef codes] = ir_odwt1(x, varargin)
%|
%| filters for discrete wavelet transform DWT
%|
%| in
%|	'wname'	char	'haar' (default) or 'db4' or 'sym2'
%|
%| option
%|	'usemat' 0|1	default: 0. (if 1 then use matlab wfilters)
%|	'ortho'	0|1	default: 1 (make normalized)
%|	'abs'	0|1	default: 0. (if 1 then take abs of filters)
%|
%| out
%|	lo	[K]	low-pass decomposition filter
%|	hi	[K]	hi-pass decomposition filter
%|
%| 2012-05-21, Jeff Fessler, Univ. of Michigan
%| 2014-01-23, Matt Muckley: added 'db4'

if nargin < 1, ir_usage, end
if streq(wname, 'test'), ir_dwt_filters_test, return, end

arg.wname = wname;
arg.ortho = true;
arg.usemat = false;
arg.abs = false;
arg = vararg_pair(arg, varargin);

if isempty(arg.wname)
	 arg.wname = 'haar';
end

% decomposition filters
if arg.usemat
	[lo hi lo_r hi_r] = wfilters(arg.wname);
	lo = lo'; lo_r = lo_r';
	hi = hi'; hi_r = hi_r';

	jf_equal(lo_r, flipud(lo))
	jf_equal(hi_r, flipud(hi))
else
	switch arg.wname
	case 'haar'
		lo = [1 1]';
		hi = [-1 1]';

	case 'db4'
		h1 = (1 + sqrt(3))/4;
		h2 = (3 + sqrt(3))/4;
		h3 = (3 - sqrt(3))/4;
		h4 = (1 - sqrt(3))/4;
		lo = [h1 h2 h3 h4]';
		hi = [h4 -h3 h2 -h1]';
		clear h1 h2 h3 h4;

	case 'sym2'
		lo = [-0.129409522550921; 0.224143868041857; ...
			0.836516303737469; 0.482962913144690];
		hi = [-0.482962913144690; 0.836516303737469; ...
			-0.224143868041857; -0.129409522550921];
	otherwise
		fail('unknown wname "%s"', arg.wname)
	end
end

if arg.ortho
	lo = lo / norm(lo);
	hi = hi / norm(hi);
end

if arg.abs
	lo = abs(lo);
	hi = abs(hi);
end

% ir_dwt_filters_test()
function ir_dwt_filters_test

list = {'haar', 'sym2', 'db4'};
for ii=1:numel(list)
	wname = list{ii};

	[lo0 hi0] = ir_dwt_filters(wname, 'usemat', 0);
	if numel(lo0) == 2
		jf_equal(lo0' * hi0, 0)
	else
		equivs(1, 1 + lo0' * hi0)
	end

	if exist('dwt', 'file') == 2 && exist('wfilters', 'file') == 2
		[lo1 hi1] = ir_dwt_filters(wname, 'usemat', 1);
	end

	[lo hi] = ir_dwt_filters(wname, 'ortho', 0);
	jf_equal(lo/norm(lo), lo0)
	jf_equal(hi/norm(hi), hi0)

	[lo hi] = ir_dwt_filters(wname, 'abs', 1);
	jf_equal(abs(lo0), lo)
	jf_equal(abs(hi0), hi)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [cl,clpos,clneg,clrpos,clrneg] = cluster_sigregions(cl,pthr,varargin)
% [cl,clpos,clneg,clrpos,clrneg] = cluster_ttest(clusters,pthr,[behavioral regressor])
% tor wager  Oct 9, 2003
%
% does t-tests on cluster timeseries and all voxels
% which should contain group data (e.g., con* data)
%
% if optional behavioral regressor is entered,
% computes cluster extent nonparametric p-values, and reports
% extent npm p-values for intercept and behavioral reg.
%
% writes output to cluster_ttest_output.txt

fid = fopen('cluster_ttest_output.txt','a');

if length(varargin) > 0, beh = varargin{1};,end

for i = 1:length(cl)

% df = n - k    primary threshold
tthr1 = abs(tinv(pthr,size(cl(i).all_data,1) - size(cl(i).ttest.t,1)));

% ----------------------------------------
% * do nonparametric permutation
% ----------------------------------------

if ~isfield(cl(i),'npm_ttest'),cl(i).npm_ttest = [];,end
if isfield(cl(i).npm_ttest,'tmax') & isfield(cl(i).npm_ttest,'textmax')
    tthr = cl(i).npm_ttest.tthr ;
    rthr = cl(i).npm_ttest.rthr ;
    tnegthr = cl(i).npm_ttest.tnegthr ;
    rnegthr = cl(i).npm_ttest.rnegthr ;
    tmax = cl(i).npm_ttest.tmax ;
    tmin = cl(i).npm_ttest.tmin ;
    rmax = cl(i).npm_ttest.rmax ;
    rmin = cl(i).npm_ttest.rmin ;
    textmax = cl(i).npm_ttest.textmax ;
    rextmax = cl(i).npm_ttest.rextmax ;
else

[tthr,rthr,tnegthr,rnegthr,tmax,rmax,tmin,rmin,textmax,rextmax] = ...
    npm_ttest(cl(i).all_data,1000,beh,pthr,cl(i).XYZ);
close

cl(i).npm_ttest.tthr = tthr;
cl(i).npm_ttest.rthr = rthr;
cl(i).npm_ttest.tnegthr = tnegthr;
cl(i).npm_ttest.rnegthr = rnegthr;
cl(i).npm_ttest.tmax = tmax;
cl(i).npm_ttest.rmax = rmax;
cl(i).npm_ttest.textmax = textmax;
cl(i).npm_ttest.rextmax = rextmax;
cl(i).npm_ttest.tmin = tmin;
cl(i).npm_ttest.rmin = rmin;

end

fprintf(fid,'\nCl.\tVoxels\tThreshold\tMask\tData\t\n');
fprintf(fid,'\tx\ty\tz\tVoxels\tCluster corr. p.\tMax. t\tCorrected p\tDirection of effect\tr\n');
fprintf(fid,'%3.0f\t%3.0f\t%3.4f\t%s\t%s\t\n',i,cl(i).numVox,pthr,cl(i).P,cl(i).imP(1,:));

% ----------------------------------------
% * get sig clusters
% ----------------------------------------

clpos(i).from_cluster = i;
clpos(i).M = cl(i).M;
clpos(i).voxSize = cl(i).voxSize;
clpos(i).threshold = tthr1;
clpos(i).title = 'Positive effects';

clrpos(i).from_cluster = i;
clrpos(i).M = cl(i).M;
clrpos(i).voxSize = cl(i).voxSize;
clrpos(i).threshold = tthr1;
clrpos(i).title = 'Positive correlations';

sig = cl(i).ttest.t > tthr1;

if any(sig(1,:))
    clpos(i).XYZ = cl(i).XYZ(:,find(sig(1,:)));
    clpos(i).XYZmm = cl(i).XYZmm(:,find(sig(1,:)));
    clpos(i).Z = cl(i).Z(:,find(sig(1,:)));
    clpos(i).t = cl(i).ttest.t(:,find(sig(1,:)));
    gopos(i) = 1;

    % print table entry
    wcl = spm_clusters(clpos(i).XYZ);
    clpos(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clpos(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clpos(i).XYZmm(:,wh)';,end
        nv = length(wh);

        if length(wh)>1,max_t = max(clpos(i).t(:,wh)'); max_t = max_t(1);,
        else, max_t = clpos(i).t(1,wh);
        end

        clpos(i).center = xyz;

        extcor_p = 1 - (sum(textmax <= nv) ./ length(textmax));
        max_p = 1 - (sum(tmax <= max_t) ./ length(tmax));

        fprintf(fid,'\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t+\n', ...
            xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p)
    end

else
    gopos(i) = 0;
end

if any(sig(2,:))
    clrpos(i).XYZ = cl(i).XYZ(:,find(sig(2,:)));
    clrpos(i).XYZmm = cl(i).XYZmm(:,find(sig(2,:)));
    clrpos(i).Z = cl(i).Z(:,find(sig(2,:)));
    clrpos(i).t = cl(i).ttest.t(:,find(sig(2,:)));
    gorpos(i) = 1;

    % print table entry
    wcl = spm_clusters(clrpos(i).XYZ);
    clrpos(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clrpos(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clrpos(i).XYZmm(:,wh)';,end
        nv = length(wh);

        tmp = cl(i).all_data(:,find(sig(2,:)));
        tmp = tmp(:,wh);
        tmp = corrcoef([tmp beh]);
        tmp = tmp(end,1:end-1);
        max_r = max(tmp);
        if max_r > .99, warning('problem?'),keyboard,end

        if length(wh)>1,max_t = max(clrpos(i).t(:,wh)'); max_t = max_t(2);,
        else, max_t = clrpos(i).t(2,wh);
        end
        clrpos(i).center = xyz;

        extcor_p = 1 - sum(rextmax <= nv) ./ length(rextmax);
        max_p = 1 - sum(rmax <= max_t) ./ length(rmax);

        fprintf(fid,'\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t+\t%3.2f\n', ...
            xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p,max_r);
    end

else
    gorpos(i) = 0;
end

% negative

clneg(i).from_cluster = i;
clneg(i).M = cl(i).M;
clneg(i).voxSize = cl(i).voxSize;
clneg(i).threshold = tthr1;
clneg(i).title = 'Negative effects';

clrneg(i).from_cluster = i;
clrneg(i).M = cl(i).M;
clrneg(i).voxSize = cl(i).voxSize;
clrneg(i).threshold = tthr1;
clrneg(i).title = 'Negative correlations';

sig = cl(i).ttest.t < -tthr1;

if any(sig(1,:))
    clneg(i).XYZ = cl(i).XYZ(:,find(sig(1,:)));
    clneg(i).XYZmm = cl(i).XYZmm(:,find(sig(1,:)));
    clneg(i).Z = cl(i).Z(:,find(sig(1,:)));
    clneg(i).t = cl(i).ttest.t(:,find(sig(1,:)));
    goneg(i) = 1;

    % print table entry
    wcl = spm_clusters(clneg(i).XYZ);
    clneg(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clneg(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clpos(i).XYZmm(:,wh)';,end
        nv = length(wh);

        if length(wh)>1,max_t = min(clneg(i).t(:,wh)'); max_t = max_t(1);,
        else, max_t = clneg(i).t(1,wh);
        end

        clneg(i).center = xyz;

        extcor_p = 1 - sum(textmax <= nv) ./ length(textmax);
        max_p = 1 - sum(tmin >= max_t) ./ length(tmin);

        fprintf(fid,'\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t-\n', ...
            xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p);
    end

else
    goneg(i) = 0;
end

if any(sig(2,:))
    clrneg(i).XYZ = cl(i).XYZ(:,find(sig(2,:)));
    clrneg(i).XYZmm = cl(i).XYZmm(:,find(sig(2,:)));
    clrneg(i).Z = cl(i).Z(:,find(sig(2,:)));
    clrneg(i).t = cl(i).ttest.t(:,find(sig(2,:)));
    gorneg(i) = 1;

    % print table entry
    wcl = spm_clusters(clrneg(i).XYZ);
    clrneg(i).clusters = wcl;
    for j = 1:max(wcl)

        wh = find(wcl==j);

        % corrected p

        xyz = mean(clrneg(i).XYZmm(:,wh),2)';
        if length(xyz) == 1, xyz = clrneg(i).XYZmm(:,wh)';,end
        nv = length(wh);

        if length(wh)>1,max_t = min(clrneg(i).t(:,wh)'); max_t = max_t(2);,
        else, max_t = clrneg(i).t(2,wh);
        end

        tmp = cl(i).all_data(:,find(sig(2,:)));
        tmp = tmp(:,wh);
        tmp = corrcoef([tmp beh]);
        tmp = tmp(end,1:end-1);
        max_r = min(tmp);
        if max_r > .99, warning('problem?'),keyboard,end

        clrneg(i).center = xyz;

        extcor_p = 1 - sum(rextmax <= nv) ./ length(rextmax);
        max_p = 1 - sum(rmin >= max_t) ./ length(rmin);

        fprintf(fid,'\t%3.0f\t%3.0f\t%3.0f\t%3.0f\t%3.4f\t%3.2f\t%3.4f\t-\t%3.2f\t\n', ...
            xyz(1),xyz(2),xyz(3),nv,extcor_p,max_t,max_p,max_r);
    end

else
    gorneg(i) = 0;
end

end % loop through clusters

fprintf(fid,'\n')

% Make montage
try

if ~isempty(clpos) & ~isempty(clneg)
    montage_clusters([],cl(unique(cat(2,[clpos.from_cluster clneg.from_cluster]))),clpos,clneg,{'k' 'r' 'b'},'nooverlap')
    %hh=get(gcf,'Children');
    %for i = 1:length(hh),axes(hh(i));, h = findobj(gca,'FaceColor',[0 0 0]);, set(h,'FaceAlpha',.5) ,end
elseif ~isempty(clpos)
    montage_clusters([],cl(cat(2,clpos.from_cluster)),clpos,{'k' 'r'},'nooverlap')
elseif ~isempty(clneg)
    montage_clusters([],cl,clneg,{'k' 'b'},'nooverlap')
else
    disp('No main contrast effects')
end

if ~isempty(clrpos) & ~isempty(clrneg)
    montage_clusters([],cl(unique(cat(2,[clrpos.from_cluster clrneg.from_cluster]))),clrpos,clrneg,{'k' 'r' 'b'},'nooverlap')
elseif ~isempty(clrpos)
    montage_clusters([],cl(cat(2,clrpos.from_cluster)),clrpos,{'k' 'r'},'nooverlap')
elseif ~isempty(clrneg)
    montage_clusters([],cl(cat(2,clrneg.from_cluster)),clrneg,{'k' 'b'},'nooverlap')
else
    disp('No correlation effects')
end

catch
    disp('error with montage')
end

fclose(fid);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function mv=v_roteucode(m)
%V_ROTEUCODE decodes a string specifying a rotation axis sequence
%     M(n)     a string of n characters from the set determining the order of rotation axes
%              as listed below. Note that the control characters 'rdoOaA' may occur anywhere in the string:
%                'x','y','z'    rotate around the given axis by the corresponding angle
%                               given in e()
%                '1','2','3'    90 degree rotation around x,y or z axis; doesn't use a value from e()
%                '4','5','6'    180 degree rotation around x,y or z axis; doesn't use a value from e()
%                '7','8','9'    270 degree rotation around x,y or z axis; doesn't use a value from e()
%                'r','d'        all angles are given in radians or degrees  [radians]
%             'o','O','a','A'   selects whether to rotate the object or the coordinate axes and
%                               whether the rotation axes remain fixed in space for consecutive
%                               rotations (extrinsic) or else move with each rotation (intrinsic).
%                                  'o' = object-extrinsic [default]
%                                  'O' = object-intrinsic
%                                  'a' = axes-extrinsic
%                                  'A' = axes-intrinsic
% Outputs:
%
%     mv(7,k)    where k-1 is the number of non-control characters in the input string m
%                    mv(1,j) = Code for the j'th rotation: 1,2,3 for x,y,z rotation and 4 to 12 for the fixed rotations listed above.
%                              All entries are in the range [1,12] except for mv(1,k)=0.
%                    mv(2,j) = index into euler angle array for x,y,z rotations. mv(2,k) gives total number of euler angles needed
%                    mv(3,j) = rotation class before rotation j. mv(3,k) is the final rotation class and equals 52 for arbitrary rotations.
%                    mv(4,j) = index into a vectorized matrix of the entry that becomes non-zero after rotation j
%                    mv(5,j) = index into a vectorized matrix of the other changing element in the same column
%                    mv(6,j) = +-1 = sign of the sine term affecting entry mv(4,j). For mv(1,j) in [1,3], mv(6,j)=0 if the rotation is unnecessary
%                    mv(7,j) = +-1 = sign of entry mv(5,j) before rotation j if known
%                Special entries:
%                    mv(7,k) = -1 to invert the rotation (i.e. transpose the matrix) or +1 otherwise
%                    mv(4,k) = scale factor for euler angles: +-1 or +-pi/180
%
%
% The string M specifies the seqeunce of axes about which the rotations are performed. There are 12
% possible 3-character sequences that avoid consecutive repetitions. These are 'Euler angles' if
% there is a repeated axis or 'Tait-Bryan angles' if not. Common choices are:
% (1) 'zxz' the most common Euler angle set
% (2) 'xyz' corresponds to 'roll, pitch, yaw' for an aeroplane heading in the x direction with y to
%     the right and z down. The intrinsic equivalent is 'Ozyx' corresponding to 'yaw, pitch, roll'.
% (3) 'z1z1z' involves 5 rotations, in which all the non-fixed rotations are around the z axis.
%

%      Copyright (C) Mike Brookes 2007-2020
%      Version: $Id: v_roteucode.m 11260 2020-07-18 20:07:58Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

persistent mes trmap zel mch mvch  jch nch
if isempty(mes)     % setup fixed arrays and initialize cache of mode strings
    nch=5;          % size of cache
    mch=cell(nch,1); % cache of input character strings
    mvch=cell(nch,1);  % cache of output mv codes
    flefch=zeros(nch,2);  % cache of output flef codes
    jch=(1:nch); % cache usage order jch(1) is the most recent, jch(nch) the oldest
    for i=1:nch
        mch{i}='';
        mvch{i}=[0;0;1;1;0;0;1];
    end
    mes=[1:3 10:12 7:9 4:6]; % sign reversal look-up table
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The trmap and zel arrays contain information about each of 52 different  %
    % patterns of -1,0,+1 that may exist in a rotation matrix as follows:      %
    %                                                                          %
    %   1-3 : identity matrix rows in order: 123, 231, 312                     %
    %   4-6 : negated identity matrix rows in order: 132, 213, 321             %
    %   7-12: As 1-6 but with rows 2,3 negated                                 %
    %  13-18: As 1-6 but with rows 1,3 negated                                 %
    %  19-24: As 1-6 but with rows 1,2 negated                                 %
    %  25-33: +1 in position (i-24) and 0's in remainder of this row and col   %
    %  34-42: -1 in position (i-24) and 0's in remainder of this row and col   %
    %  43-51: 0 in position (i-42)                                             %
    %  52: no special symmetry                                                 %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % trmap(i,j) gives the pattern that i is transformed into by rotation j    %
    % where j=1:3 corresponds to x,y,z and j=4:12 corresponds to the 9         %
    % multiples of 90 degree rotations listed in the main comments.            %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    trmap=[ 25 29 33 16 24 11  7 13 19 22 12 17;
        28 32 27 17 22 12  8 14 20 23 10 18;
        31 26 30 18 23 10  9 15 21 24 11 16;
        34 41 39 13 20  9 10 16 22 19  8 15;
        37 35 42 14 21  7 11 17 23 20  9 13;
        40 38 36 15 19  8 12 18 24 21  7 14;
        25 38 42 22  6 23  1 19 13 16 18  5;
        28 41 36 23  4 24  2 20 14 17 16  6;
        31 35 39 24  5 22  3 21 15 18 17  4;
        34 32 30 19  2 21  4 22 16 13 14  3;
        37 26 33 20  3 19  5 23 17 14 15  1;
        40 29 27 21  1 20  6 24 18 15 13  2;
        34 29 42 10 12  5 19  1  7  4 24 23;
        37 32 36 11 10  6 20  2  8  5 22 24;
        40 26 39 12 11  4 21  3  9  6 23 22;
        25 41 30  7  8  3 22  4 10  1 20 21;
        28 35 33  8  9  1 23  5 11  2 21 19;
        31 38 27  9  7  2 24  6 12  3 19 20;
        34 38 33  4 18 17 13  7  1 10  6 11;
        37 41 27  5 16 18 14  8  2 11  4 12;
        40 35 30  6 17 16 15  9  3 12  5 10;
        25 32 39  1 14 15 16 10  4  7  2  9;
        28 26 42  2 15 13 17 11  5  8  3  7;
        31 29 36  3 13 14 18 12  6  9  1  8;
        25 44 45 25 36 26 25 34 34 25 27 35;
        43 26 45 27 26 34 35 26 35 36 26 25;
        43 44 27 35 25 27 36 36 27 26 34 27;
        28 47 48 28 39 29 28 37 37 28 30 38;
        46 29 48 30 29 37 38 29 38 39 29 28;
        46 47 30 38 28 30 39 39 30 29 37 30;
        31 50 51 31 42 32 31 40 40 31 33 41;
        49 32 51 33 32 40 41 32 41 42 32 31;
        49 50 33 41 31 33 42 42 33 32 40 33;
        34 44 45 34 27 35 34 25 25 34 36 26;
        43 35 45 36 35 25 26 35 26 27 35 34;
        43 44 36 26 34 36 27 27 36 35 25 36;
        37 47 48 37 30 38 37 28 28 37 39 29;
        46 38 48 39 38 28 29 38 29 30 38 37;
        46 47 39 29 37 39 30 30 39 38 28 39;
        40 50 51 40 33 41 40 31 31 40 42 32;
        49 41 51 42 41 31 32 41 32 33 41 40;
        49 50 42 32 40 42 33 33 42 41 31 42;
        43 52 52 43 45 44 43 43 43 43 45 44;
        52 44 52 45 44 43 44 44 44 45 44 43;
        52 52 45 44 43 45 45 45 45 44 43 45;
        46 52 52 46 48 47 46 46 46 46 48 47;
        52 47 52 48 47 46 47 47 47 48 47 46;
        52 52 48 47 46 48 48 48 48 47 46 48;
        49 52 52 49 51 50 49 49 49 49 51 50;
        52 50 52 51 50 49 50 50 50 51 50 49;
        52 52 51 50 49 51 51 51 51 50 49 51;
        52 52 52 52 52 52 52 52 52 52 52 52];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Each Euler angle is chosen so that the inverse rotation forces a specific element  %
    % of the rotation matrix to zero. zel(k,j,i) gives information about which element   %
    % ceases to be zero when a rotation around axis j is applied to pattern i.           %
    %    k=1 gives the index into a vectorized matrix of the entry that becomes non-zero %
    %    k=2 gives the index of the other element in the same column that changes        %
    %    k=3 gives the sign of the sine term affecting the first of these entries        %
    %    k=4 gives the sign of the initial value of the second of these entries if known %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    zel=reshape([  6  5  1  1  3  1 -1  1  2  1  1  1;
        2  3 -1  1  1  3  1  1  5  4  1  1;
        3  2  1  1  4  6  1  1  1  2 -1  1;
        5  6 -1 -1  3  1 -1 -1  2  1  1 -1;
        3  2  1 -1  6  4 -1 -1  1  2 -1 -1;
        2  3 -1 -1  1  3  1 -1  4  5 -1 -1;
        6  5  1 -1  3  1 -1  1  2  1  1  1;
        2  3 -1 -1  1  3  1 -1  5  4  1  1;
        3  2  1 -1  4  6  1 -1  1  2 -1 -1;
        5  6 -1  1  3  1 -1 -1  2  1  1 -1;
        3  2  1  1  6  4 -1 -1  1  2 -1  1;
        2  3 -1  1  1  3  1  1  4  5 -1  1;
        6  5  1  1  3  1 -1 -1  2  1  1 -1;
        2  3 -1 -1  1  3  1 -1  5  4  1 -1;
        3  2  1  1  4  6  1 -1  1  2 -1  1;
        5  6 -1  1  3  1 -1  1  2  1  1  1;
        3  2  1 -1  6  4 -1  1  1  2 -1 -1;
        2  3 -1  1  1  3  1  1  4  5 -1 -1;
        6  5  1 -1  3  1 -1 -1  2  1  1 -1;
        2  3 -1  1  1  3  1  1  5  4  1 -1;
        3  2  1 -1  4  6  1  1  1  2 -1 -1;
        5  6 -1 -1  3  1 -1  1  2  1  1  1;
        3  2  1  1  6  4 -1  1  1  2 -1  1;
        2  3 -1 -1  1  3  1 -1  4  5 -1  1;
        0  0  0  0  3  1 -1  1  2  1  1  1;
        3  2  1  1  0  0  0  0  1  2 -1  1;
        2  3 -1  1  1  3  1  1  0  0  0  0;
        0  0  0  0  6  4 -1  1  5  4  1  1;
        6  5  1  1  0  0  0  0  4  5 -1  1;
        5  6 -1  1  4  6  1  1  0  0  0  0;
        0  0  0  0  9  7 -1  1  8  7  1  1;
        9  8  1  1  0  0  0  0  7  8 -1  1;
        8  9 -1  1  7  9  1  1  0  0  0  0;
        0  0  0  0  3  1 -1 -1  2  1  1 -1;
        3  2  1 -1  0  0  0  0  1  2 -1 -1;
        2  3 -1 -1  1  3  1 -1  0  0  0  0;
        0  0  0  0  6  4 -1 -1  5  4  1 -1;
        6  5  1 -1  0  0  0  0  4  5 -1 -1;
        5  6 -1 -1  4  6  1 -1  0  0  0  0;
        0  0  0  0  9  7 -1 -1  8  7  1 -1;
        9  8  1 -1  0  0  0  0  7  8 -1 -1;
        8  9 -1 -1  7  9  1 -1  0  0  0  0;
        0  0  0  0  1  3  1  1  1  2 -1  1;
        2  3 -1  1  0  0  0  0  2  1  1  1;
        3  2  1  1  3  1 -1  1  0  0  0  0;
        0  0  0  0  4  6  1  1  4  5 -1  1;
        5  6 -1  1  0  0  0  0  5  4  1  1;
        6  5  1  1  6  4 -1  1  0  0  0  0;
        0  0  0  0  7  9  1  1  7  8 -1  1;
        8  9 -1  1  0  0  0  0  8  7  1  1;
        9  8  1  1  9  7 -1  1  0  0  0  0;
        0  0  0  0  0  0  0  0  0  0  0  0]',4,3,52);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert the m string
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~ischar(m) % lecacy call with integer m argument
    m=char(m+'w'); % convert to characters
end
ich=find(strcmp(m,mch),1);      % check if already in the cache
if isempty(ich)                 % not yet in the cache
    mm=m-'w';                   % convert to integers with x -> 1
    mi=mm>=-31 & mm<=-29;       % find characters XYZ
    mm(mi)=mm(mi)+32;           % convert XYZ to xyz (for compatibility)
    mi=mm>=-70 & mm<=-62;       % find digits 1:9
    mm(mi)=mm(mi)+74;           % convert to 4:12
    mi=mm<=0;                   % select control characters
    mc=mm(mi);                  % controls
    mm=mm(~mi);                 % rotations
    ef=1;                       % angle scale factor
    es=1;                       % angle sign
    fl=1;                       % default to no rotation matrix tranposing
    for i=1:length(mc)
        switch mc(i)
            case -5             % 'r' = radians
            case -19            % 'd' = degrees
                ef=pi/180;      % scale factor to convert to radians
            case -37            % 'R' = negated radians
                ef=-1;
            case -51            % 'D' = negated degrees
                ef=-pi/180;      % scale factor to convert to radians
            case -8             % 'o' = object-extrinsic
            case -40            % 'O' = object-intrinsic
                fl=-1;
                es=-1;
            case -22            % 'a' = axes-extrinsic
                fl=-1;
            case -54            % 'A' = axes-intrinsic
                es=-1;
            otherwise
                error('Invalid character: %s',mc(i)+'w')
        end
    end
    ef=ef*es;               % change sign of scale factor if necessary
    if es<0
        mm=mes(mm);         % sign-reverse: interchage 4,5,6 with 10,11,12
    end
    nm=length(mm);
    mv=zeros(7,nm+1);
    mv(1,:)=[mm 0];
    mv(2,:)=cumsum([mm<=3 0]);      % index into euler angle array
    mv(3,1)=1; % initial pattern is the identity matrix
    for i=1:nm                % loop for each rotation
        mmi=mm(i); % rotation code
        mv(3,i+1)=trmap(mv(3,i),mmi); % pattern ID after rotation
        if mmi<4
            mv(4:7,i)=zel(:,mmi,mv(3,i)); % information about which matrix elements change from zero
        end
    end
    mv(end)=fl;
    mv(end-3)=ef;
    % now save in the cache
    ich=jch(nch);       % find oldest cache entry
    mch{ich}=m;                 % save input string
    mvch{ich}=mv;               % save parameters
    jch=[ich jch(1:nch-1)];     % age all the other cache entries
else                            % already in the cache
    kch=find(jch==ich,1);       % find existing ich entry
    jch(1:kch)=[ich jch(1:kch-1)];
    mv=mvch{ich};               % retrieve from cache
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function voxel2FNpair = fgFNPairs(coords,node2voxel)
% Compute the (fiber,node) pairs in each of the coords
%
%   voxel2FNpair = fgFNPairs(coords,node2voxel)
%
% Explain me.  Part of the diffusion prediction code.
%
% See also: t_arcuatePlotXXX
%
% (c) Stanford VISTA Team

nFiber = length(node2voxel);

nCoords = size(coords,1);
voxel2FNpair = cell(nCoords,1);

for thisFiber=1:nFiber
    % We need to know both the nodes that pass through a voxel
    lst = (node2voxel{thisFiber}~=0);
    nodes = find(lst);

    % And we need to know which voxels they pass through
    voxelsInFiber = node2voxel{thisFiber}(lst);

    % Then we store for each voxel the (fiber,node) pairs that pass through
    % it.
    for jj=1:length(voxelsInFiber)
        voxel2FNpair{voxelsInFiber(jj)} = ...
            [voxel2FNpair{voxelsInFiber(jj)},[thisFiber,nodes(jj)]];
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [RSS, XYproj] = Residuals_ellipse(XY,ParG)
%
%   Projecting a given set of points onto an ellipse
%   and computing the distances from the points to the ellipse
%
%   This is a modified version of an iterative algorithm published by D. Eberly
%     Internet publication: "Distance from a point to an ellipse in 2D" (2004)
%                           Geometric Tools, LLC, www.geometrictools.com
%     Book publication: "3D Game Engine Design", 2nd edition.
%                       Morgan Kaufmann Publishers, San Francisco, CA, 2007.
%                              (see Section 14.13.1)
%
%   Input:  XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)
%           ParG is a vector 5x1 of the ellipse parameters
%           ParG =  [Center(1:2), Axes(1:2), Angle]'
%                Center - the coordinates of the ellipse's center
%                Axes   - the axes (major, minor)
%                Angle  - the angle of tilt of the ellipse
%
%   Output:  RSS is the Residual Sum of Squares (the sum of squares of the distances)
%            XYproj is the array of coordinates of projections
%
%   The algorithm is proven to converge and reaches an accuracy of 7-8 significant digit
%   It takes 4-5 iterations per point, on average.

Center = ParG(1:2);   Axes = ParG(3:4);  Angle = ParG(5);

n = size(XY,1);
XYproj = zeros(n,2);
tolerance = 1e-9;

%  First handling the circle case

if abs((Axes(1)-Axes(2))/Axes(1))<tolerance
    phiall = angle(XY(:,1)-Center(1) + sqrt(-1)*(XY(:,2)-Center(2)));
    XYproj = [Axes(1)*cos(phiall)+Center(1) Axes(2)*sin(phiall)+Center(2)];
    RSS = norm(XY-XYproj,'fro')^2;
    return;
end

%  Now dealing with proper ellipses

a = Axes(1); b = Axes(2);
aa = a^2;  bb = b^2;

tol_a = tolerance*a;
tol_b = tolerance*b;
tol_aa = tolerance*aa;

%  Matrix Q for rotating the points and the ellipse to the canonical system

s = sin(Angle); c = cos(Angle);
Q = [c -s; s c];

%  data points in canonical coordinates

XY0  = [XY(:,1)-Center(1) XY(:,2)-Center(2)]*Q;
XYA = abs(XY0);
Tini = max(a*(XYA(:,1)-a),b*(XYA(:,2)-b));

%  main loop over the data points

for i=1:n
    u = XYA(i,1);  v = XYA(i,2);
    ua = u*a;      vb = v*b;

    if (u == 0)
        z1 = 1;
    else
        z1 = sign(XY0(i,1));
    end

    if (v == 0)
        z2 = 1;
    else
        z2 = sign(XY0(i,2));
    end

    %       does the point lie on the minor axis?

    if u<tol_a
        if XY0(i,2)<0
            XYproj(i,:) = [0 -b];
        else
            XYproj(i,:) = [0 b];
        end
        continue;
    end

    %       does the point lie on the major axis?

    if v<tol_b
        if u<a-bb/a
            xproj = aa*u/(aa-bb);
            XYproj(i,:) = [z1*xproj z2*b*sqrt(1-(xproj/a)^2)];
        else
            XYproj(i,:) = [z1*a 0];
        end
        continue;
    end

    %      generic case: start the iterative procedure

    T = Tini(i);
    for iter=1:100
        Taa = T + aa;
        Tbb = T + bb;
        PP1 = (ua/Taa)^2;
        PP2 = (vb/Tbb)^2;
        F  = PP1 + PP2 - 1;
        if F<0; break; end;
        Fder = 2*(PP1/Taa + PP2/Tbb);
        Ratio = F/Fder;
        if (Ratio<tol_aa); break; end;
        T = T + Ratio;
    end

    %       compute the projection of the point onto the ellipse

    xproj = XY0(i,1)*aa/Taa;
    XYproj(i,:) = [xproj sign(XY0(i,2))*b*sqrt(1-(xproj/a)^2)];

end % end the main loop

%    rotate back to the original system
XYproj = XYproj*Q';
XYproj = [XYproj(:,1)+Center(1) XYproj(:,2)+Center(2)];
RSS = norm(XY-XYproj,'fro')^2;
end   % Residuals_ellipse

author: Sable

repo: mcbench-benchmarks


text:
function [basis,gradbasis]=legs(x,dir,n,scale)
% usage: [basis,gradbasis]=legs(x,dir,n,scale)
%
% returns the values and directional derivatives  of (n+1)^2-1 basis functions
% constructed from spherical harmonics at locations given in x and, for the
% gradients, for (in general non-normalized) directions given in dir.
%
% input: x      set of N locations given as an Nx3 matrix
%        dir    set of N direction vectors given as an Nx3 matrix
%                  (dir is not normalized (it hence can be a dipole moment))
%        n       order of spherical harmonics
%
% output: basis: Nx((n+1)^2-1)  matrix containing in the j.th  row the real
%                and imaginary parts of r^kY_{kl}(theta,Phi)/(N_{kl}*scale^k) ( (r,theta,phi)
%                are the spherical coordinates corresponding to  the j.th row in x)
%                for k=1 to n and l=0 to k
%                the order is:
%                          real parts for k=1 and l=0,1 (2 terms)
%                  then    imaginary parts for k=1 and l=1 (1 term)
%                  then    real parts for k=2 and l=0,1,2 (3 terms)
%                  then    imaginary parts for k=2 and l=1,2 (2 term)
%                              etc.
%                   the spherical harmonics are normalized with
%                   N_{kl}=sqrt(4pi (k+l)!/((k-l)!(2k+1)))
%                    the phase does not contain the usual (-1)^l term !!!
%                   scale is constant preferably set to the avererage radius
%
%         gradbasis: Nx((n+1)^2-1) matrix containing in the j.th row the scalar
%                     product of the gradient of the former with the j.th row of dir

% Copyright (C) 2003, Guido Nolte
%
% 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$

[n1,n2]=size(x);

comi=sqrt(-1);

normalize=ones(n,n+1);
facto=ones(2*n+2,1);for i=3:2*n+2,facto(i)=facto(i-1)*(i-1);end;
for i=1:n
    for j=1:i+1
        if j==1
            normalize(i,j)=scale^i*i*sqrt(2*facto(i+j-1+1)/facto(i-j+1+1)/(2*i+1));
        else
            normalize(i,j)=scale^i*i*sqrt(facto(i+j-1+1)/facto(i-j+1+1)/(2*i+1));
        end
    end
end
normalize=reshape(repmat(reshape(normalize,n*(n+1),1),1,n1),n,n+1,n1);

rad=sqrt(x(:,1).^2+x(:,2).^2+x(:,3).^2);
phi=angle(x(:,1)+comi*x(:,2));
costheta=x(:,3)./(rad+eps);

ms=0:n;
ns=1:n;

shiftfactors=zeros(n,n+1);
shiftminusfactors=zeros(n,n+1);
for in=1:n; for im=1:n+1
        shiftfactors(in,im)=in+im-1;
        shiftminusfactors(in,im)=(in+im-2)*(in+im-1);
    end;end;
for in=1:n
    shiftminusfactors(in,1)=1;
end

leg0=zeros(n,n+1,n1);

for i=1:n
 p=legendre(i,costheta);
 leg0(i,1:i+1,:)=p;
end
 leg0=leg0.*reshape(repmat((-1).^ms,n,n1),n,n+1,n1);
ephi=exp(comi*ms'*phi');

leg0=leg0.*reshape(repmat(reshape(ephi,1,(n+1)*n1),n,1),n,n+1,n1)...
         .*reshape(repmat((repmat(rad,1,n).^repmat(ns,n1,1))',n+1,1),n,n+1,n1);

     onesx=zeros(1,n+1,n1);
     onesx(1,1,:)=ones(1,1,n1);

legshift=[onesx;leg0(1:n-1,:,:)];
legshiftminus=[-conj(legshift(:,2,:)),legshift(:,1:n,:)].*reshape(repmat(reshape(shiftminusfactors,n*(n+1),1),n1,1),n,n+1,n1);
legshiftplus=-[legshift(:,2:n+1,:),zeros(n,1,n1)];
legshift=legshift.*reshape(repmat(reshape(shiftfactors,n*(n+1),1),n1,1),n,n+1,n1);

dirp=[(dir(:,1)+dir(:,2)/comi)/2,(dir(:,1)-dir(:,2)/comi)/2,dir(:,3)];

gradleg= reshape(repmat(transpose(dirp(:,1)),n*(n+1),1),n,n+1,n1).*legshiftplus...
        +reshape(repmat(transpose(dirp(:,2)),n*(n+1),1),n,n+1,n1).*legshiftminus...
        +reshape(repmat(transpose(dirp(:,3)),n*(n+1),1),n,n+1,n1).*legshift;

leg0=leg0./normalize;
gradleg=gradleg./normalize;

    basis=zeros(n1,(n+1)^2-1);gradbasis=zeros(n1,(n+1)^2-1);
 for i=1:n
     basis(:,i^2:i^2+i)=(reshape(real(leg0(i,1:i+1,:)),i+1,n1))';
     basis(:,i^2+i+1:(i+1)^2-1)=(reshape(imag(leg0(i,2:i+1,:)),i,n1))';
     gradbasis(:,i^2:i^2+i)=(reshape(real(gradleg(i,1:i+1,:)),i+1,n1))';
     gradbasis(:,i^2+i+1:(i+1)^2-1)=(reshape(imag(gradleg(i,2:i+1,:)),i,n1))';
 end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function qwv_2d_test01 ( )

%*****************************************************************************80
%
%% QWV_2D_TEST01 tests QWV_2D for a trio of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 May 2014
%
%  Author:
%
%    John Burkardt
%
  a = -1.0;
  b = +1.0;
  c = -1.0;
  d = +1.0;

  t = 1;
  n = ( ( t + 1 ) * ( t + 2 ) ) / 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QWV_2D_TEST01:\n' );
  fprintf ( 1, '  Compute the weights associated with an interpolatory\n' );
  fprintf ( 1, '  quadrature rule defined by N=(T+1)*(T+2)/2 points,\n' );
  fprintf ( 1, '  exact for polynomials of total degree T or less.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Degree T = %d\n', t );
  fprintf ( 1, '  Number of points N = %d\n', n );
  fprintf ( 1, '  X Interval = [%g,%g]\n', a, b );
  fprintf ( 1, '  Y Interval = [%g,%g]\n', c, d );
%
%  Set the points.
%
  x(1) =  1.0;
  y(1) = -1.0;
  x(2) =  1.0;
  y(2) =  1.0;
  x(3) = -1.0;
  y(3) =  1.0;

  r8vec2_print ( n, x, y, '  Abscissas:' )
%
%  Compute the weights.
%
  w = qwv_2d ( t, n, a, b, c, d, x, y );

  r8vec_print ( n, w, '  Weights:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x]=grandn(a,b)

% function [x]=grandn(a,b)
% random number generator from the Gamma (a,b) distribution

% part 1: obtain a random number from Gamma(a,1)

if a>=1
x=gammadrawover1(a);
else
xtilde=gammadrawover1(a+1);
u=rand;
x=xtilde*u^(1/a);
end

% part 2: transform into a draw from G(a,b)
x=x*b;

% auxiliary nested function to draw a Gamma random number when argument 'a' is greater than or equal to 1
function x=gammadrawover1(a)
d=a-1/3;
c=1/((9*d)^0.5);
check=0;
   while check==0
   z=randn;
   u=rand;
   v=(1+c*z)^3;
      if (v>0 && log(u)<0.5*z^2+d-d*v+d*log(v))
      x=d*v;
      check=1;
      else
      check=0;
      end
   end
end

% declare end of general function
end

author: european-central-bank

repo: BEAR-toolbox


text:
function en = end(a,k,n)
% Overloaded end function for file_array objects.
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

dim = size(a);
if k>length(dim)
    en = 1;
else
    if n<length(dim),
    dim = [dim(1:(n-1)) prod(dim(n:end))];
    end;
    en = dim(k);
end;


text:
% Local Regression and Likelihood, Figure 7.3.
%
% Censored local regression.
%
% Author: Catherine Loader
%
% NEEDS: Kaplan-Meier based estimate. Identify curves and cens/uncens. data

load heart;
fit = locfit(age,log(0.5+surv));
figure('Name','fig7_3: censored local regression Kaplan-Meier based;' );
lfplot(fit);
xlabel('Age at Transplant (years)');
ylabel('0.5+Survival Time (Days)');

fit = lf_censor(age,log(0.5+surv),cens);
hold on;
lfplot(fit,'nodata');
hold off;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function y = ismatrix(X)
    % True if X is a 2-Dimensional list

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % evaluate the operation in Mathematica and return the
    % expression string
    ret = eval_math(['Boole@MatrixQ[',X.s,']']);

    y = logical(eval(ret));
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function msfun_odometry(block)

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

%%
%% The setup method is used to set up the basic attributes of the
%% S-function such as ports, parameters, etc. Do not add any other
%% calls to the main body of the function.
%%
myPrevTicks = struct('right', 0, 'left', 0);

setup(block);

%endfunction

%% Function: setup ===================================================
%% Abstract:
%%   Set up the basic characteristics of the S-function block such as:
%%   - Input ports
%%   - Output ports
%%   - Dialog parameters
%%   - Options
%%
%%   Required         : Yes
%%   C-Mex counterpart: mdlInitializeSizes
%%
    function setup(block)

        % Register number of ports
        block.NumInputPorts  = 4;
        block.NumOutputPorts = 3;

        % Setup port properties to be inherited or dynamic
        % block.SetPreCompInpPortInfoToDynamic;
        % block.SetPreCompOutPortInfoToDynamic;
        block.SetPreCompPortInfoToDefaults;

        % Override input port properties
        for i=1:3
            block.InputPort(i).Dimensions        = 1;
            block.InputPort(i).DatatypeID  = 0;  % double
            block.InputPort(i).Complexity  = 'Real';
            block.InputPort(i).DirectFeedthrough = true;
        end

        block.InputPort(4).Dimensions        = 2;
        block.InputPort(4).DatatypeID  = 0;  % double
        block.InputPort(4).Complexity  = 'Real';
        block.InputPort(4).DirectFeedthrough = true;

        block.InputPort(1).SamplingMode = 0;
        block.InputPort(2).SamplingMode = 0;
        block.InputPort(3).SamplingMode = 0;
        block.InputPort(4).SamplingMode = 0;

        % Override output port properties
        block.OutputPort(1).Dimensions       = 1;
        block.OutputPort(1).DatatypeID  = 0; % double
        block.OutputPort(1).Complexity  = 'Real';

        block.OutputPort(2).Dimensions       = 1;
        block.OutputPort(2).DatatypeID  = 0; % double
        block.OutputPort(2).Complexity  = 'Real';

        block.OutputPort(2).Dimensions       = 1;
        block.OutputPort(2).DatatypeID  = 0; % double
        block.OutputPort(2).Complexity  = 'Real';

        % Set up the states
        block.NumContStates = 0;
        block.NumDworks = 0;

        % Register parameters
        block.NumDialogPrms     = 3;
        block.DialogPrmsTunable = {'NonTunable', 'NonTunable', 'NonTunable'};

        % Register sample times
        %  [0 offset]            : Continuous sample time
        %  [positive_num offset] : Discrete sample time
        %
        %  [-1, 0]               : Inherited sample time
        %  [-2, 0]               : Variable sample time
        block.SampleTimes = [0.033 0];

        % Specify the block simStateCompliance. The allowed values are:
        %    'UnknownSimState', < The default setting; warn and assume DefaultSimState
        %    'DefaultSimState', < Same sim state as a built-in block
        %    'HasNoSimState',   < No sim state
        %    'CustomSimState',  < Has GetSimState and SetSimState methods
        %    'DisallowSimState' < Error out when saving or restoring the model sim state
        block.SetAccelRunOnTLC(false); % run block in interpreted mode even w/ Acceleration
        block.SimStateCompliance = 'DefaultSimState';

        %% -----------------------------------------------------------------
        %% The MATLAB S-function uses an internal registry for all
        %% block methods. You should register all relevant methods
        %% (optional and required) as illustrated below. You may choose
        %% any suitable name for the methods and implement these methods
        %% as local functions within the same file. See comments
        %% provided for each function for more information.
        %% -----------------------------------------------------------------

        block.RegBlockMethod('PostPropagationSetup',    @DoPostPropSetup);
        block.RegBlockMethod('InitializeConditions', @InitializeConditions);
        block.RegBlockMethod('Start', @Start);
        block.RegBlockMethod('Outputs', @Outputs);     % Required
        block.RegBlockMethod('Update', @Update);
        block.RegBlockMethod('Terminate', @Terminate); % Required

    end %setup

%%
%% PostPropagationSetup:
%%   Functionality    : Setup work areas and state variables. Can
%%                      also register run-time methods here
%%   Required         : No
%%   C-Mex counterpart: mdlSetWorkWidths
%%
    function DoPostPropSetup(block)

    end

%%
%% InitializeConditions:
%%   Functionality    : Called at the start of simulation and if it is
%%                      present in an enabled subsystem configured to reset
%%                      states, it will be called when the enabled subsystem
%%                      restarts execution to reset the states.
%%   Required         : No
%%   C-MEX counterpart: mdlInitializeConditions
%%
    function InitializeConditions(block)

        myPrevTicks = struct('right', 0, 'left', 0);

    end

%%
%% Start:
%%   Functionality    : Called once at start of model execution. If you
%%                      have states that should be initialized once, this
%%                      is the place to do it.
%%   Required         : No
%%   C-MEX counterpart: mdlStart
%%
    function Start(block)

    end %endfunction

%%
%% Outputs:
%%   Functionality    : Called to generate block outputs in
%%                      simulation step
%%   Required         : Yes
%%   C-MEX counterpart: mdlOutputs
%%
    function Outputs(block)

        % Get wheel encoder ticks from the robot
        right_ticks = block.InputPort(4).Data(1);
        left_ticks = block.InputPort(4).Data(2);

        % Recal the previous wheel encoder ticks
        prev_right_ticks = myPrevTicks.right;
        prev_left_ticks = myPrevTicks.left;

        % Previous estimate
        x = block.InputPort(1).Data;
        y = block.InputPort(2).Data;
        theta = block.InputPort(3).Data;

        % Compute odometry here
        R = block.DialogPrm(1).Data;
        L = block.DialogPrm(2).Data;
        ticks_per_rev = block.DialogPrm(3).Data;
        m_per_tick = (2*pi*R)/ticks_per_rev;

        d_right = (right_ticks-prev_right_ticks)*m_per_tick;
        d_left = (left_ticks-prev_left_ticks)*m_per_tick;
        d_center = (d_right + d_left)/2;

        x_dt = d_center*cos(theta);
        y_dt = d_center*sin(theta);
        theta_dt = (d_right - d_left)/L;

        theta_new = theta + theta_dt;
        x_new = x + x_dt;
        y_new = y + y_dt;
        %             fprintf('Estimated pose (x,y,theta): (%0.3g,%0.3g,%0.3g)\n', x_new, y_new, theta_new);

        % Save the wheel encoder ticks for the next estimate
        myPrevTicks.right = right_ticks;
        myPrevTicks.left = left_ticks;

        % Update your estimate of (x,y,theta)
        block.OutputPort(1).Data = x_new;
        block.OutputPort(2).Data = y_new;
        block.OutputPort(3).Data = theta_new;

    end %end Outputs

%%
%% Update:
%%   Functionality    : Called to update discrete states
%%                      during simulation step
%%   Required         : No
%%   C-MEX counterpart: mdlUpdate
%%
    function Update(block)

    end %end Update

%%
%% Terminate:
%%   Functionality    : Called at the end of simulation for cleanup
%%   Required         : Yes
%%   C-MEX counterpart: mdlTerminate
%%
    function Terminate(block)

    end %end Terminate

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function simplex_model = find_simplex_models(p);

for i = 1:length(p)
    simplex_model(i)= 0;
    if p{i}.K.f == 0
        continue
    elseif any(p{i}.K.q > 0) | any(p{i}.K.s > 0)
        continue
    elseif p{i}.K.f ~= 1
        continue
    else
        aux = p{i};
        b = aux.F_struc(1,1);
        a = aux.F_struc(1,2:end);
        if all(abs(a) == 1)
            b = b/sign(-a(1));
            aux.F_struc(1:p{i}.K.f,:) = [];
            aux.K.f = 0;
            [aux,lower,upper] = find_simple_variable_bounds(aux);
            if all(lower == 0) & aux.K.l == 0 & all((upper == b) | isinf(upper))
                simplex_model(i)=1;
            end
        else
            continue
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef StiffnessTensor2VoigtConverterPS < FourthOrderTensor2VoigtConverterPS

    properties
    end

    methods (Access = public)

        function obj = StiffnessTensor2VoigtConverterPS(tensor)
           obj.computeConversion(tensor)
        end

    end

    methods (Access = protected)

        function selectVoigtTensorClass(obj)
            obj.voigtTensor = StiffnessPlaneStressVoigtTensor();
        end

    end

    methods (Access = protected,Static)

        function f = getVoigtFactor(iv,jv)
           f = 1;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function d = grdinfo2(file);
%GRDINFO2  Print information about a GMT grdfile (netCDF format, GMT v3 or v4)
%
% Uses NetCDF libraries to display information about a GMT grid file.
% Duplicates (some) functionality of the program grdinfo (which requires
% compilation as a mexcdf function on each architecture) using
% Matlab 2008b (and later) built-in NetCDF functionality
% instead of GMT libraries.
%
% GRDINFO2('file.grd') will display information about the GMT grid
% file 'file.grd' in a format similar to the gmt command grdinfo.
%
% D = GRDINFO('file.grd') will in addition return a vector containing
% (xmin, xmax, ymin, ymax, zmin, zmax, format, xinc, yinc). Format is
% 1 for pixel registration and 0 for grid node registration.
%
% See also GRDREAD2, GRDWRITE2

% This program is expected to work on any GMT netCDF format file,
% but it does not duplicate all the functionality of GMT's I/O
% library, so it will not work on all files supported by GMT.
% In particular, it will fail on binary format grdfiles.
% It is the responsibility of the user to determine whether this
% program is appropriate for any given task.
%
% For more information on GMT grid file formats, see:
% http://www.soest.hawaii.edu/gmt/gmt/doc/gmt/html/GMT_Docs/node70.html
% Details on Matlab's native netCDF capabilities are at:
% http://www.mathworks.com/access/helpdesk/help/techdoc/ref/netcdf.html

% GMT (Generic Mapping Tools, <http://gmt.soest.hawaii.edu>)
% was developed by Paul Wessel and Walter H. F. Smith

% Kelsey Jordahl
% Marymount Manhattan College
% http://marymount.mmm.edu/faculty/kjordahl/software.html

% Time-stamp: <Wed Jan  6 16:26:46 EST 2010>

% Version 1.1.1, 6-Jan-2010
% first released on MATLAB Central
% modification
% DB 	05/2014 	Suppress command window output

if nargin < 1,
  help(mfilename);
  return,
end

% check for appropriate Matlab version (>=7.7)
V=regexp(version,'[ \.]','split');
if (str2num(V{1})<7) | (str2num(V{1})==7 & str2num(V{2})<7),
  ver
  error('grdread2: Requires Matlab R2008b or later!');
end

ncid = netcdf.open(file, 'NC_NOWRITE');
if isempty(ncid),
  return,
end

[ndims,nvars,ngatts,unlimdimid] = netcdf.inq(ncid);

if (nvars==3),                        % new (v4) GMT netCDF grid file
%  for attnum=0:(ngatts-1),
  title=getatt_clean(ncid,'title');
  conv=getatt_clean(ncid,'Conventions');
  pixel=getatt_clean(ncid,'node_offset');
  if isempty(pixel), pixel=0; end
  desc=getatt_clean(ncid,'description');
  command=getatt_clean(ncid,'history');
  vers = getatt_clean(ncid,'GMT_version');
  xname = netcdf.getAtt(ncid,0,'long_name');
  yname = netcdf.getAtt(ncid,1,'long_name');
  zname = netcdf.getAtt(ncid,2,'long_name');
  xrange = netcdf.getAtt(ncid,0,'actual_range');
  yrange = netcdf.getAtt(ncid,1,'actual_range');
  zrange = netcdf.getAtt(ncid,2,'actual_range');
  [dimname,nx]=netcdf.inqDim(ncid,0);
  [dimname,ny]=netcdf.inqDim(ncid,1);
else
  if (nvars==6),                        % old (v3) GMT netCDF grid file
    [dimname, dimlen] = netcdf.inqDim(ncid,1);
    if (dimname=='xysize'),             % make sure it really is v3 netCDF
      title=getatt_clean(ncid,'title');
      command=getatt_clean(ncid,'source');
      conv=[];
      desc=[]; zrange=[0 0]; vers='3.x format file';
      xrange=netcdf.getVar(ncid,0)';
      xname=netcdf.getAtt(ncid,0,'units');
      yrange=netcdf.getVar(ncid,1)';
      yname=netcdf.getAtt(ncid,1,'units');
      zrange=netcdf.getVar(ncid,2)';
      zname=netcdf.getAtt(ncid,2,'units');
%      z=netcdf.getVar(ncid,5);
      dim=netcdf.getVar(ncid,4)';
      nx=dim(1); ny=dim(2);
      zname='z';
      pixel=netcdf.getAtt(ncid,5,'node_offset');
    else
      error('Apparently not a GMT netCDF grid');
    end
  else
    error('Wrong number of variables in netCDF file!');
  end
end

if pixel,                         % pixel node registered
  dx=diff(xrange)/double(nx); % convert int to double for division
  dy=diff(yrange)/double(ny);
else                              % gridline registered
  dx=diff(xrange)/double(nx-1); % convert int to double for division
  dy=diff(yrange)/double(ny-1);
end

%disp(['Title: ' title]);
%disp(['Conventions: ' conv]);
%disp(['GMT version: ' vers]);
%disp(['Command: ' command]);
%disp(['Remark: ' desc]);
%if pixel,
%  disp('Pixel node registration used');
%else
%  disp('Gridline node registration used');
%end
%disp(['x_min: ' num2str(xrange(1)) ' x_max: ' num2str(xrange(2)) ...
%                   ' name: ' xname])
%disp(['y_min: ' num2str(yrange(1)) ' y_max: ' num2str(yrange(2)) ...
%                   ' name: ' yname])
%disp(['z_min: ' num2str(zrange(1)) ' z_max: ' num2str(zrange(2)) ...
%                   ' name: ' zname])

switch nargout
  case 1,
   % need to make sure that zrange and pixel are cast into double precision
   d=[xrange(1) xrange(2) yrange(1) yrange(2) double(zrange(1)) double(zrange(2)) double(pixel) dx dy];
end

function val = getatt_clean(ncid,attname)
%
% Call netcdf.getAtt when not sure whether global attribute 'attname' exists.
% Trap error and return empty vector in that case.
%
eval(['val = netcdf.getAtt(ncid,netcdf.getConstant(''NC_GLOBAL''),''' attname ''');'],'val=[];');

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
% phd_test.m
% ====================================================>
% This is a test script which demonstrates the usage of the "SMC_PHD" class.
%
% SETUP:
%  * Before running the simulation, open "2_crossing_targets.mat" or "3_roaming_targets.mat" datasets, from the "datasets" folder
%  * The datasets have been extracted by simulating the motion of differential robots in a 2D-plane (x,y)
%  * The "gen_obs_cluttered_multi3" function takes as an input the ground truth data, including information about the measurement noise and clutter rate
%     and then produces 1xNk cell array of corrupted and cluttered measurements, Nk being the total number of timesteps

% Load dataset
load('multiple-robot-tracking.mat');

tot_ellapsed = 0;
% Plot settings
ShowPlots = 1;              % Set to 0 to hide plots
ShowPrediction = 0;         % Set to 0 to skip showing prediction
ShowUpdate = 1;             % Set to 0 to skip showing update
SmoothTrajectories = 0;

% Recording settings
clear F;
Record = 1;                 % Set to (0|1) to turn video recording (off|on)
FrameRate = 10;            % Number of frames per second
VideoQuality = 100;         % Set to desired quality percentage
VideoPathName = 'tomb_tracks_only.avi'; % Set to the desired path and name of produced recording

lambdaV = 1; % Expected number of clutter measurements over entire surveillance region
V = 10^2;     % Volume of surveillance region (10x10 2D-grid)
V_bounds = [0 10 0 10]; % [x_min x_max y_min y_max]

% Instantiate a Transitionamic model
transition_model = ConstantVelocityX('NumDims',2,'VelocityErrVariance',0.0001);

% Instantiate a Measurement model
%measurement_model = LinearGaussianX('NumMeasDims',2,'NumStateDims',4,'MeasurementErrVariance',0.02,'Mapping',[1 3]);
measurement_model = RangeBearing2CartesianX('NumStateDims',4,'MeasurementErrVariance',[(pi/50)^2,0.02],'Mapping',[1 3]);

% Instantiate a clutter model
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,'Limits',[V_bounds(1:2);V_bounds(3:4)]);

% Instantiate birth model
numBirthComponents = 10;
BirthComponents.Means = [ 1 9 9 1; 0 0 0 0; 1 1 9 9; 0 0 0 0];
BirthComponents.Covars = repmat(diag([2,0.1,2,0.1]),1,1,4);
BirthComponents.Weights = [.25, .25, .25, .25];
birth_distribution = GaussianMixtureX(BirthComponents.Means,BirthComponents.Covars, BirthComponents.Weights);
birth_model = DistributionBasedBirthModelX('Distribution', birth_distribution,...
                                           'BirthIntensity', 0.000001);

% Compile the State-Space model
ssm = StateSpaceModelX(transition_model,measurement_model,'Clutter',clutter_model, 'Birth', birth_model);

% Extract the ground truth data from the example workspace
load('example.mat');
NumIter = size(GroundTruth,2);

% Set BirthIntensity
NumTracks = 3;

% Generate DataList
meas_simulator = MultiTargetMeasurementSimulatorX('Model',ssm);
%meas_simulator.DetectionProbability = 1;
DataList = meas_simulator.simulate(GroundTruthStateSequence);

% Assign PHD parameter values
config.Model = ssm;
config.SurvivalProbability = 0.9;
config.DetectionProbability = 0.9;

% Instantiate PHD filter
filter = TrackOrientedMeMBerPoissonGMFilterX(config);
filter.Poisson.StatePosterior = copy(birth_distribution);

% Create figure windows
if(ShowPlots)
    img = imread('maze.png');

    % set the range of the axes
    % The image will be stretched to this.
    min_x = 0;
    max_x = 10;
    min_y = 0;
    max_y = 10;

    % make data to plot - just a line.
    x = min_x:max_x;
    y = (6/8)*x;

    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
    axis(ax(1), 'manual');
    figure('units','normalized','outerposition',[.5 0 .5 1])
    ax(2) = gca;

    axis
end

% START OF SIMULATION
% ===================>
for k=1:NumIter
    fprintf('Iteration = %d/%d\n================>\n',k,NumIter);

%     % Extract DataList at time k
%     tempDataList = DataList{k}(:,:);
%     tempDataList( :, ~any(tempDataList,1) ) = [];

    % Change PHD filter parameters
    filter.MeasurementList = DataList(k); % New observations

    tic;
    % Predict PHD filter
    filter.predict();

    % Update PHD filter
    filter.update();
    ellapsed = toc;
    tot_ellapsed = tot_ellapsed + ellapsed;
    %fprintf("Probability of existence: %f\n", filter.ProbOfExistence);
    fprintf("Ellapsed time: %f\n\n", ellapsed);
    % Plot update step results
    if(ShowPlots && ShowUpdate)
        % Plot data
        cla(ax(1));
         % Flip the image upside down before showing it

        % NOTE: if your image is RGB, you should use flipdim(img, 1) instead of flipud.
        hold on;
        axis(ax(1),V_bounds)
        hold(ax(1),'on');

        hold(ax(1),'on');
        for i=1:filter.Bernoulli.StatePosterior.NumComponents
            if(filter.Bernoulli.StatePosterior.Weights(i)>0.8)
                means = filter.Bernoulli.StatePosterior.Trajectories{i}.StateMean;
                if SmoothTrajectories
                    means = smoothdata(means','gaussian',10)';
                end
                plot(ax(1),means(1,:),means(3,:),'.-');
                plot_gaussian_ellipsoid(filter.Bernoulli.StatePosterior.Trajectories{i}.StateMean([1,3],end),filter.Bernoulli.StatePosterior.Trajectories{i}.StateCovar([1,3],[1,3],end),'r',1,50,ax(1));
            end
        end

        h2 = plot(ax(1), DataList(k).Vectors(1,:),DataList(k).Vectors(2,:),'k*','MarkerSize', 10);
        str = sprintf('Robot positions (Update)');
        title(ax(1),str)
        xlabel('X position (m)')
        ylabel('Y position (m)')

        % Plot PHD
        cla(ax(2), 'reset');
        p = filter.Bernoulli.StatePosterior.random(100000);
        [bandwidth,density,X,Y]=kde2d(p([1,3],:)');
        %contour3(X,Y,density,50);
        h = surf(ax(2),X,Y,density);
        shading interp
        colormap(ax(2), jet(3000))
        %set(h, 'edgecolor','none')
        hold on;
        plot(ax(2), filter.MeasurementList.Vectors(1,:), filter.MeasurementList.Vectors(2,:), 'y*');
        axis(ax(2), [V_bounds]);
        str = sprintf('PHD intensity (Update)');
        xlabel(ax(2),'X position (m)')
        ylabel(ax(2),'Y position (m)')
        zlabel(ax(2),'Intensity')
        title(ax(2),str)
        pause(0.01)

        % Store video frame
        if(Record)
            F(k) = getframe(ax(1));
        end
    end
end

% Create video file and write to it
if(Record)
    %F = F(2:end);
    vidObj = VideoWriter(char(VideoPathName));
    vidObj.Quality = VideoQuality;
    vidObj.FrameRate = FrameRate;
    open(vidObj);
    writeVideo(vidObj, F);
    close(vidObj);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function s = char(ker)

% CHAR
%
% Return a textual representation of a polynomial kernel object.
%
%    str = char(ker);

%
% File        : @polynomial/char.m
%
% Date        : Tuesday 12th Spetember 2000
%
% Author      : Dr Gavin C. Cawley
%
% Description : Cast a polynomial kernel object to a character string.
%               Part of an object-oriented implementation of Vapnik's Support
%               Vector Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 07/07/2000 - v1.00
%               12/09/2000 - v1.01 minor improvements to comments and help
%                                  messages
%
% Copyright   : (c) Dr Gavin C. Cawley, September 2000.
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

s = ['polynomial (d = ' num2str(ker.d) ')'];

% bye bye...