text:
function [b] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf)
%function [b] = compute_rf_spsp_mgh(kp,rfp,gz,kz,kf)
%Function that computes the complex-valued SPSP RF pulse waveform,
%using parameters specified in structure rfp.
%Inputs:
%kp: k-space trajectory parameter structure, generated by kparameterSPSP.m
%rfp:RF waveform design parameter structure, generated by rfparameterSPSP.m
%gz: z gradient waveform vector
%kz,kf: SPSP k-space (kf-kz space) trajectory vectors
%Output:
%b:complex-value SPSP RF pulse waveform in units of Gauss (g).
%
%
%Chun-yu Yip, 4/1/2009

%Physical parameters
gam = 26751;                       %rad/s/g; gyromagnetic ratio
gambar = gam/2/pi;                 %Hz/g

%Create grid in SPSP space
dresz = rfp.dfovz/rfp.ddimz;       %cm
dresf = rfp.dfovf/rfp.ddimf;       %Hz
z = [-rfp.dfovz/2:dresz:rfp.dfovz/2-dresz];
f = [-rfp.dfovf/2:dresf:rfp.dfovf/2-dresf];
[Z,F] = ndgrid(z,f);

disp('Creating desired pattern for SPSP pulse...');

%phase pattern
TD = rfp.TE-kp.pointtime*kp.npnts/2; %pulse end to acquisition of DC sample
dphs = gam * rfp.alpha * TD * F .* (Z-rfp.zshift); %Eq. 5

%magnitude pattern
if streq(rfp.profileshape,'rect')

    dmag = zeros(size(Z));
    dmag(find((Z<(rfp.slthickz/2)) & (Z>(-rfp.slthickz/2)))) = 1;

elseif streq(rfp.profileshape,'gaussian')

    a = -4*log(2)/rfp.slthickz/rfp.slthickz;
    dmag = exp(a*(Z.^2));

else

    disp('Slice profile shape unknown.');

end

if rfp.smoothprofile

    kernel = gaussian(z,0,rfp.kernelstdv)';
    dmag=convn_fft(dmag,kernel);     %Smoothing implemented in Fourier domain

end;

dmag = dmag/max(max(dmag));
dmag = sin(rfp.flipangle/180 * pi) * dmag;
d = exp(1i*dphs).*dmag;

%Define matrices and vectors for conjugate gradient
b0 = zeros(length(kz),1);            %Initiation of CG
mask = ones(size(d));                %ROI in f-z space
W = spdiag(mask);                    %weighting based on ROI definition

%total energy penalty
beta = rfp.beta;
C = beta*speye(length(kz));

M = length(Z(:));
precon =speye(length(kz));           %CG Preconditioner
Niter = rfp.Niter;                   %Number of CG iterations

%Define matlab object "A" to represent the gigantic system matrix.
%"A" is of type "Fatrix" from Jeff Fessler's toolbox.
disp('Forming spatial-spectral fast object...');
A=spsp_Af(z,f,kz,kf,kp.pointtime,rfp.dfovf,rfp.dfovz);

%%%%%%%%%%%Conjugate Gradient%%%%%%%%%%%%%%%%%%%
disp('Running preconditioned conjugate gradient...');

%keyboard
[b_all, info] = qpwls_pcg(b0, A, W, d(:), zeros(M,1), sqrt(C), precon, Niter);
b = b_all(:,Niter);

if rfp.hamming
scale = 2;          %This parameter controls how wide the window is.
hamm=(0.54-0.46.*cos(pi*((kz/scale/max(kz))+1)));
b= b.*hamm;         %Intended to smooth the excited slice profile.

end;

m = A*b;            %Small-tip-angle approximation of excitation pattern
mm = reshape(m,size(d));

disp('Displaying SPSP pulse design results...');

t = kp.pointtime*[0:1:length(gz)-1].';

figure;
subplot(2,1,1);
plot(t*1000,gz);
grid;
x1=xlabel('Time (ms)');
y1=ylabel('g/cm');
t1=title('z gradient');
set(t1,'Fontsize',14);
set(x1,'Fontsize',12);
set(y1,'Fontsize',12);

subplot(2,1,2);
hold on;
plot(t*1000,real(b),'b');
plot(t*1000,imag(b),'r');
grid;
x2=xlabel('Time (ms)');
y2=ylabel('g');
t2=title('SPSP pulse');
l2=legend('Real','Imaginary');
set(t2,'Fontsize',14);
set(x2,'Fontsize',12);
set(y2,'Fontsize',12);
set(l2,'Fontsize',12);

%Display desired pattern of the SPSP recovery pulse.
figure;
subplot(2,1,1);
imagesc(f,z,abs(d));
colormap default;colorbar;
x1 = xlabel('Frequency offset (Hz)');
y1 = ylabel('z (cm)');
t1 = title('Desired SPSP pattern (magnitude)');

set(x1,'Fontsize',12);
set(y1,'Fontsize',12);
set(t1,'Fontsize',12);

subplot(2,1,2);
imagesc(f,z,angle(d).*abs(d),[-pi,pi]);colorbar;
x2 =xlabel('Frequency offset (Hz)');
y2=ylabel('z (cm)');
t2=title('Desired SPSP pattern (phase)');
set(x2,'Fontsize',12);
set(y2,'Fontsize',12);
set(t2,'Fontsize',12);

figure;
subplot(2,1,1);
imagesc(f,z,abs(mm));
xlabel('frequency (Hz)');
ylabel('z (cm)');
title('Resulting pattern (linear regime prediction)');

subplot(2,1,2);
imagesc(f,z,angle(mm));
xlabel('frequency (Hz)');
ylabel('z (cm)');
title('Resulting phase pattern (linear regime prediction)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end


text:
function [] = checkMlpCost()

% Check the gradients for the MLP
%

addpath '../library/'

%% Setup random data / small model

visibleDim = 9;  % width/height of input
visibleSize = visibleDim^2;

hiddenDimL1 = 4;  % width/height of hidden layer 1
hiddenSizeL1 = hiddenDimL1^2;
hiddenViewDimL1 = 3;  % width/height of local receptive area for hidden layer 1
hiddenViewSizeL1 = hiddenViewDimL1^2;  % local receptive area size for hidden layer 1
hiddenViewStepL1 = 2;  % step for receptive area moving across input
assert (((hiddenDimL1-1) * hiddenViewStepL1 + hiddenViewDimL1) == visibleDim);

hiddenDimL2 = 2;  % width/height of hidden layer 2
hiddenSizeL2 = hiddenDimL2^2;

hiddenViewDimL2 = 3;  % width/height of local receptive area for hidden layer 2
hiddenViewSizeL2 = hiddenViewDimL2^2;  % local receptive area size for hidden layer 2
hiddenViewStepL2 = 1;  % step for receptive area moving across hidden layer 1
assert (((hiddenDimL2-1) * hiddenViewStepL2 + hiddenViewDimL2) == hiddenDimL1);

[W1, b1] = initializeOneLayerParams(hiddenSizeL1, hiddenViewSizeL1);
[W2, b2] = initializeOneLayerParams(hiddenSizeL2, hiddenViewSizeL2);
[W3, b3] = initializeOneLayerParams(1, hiddenSizeL2);

W1Indices = buildIndices(visibleDim, hiddenDimL1, hiddenViewDimL1, hiddenViewStepL1);
W2Indices = buildIndices(hiddenDimL1, hiddenDimL2, hiddenViewDimL2, hiddenViewStepL2);

assert(all(size(W1Indices) == [hiddenViewSizeL1, hiddenSizeL1]));
assert(all(size(W2Indices) == [hiddenViewSizeL2, hiddenSizeL2]));

theta = [W1(:);b1(:);W2(:);b2(:);W3(:);b3(:)];

lambda = 0.01;
data   = randn(visibleSize, 7);
labels = [ 0 1 0 1 0 1 0];

[cost, grad] = mlpCost(theta, visibleSize, ...
                       hiddenSizeL1, hiddenViewSizeL1, ...
                       hiddenSizeL2, hiddenViewSizeL2, ...
                       W1Indices, W2Indices, ...
                       lambda, data, labels);

% Check that the numerical and analytic gradients are the same
numgrad = computeNumericalGradient( @(x) mlpCost(x, visibleSize, ...
                                                 hiddenSizeL1, hiddenViewSizeL1, ...
                                                 hiddenSizeL2, hiddenViewSizeL2, ...
                                                 W1Indices, W2Indices, ...
                                                 lambda, data, labels), ...
                                   theta);

% Use this to visually compare the gradients side by side
% disp([numgrad grad (numgrad-grad)]);

% Compare numerically computed gradients with the ones obtained from backpropagation
% disp('Norm between numerical and analytical gradient (should be less than 1e-9)');
diff = norm(numgrad-grad)/norm(numgrad+grad);

% disp(diff); % Should be small. In our implementation, these values are
              % usually less than 1e-9.
assert (diff < 1e-9);

end

author: zellyn


text:
function B = mirror_expand(A)
    [row,col,channel] = size(A);
    yi = 2:row+1;
    xi = 2:col+1;
    B = zeros(row+2,col+2,channel);
    B(yi,xi,:) = A;
    B([1 row+2],[1 col+2],:) = B([3 row],[3 col],:);  % mirror corners
    B([1 row+2],xi,:) = B([3 row],xi,:);          % mirror left and right boundary
    B(yi,[1 col+2],:) = B(yi,[3 col],:);          % mirror top and bottom boundary
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [mask3M] = getDown3Mask(mask3M, sampleTrans, sampleAxis)
%Downsample a 3D matrix: but return a mask the same size as the original,
%showing where the samples are taken.
%JOD. 16 Nov 03.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

sV  = size(mask3M);

sampleSlices = ceil(sV(3));

mask3M = repmat(logical(0), sV);

indV = 1 : sampleAxis: sampleSlices;
for i = 1 : length(indV)
  maskM = getDown2Mask(mask3M(:,:,indV(i)),sampleTrans);
  mask3M(:,:,indV(i)) = maskM;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, 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 dyn(r)
%SerialLink.dyn Display inertial properties
%
% R.dyn() displays the inertial properties of the SerialLink object in a multi-line
% format.  The properties shown are mass, centre of mass, inertia, gear ratio,
% motor inertia and motor friction.
%
% See also Link.dyn.
	for j=1:r.n
		fprintf('----- link %d\n', j);
		l = r.links(j);
		l.dyn()
	end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ ALL_HYPS ] = globalSampling( hyps_obj, number, typenum )
%GLOBALSAMPLING Summary of this function goes here
%   Detailed explanation goes here
% room_xyz = [min(room_obj.out_points_w,[],1) max(room_obj.out_points_w,[],1)];
% typenum = 12;
global config

%% rotate pairwise constraints
load(config.pairwiseName);
angleset = [0 -pi/2 pi +pi/2];
rotate_pairwise = cell(size(pairwise,1), size(pairwise,2), 4);
rotate_pairwise_norm = cell(size(pairwise_norm,1), size(pairwise_norm,2), 4);
rotate_pairwise_wall = cell(size(pairwise_wall,1), size(pairwise_wall,2), 4);

for a = 1:4
    for i = 1:size(pairwise,1)
        for j = 1:size(pairwise,2)
            if ~isempty(pairwise{i,j})
                p = pairwise{i,j};
                r = p;
                r(:,1) = p(:,1).*cos(angleset(a)) + p(:,2).*sin(angleset(a));
                r(:,2) = -p(:,1).*sin(angleset(a)) + p(:,2).*cos(angleset(a));
                rotate_pairwise{i,j,a} = r;

                p = pairwise_norm{i,j};
                r = p;
                r(:,1) = p(:,1).*cos(angleset(a)) + p(:,2).*sin(angleset(a));
                r(:,2) = -p(:,1).*sin(angleset(a)) + p(:,2).*cos(angleset(a));
                rotate_pairwise_norm{i,j,a} = r;

                p = pairwise_wall{i,j};
                r = p;
                r(:,1) = p(:,1).*cos(angleset(a)) + p(:,2).*sin(angleset(a));
                r(:,2) = -p(:,1).*sin(angleset(a)) + p(:,2).*cos(angleset(a));
                rotate_pairwise_wall{i,j,a} = r;
            end
        end
    end
end

%% convert every score to 0~1 by sigmoid function
OBJFEA = hyps_obj.objfea;
unaryscore = vertcat(OBJFEA.unaryscore);
rectscore = vertcat(OBJFEA.rectscore);
segmscore = vertcat(OBJFEA.segmscore);
randomforest = vertcat(OBJFEA.randomforest);

unaryscore = sigmoidFunc_flex(unaryscore, 33.67, 0.15);
rectscore = sigmoidFunc_flex(rectscore, -4, 0);
segmscore = ones(length(segmscore),1);
% segmscore = sigmoidFunc_flex(segmscore, -4, 0.5);
randomforest = sigmoidFunc_flex(randomforest, -4, 0.5);

dpmvalid = true(length(rectscore),typenum);
% dpmscore = vertcat(OBJFEA.dpmscore);
% dpmvalid(:,10) = dpmscore(:,1)>-0.9;
% dpmvalid(:,12) = dpmscore(:,2)>-0.95;

imgvalid = rectscore.*segmscore>0.1;
semvalid = unaryscore(:,1:typenum).*randomforest(:,1:typenum)>0.3;

totalvalid = dpmvalid & semvalid & repmat(imgvalid, 1, typenum);
totalid = cell(typenum,1);
for i = 1:typenum
    totalid{i} = find(totalvalid(:,i));
end

% rectscore = sqrt(rectscore);
% scr10 =  sigmoidFunc_flex(dpmscore(:,1), -50, -0.9);
% scr12 =  sigmoidFunc_flex(dpmscore(:,2), -15, -0.9);
% unaryscore(:,10) = unaryscore(:,10).*scr10;
% unaryscore(:,12) = unaryscore(:,12).*scr12;

%% decide object sequence and number
numberbias = exp(0.02*(-1:-1:-10));
load(config.numberName);
ALL_TYPERANK = cell(number,1);
for setid = 1:number
    sample_score = zeros(1,0);
    sample_ids = zeros(1,0);
    num_of_type = zeros(1,typenum);
    for oid = 1:typenum
        num_of_type(oid) = randp( prob(oid,1:10).*numberbias, 1) - 1;
        likelihood = rectscore.*segmscore.*unaryscore(:,oid).*randomforest(:,oid);
        sample_score = [sample_score likelihood(randp(likelihood, 1, num_of_type(oid)))'];
        sample_ids = [sample_ids oid*ones(1, num_of_type(oid))];
    end
    [~, I] = sort(sample_score, 'descend');
    ALL_TYPERANK{setid} = sample_ids(I);
end

%% start sampling
occ_thres = 0.05;
par_thres = [50 0.3 50];
ori_sel_hyps = repmat(struct('fixed',false,'selID',[]), typenum, 1);

ALL_HYPS = repmat(struct('sel_hyps',[]), number, 1);
for count = 1:number % parfor
    fprintf('Whole Room Hypothesis: %d/%d\n', count, number);
    typerank = ALL_TYPERANK{count};
    sel_hyps = ori_sel_hyps;

    degree = 2*ones(1,length(typerank));

    for samid = 1:length(typerank)
        tid = typerank(samid);
        unary_semantic = unaryscore(totalvalid(:,tid),tid).*randomforest(totalvalid(:,tid),tid); % type_hyps(tid).unary_score; %type_hyps(tid).h_scores(:, tid); %
        unary_imageevi = rectscore(totalvalid(:,tid)).*segmscore(totalvalid(:,tid)); % type_hyps(tid).image_score; %ones(length(unary_semantic),1); %
        pairwise_abso = ones(length(unary_semantic),1);
        pairwise_norm = ones(length(unary_semantic),1);
        pairwise_wall = ones(length(unary_semantic),1);
        pairwise_occl = ones(length(unary_semantic),1);

        % check fixed objects and remove invalid hypothesis
        xyz = hyps_obj.obj_xyz(totalvalid(:,tid),:);
        for i = 1:typenum
            if sel_hyps(i).fixed
                for j = sel_hyps(i).selID
                    if ~((tid==7 && i==10) || (i==7 && tid==10))
                        o = occlusionTest(hyps_obj.obj_xyz(j,:), xyz);
                    else
                        o = zeros(size(xyz,1),1);
                    end
%                     locdata = data;
                    locdata.seedangle = hyps_obj.anglesid(j);
                    locdata.testangle = hyps_obj.anglesid(totalvalid(:,tid));
                    p1 = pairwiseHyperTest( hyps_obj.obj_xyz(j,:), xyz, ...
                        rotate_pairwise{i,tid,hyps_obj.anglesid(j)}, ...
                        1, i, tid, locdata);
%                     p2 = pairwiseHyperTest( type_hyps(i).h_obj_xyz(j,:), xyz, ...
%                         rotate_pairwise_norm{i,tid,type_hyps(i).anglesid(j)}, ...
%                         2, i, tid, data);
                    p3 = pairwiseHyperTest( hyps_obj.obj_xyz(j,:), xyz, ...
                        rotate_pairwise_wall{i,tid,hyps_obj.anglesid(j)}, ...
                        3, i, tid, locdata);

                    pairwise_occl = pairwise_occl .* (o<occ_thres);
                    pairwise_abso = pairwise_abso .* sigmoidFunc(p1, par_thres(1), 5);
%                     pairwise_norm = pairwise_norm .* sigmoidFunc(p2, 0.3, 5);
                    pairwise_wall = pairwise_wall .* sigmoidFunc(p3, par_thres(3), 5);
%                     validID = validID & o<occ_thres & p1<par_thres(1) & p2<par_thres(2) & p3<par_thres(3);
                end
            end
        end

        % randomly pick one
        sub_prob = unary_semantic.*unary_imageevi.*pairwise_abso.*pairwise_wall.*pairwise_occl;
        if sum(sub_prob.^degree(samid))==0
            continue;
        end

        try
            sub_sel = totalid{tid}(randp(sub_prob.^degree(samid), 1));
        catch
            continue;
        end

        sel_hyps(tid).fixed = true;
        sel_hyps(tid).selID = [sel_hyps(tid).selID sub_sel];

    end

    ALL_HYPS(count).sel_hyps = sel_hyps;
end

end

%%
function p = sigmoidFunc(v, m, a)
p = 1./(1+exp(a/m*(v-m)));
end

function p = sigmoidFunc_flex(v, a, b)
% p = 1./(1+exp(a*v+b));
p = 1./(1+exp(a*(v-b)));
end

author: zouchuhang

repo: LayoutNet


text:
function varargout = selStream(varargin)
% SELSTREAM MATLAB code for selStream.fig
%      SELSTREAM, by itself, creates a new SELSTREAM or raises the existing
%      singleton*.
%
%      H = SELSTREAM returns the handle to a new SELSTREAM or the handle to
%      the existing singleton*.
%
%      SELSTREAM('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SELSTREAM.M with the given input arguments.
%
%      SELSTREAM('Property','Value',...) creates a new SELSTREAM or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before selStream_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to selStream_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help selStream

% Last Modified by GUIDE v2.5 28-Jan-2015 13:19:58

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @selStream_OpeningFcn, ...
                   'gui_OutputFcn',  @selStream_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before selStream is made visible.
function selStream_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to selStream (see VARARGIN)

% Choose default command line output for selStream
handles.output = hObject;
if length(varargin) == 1
    handles.streamnames = varargin{1};
else
    errordlg('Input correct streamnames sturcture.');
end

nStream = length(handles.streamnames);
menuOptions = [];
for it=1:nStream
   menuOptions = [menuOptions;handles.streamnames(it)];
end
set(handles.popupmenu1,'String',menuOptions);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes selStream wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = selStream_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
% varargout{1} = handles.output;

% --- Executes on selection change in popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1

% items = get(hObject,'String');
% idx_sel = get(hObject,'Value');
% item_sel = items(idx_sel);
%
% assignin('base','streamname',item_sel);

% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

items = get(handles.popupmenu1,'String');
idx_sel = get(handles.popupmenu1,'Value');
item_sel = items(idx_sel);
assignin('base','streamname',item_sel);

% Hint: delete(hObject) closes the figure
if isequal(get(hObject, 'waitstatus'), 'waiting')
    % The GUI is still in UIWAIT, us UIRESUME
    uiresume(hObject);
%     guidata(hObject,handles);
end
delete(hObject);

repo: REST


function [ occludingVertices ] = occludingBoundaryVertices( FV,Ef,Ev,R )
%OCCLUDINGBOUNDARYVERTICES Return vertices in FV that lie on boundary edges
%   FV - mesh structure containing vertices and faces
%   Ef - faces adjacent to each edge in the mesh
%   Ev - vertices adjacent to each face in the mesh
%
% This runs during edge fitting. As shape/pose changes, so do the vertices
% that lie on the occluding boundary. So these probably want to be updated
% every iteration.

if nargin>3
    FV.vertices = (R*FV.vertices')';
end

% Compute face normals
fn = facenormals( FV.vertices,FV.faces );

% Edges with a zero index lie on the mesh boundary, i.e. they are only
% adjacent to one face
boundaryEdges = Ef(:,1)==0;

% To avoid matlab indexing problems, replace zero index with an arbitrary
% value (these edges are going to be ignored anyway)
idx1 = Ef(:,1);
idx1(boundaryEdges)=1;

% Compute the occluding edges as those where the two adjacent face normals
% differ in the sign of their Z component
occludingEdges = (sign(fn(idx1,3)) ~= sign(fn(Ef(:,2),3))) & ~boundaryEdges;

% Select the vertices lying at the end of the occluding edges and remove
% duplicates
occludingVertices = Ev(occludingEdges,:);
occludingVertices = unique(occludingVertices(:));

% Remove vertices from occluding boundary list that are not visible
occludingVertices = intersect(occludingVertices,visiblevertices(FV, eye(3)));

end


text:
function db=lpccc2db(cc,np,nc)
%LPCCC2DBf Convert complex cepstrum to dB power spectrum DB=(CC,NP,NC)
% cc = complex cepstrum coefs
% np+2 = number of frequency values 0 to nyquist (default: np=size(cc,2))
% nc = number of cepstral coefs to use (default: nc=np);
% for high speed make np one less than a power of 2

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpccc2db.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,mc]=size(cc);
if nargin<2 np=mc+2; end
if nargin<3 nc=np; end
k=10/log(10);
if nc==mc
   db=k*rsfft([zeros(nf,1) cc].',2*np+2).';
else
   db=k*rsfft([zeros(nf,1) lpccc2cc(cc,nc)].',2*np+2).';
end
if nargout==0
   plot((0:np+1)/(2*np+2),db.');
   ylabel('dB');
end

author: decouples


text:
function sD = som_normalize(sD,method,comps)

%SOM_NORMALIZE (Re)normalize data or add new normalizations.
%
% sS = som_normalize(sS,[method],[comps])
%
%   sS = som_normalize(sD)
%   sS = som_normalize(sS,sNorm)
%    D = som_normalize(D,'var')
%   sS = som_normalize(sS,'histC',[1:3 10])
%
%  Input and output arguments ([]'s are optional):
%   sS                The data to which the normalization is applied.
%                     The modified and updated data is returned.
%            (struct) data or map struct
%            (matrix) data matrix (a matrix is also returned)
%   [method]          The normalization method(s) to add/use. If missing,
%                     or an empty variable ('') is given, the
%                     normalizations in sS are used.
%            (string) identifier for a normalization method to be added:
%                     'var', 'range', 'log', 'logistic', 'histD' or 'histC'.
%            (struct) Normalization struct, or an array of such.
%                     Alternatively, a map/data struct can be given
%                     in which case its '.comp_norm' field is used
%                     (see below).
%            (cell array) Of normalization structs. Typically, the
%                     '.comp_norm' field of a map/data struct. The
%                     length of the array must be equal to data dimension.
%            (cellstr array) norm and denorm operations in a cellstr array
%                     which are evaluated with EVAL command with variable
%                     name 'x' reserved for the variable.
%   [comps]  (vector) the components to which the normalization is
%                     applied, default is [1:dim] ie. all components
%
% For more help, try 'type som_normalize' or check out online documentation.
% See also SOM_DENORMALIZE, SOM_NORM_VARIABLE, SOM_INFO.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_normalize
%
% PURPOSE
%
% Add/apply/redo normalization on data structs/sets.
%
% SYNTAX
%
%  sS = som_normalize(sS)
%  sS = som_normalize(sS,method)
%   D = som_normalize(D,sNorm)
%  sS = som_normalize(sS,csNorm)
%  sS = som_normalize(...,comps)
%
% DESCRIPTION
%
% This function is used to (initialize and) add, redo and apply
% normalizations on data/map structs/sets. If a data/map struct is given,
% the specified normalizations are added to the '.comp_norm' field of the
% struct after ensuring that all normalizations specified therein have
% status 'done'. SOM_NORMALIZE actually uses function SOM_NORM_VARIABLE
% to handle the normalization operations, and only handles the data
% struct/set specific stuff itself.
%
% The different normalization methods are listed below. For more
% detailed descriptions, see SOM_NORM_VARIABLE.
%
%   method     description
%   'var'      Variance is normalized to one (linear operation).
%   'range'    Values are normalized between [0,1] (linear operation).
%   'log'      Natural logarithm is applied to the values:
%                xnew = log(x-m+1)
%              where m = min(x).
%   'logistic' Logistic or softmax trasformation which scales all
%              possible values between [0,1].
%   'histD'    Histogram equalization, values scaled between [0,1].
%   'histC'    Approximate histogram equalization with partially
%              linear operations. Values scaled between [0,1].
%   'eval'     freeform operations
%
% To enable undoing and applying the exactly same normalization to
% other data sets, normalization information is saved into a
% normalization struct, which has the fields:
%
%   .type   ; struct type, ='som_norm'
%   .method ; normalization method, a string
%   .params ; normalization parameters
%   .status ; string: 'uninit', 'undone' or 'done'
%
% Normalizations are always one-variable operations. In the data and map
% structs the normalization information for each component is saved in the
% '.comp_norm' field, which is a cell array of length dim. Each cell
% contains normalizations for one vector component in a struct array of
% normalization structs. Each component may have different amounts of
% different kinds of normalizations. Typically, all normalizations are
% either 'undone' or 'done', but in special situations this may not be the
% case. The easiest way to check out the status of the normalizations is to
% use function SOM_INFO, e.g. som_info(sS,3)
%
% REQUIRED INPUT ARGUMENTS
%
%   sS                The data to which the normalization is applied.
%            (struct) Data or map struct. Before adding any new
%                     normalizations, it is ensured that the
%                     normalizations for the specified components in the
%                     '.comp_norm' field have status 'done'.
%            (matrix) data matrix
%
% OPTIONAL INPUT ARGUMENTS
%
%   method            The normalization(s) to add/use. If missing,
%                     or an empty variable ('' or []) is given, the
%                     normalizations in the data struct are used.
%            (string) Identifier for a normalization method to be added:
%                     'var', 'range', 'log', 'logistic', 'histD' or 'histC'. The
%                     same method is applied to all specified components
%                     (given in comps). The normalizations are first
%                     initialized (for each component separately, of
%                     course) and then applied.
%            (struct) Normalization struct, or an array of structs, which
%                     is applied to all specified components. If the
%                     '.status' field of the struct(s) is 'uninit',
%                     the normalization(s) is initialized first.
%                     Alternatively, the struct may be map or data struct
%                     in which case its '.comp_norm' field is used
%                     (see the cell array option below).
%            (cell array) In practice, the '.comp_norm' field of
%                     a data/map struct. The length of the array
%                     must be equal to the dimension of the given
%                     data set (sS). Each cell contains the
%                     normalization(s) for one component. Only the
%                     normalizations listed in comps argument are
%                     applied though.
%            (cellstr array) norm and denorm operations in a cellstr array
%                     which are evaluated with EVAL command with variable
%                     name 'x' reserved for the variable.
%
%   comps    (vector) The components to which the normalization(s) is
%                     applied. Default is to apply to all components.
%
% OUTPUT ARGUMENTS
%
%   sS                Modified and/or updated data.
%            (struct) If a struct was given as input argument, the
%                     same struct is returned with normalized data and
%                     updated '.comp_norm' fields.
%            (matrix) If a matrix was given as input argument, the
%                     normalized data matrix is returned.
%
% EXAMPLES
%
%  To add (initialize and apply) a normalization to a data struct:
%
%    sS = som_normalize(sS,'var');
%
%  This uses 'var'-method to all components. To add a method only to
%  a few selected components, use the comps argument:
%
%    sS = som_normalize(sS,'log',[1 3:5]);
%
%  To ensure that all normalization operations have indeed been done:
%
%    sS = som_normalize(sS);
%
%  The same for only a few components:
%
%    sS = som_normalize(sS,'',[1 3:5]);
%
%  To apply the normalizations of a data struct sS to a new data set D:
%
%    D = som_normalize(D,sS);
%  or
%    D = som_normalize(D,sS.comp_norm);
%
%  To normalize a data set:
%
%    D = som_normalize(D,'histD');
%
%  Note that in this case the normalization information is lost.
%
%  To check out the status of normalization in a struct use SOM_INFO:
%
%    som_info(sS,3)
%
%
% SEE ALSO
%
%  som_denormalize    Undo normalizations of a data struct/set.
%  som_norm_variable  Normalization operations for a set of scalar values.
%  som_info           User-friendly information of SOM Toolbox structs.

% Copyright (c) 1998-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 151199 150500

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

error(nargchk(1, 3, nargin));  % check no. of input arguments is correct

% sD
struct_mode = isstruct(sD);
if struct_mode,
  switch sD.type
   case 'som_map', D = sD.codebook;
   case 'som_data', D = sD.data;
   otherwise, error('Illegal struct.')
  end
else
  D = sD;
end
[dlen dim] = size(D);

% comps
if nargin<3 || (ischar(comps) && strcmp(comps,'all')),
  comps = [1:dim];
end
if isempty(comps), return; end
if size(comps,1)>1, comps = comps'; end  % make it a row vector

% method
csNorm = cell(dim,1);
if nargin<2 || isempty(method),
  if ~struct_mode,
    warning('No normalization method given. Data left unchanged.');
    return;
  end
  method = '';
else
  % check out the given method
  % (and if necessary, copy it for each specified component)
  if ischar(method),
    switch method,
    case {'var','range','log','histD','histC','logistic'},
      sN = som_set('som_norm','method',method);
    otherwise,
      error(['Unrecognized method: ' method]);
    end
    for i=comps, csNorm{i} = sN; end
  elseif isstruct(method),
    switch method(1).type,
    case {'som_map','som_data'}, csNorm = method(1).comp_norm;
    case {'som_norm'}, for i=comps, csNorm{i} = method; end
    otherwise,
      error('Invalid struct given as normalization method.')
    end
  elseif iscellstr(method),
    [dummy,sN] = som_norm_variable(1,method,'init');
    for i=comps, csNorm{i} = sN; end
  elseif iscell(method),
    csNorm = method;
  else
    error('Illegal method argument.')
  end
  % check the size of csNorm is the same as data dimension
  if length(csNorm) ~= dim,
    error('Given number of normalizations does not match data dimension.')
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize

% make sure all the current normalizations for current
% components have been done
if struct_mode,
  alldone = 1;
  for i = comps,
    for j=1:length(sD.comp_norm{i}),
      sN = sD.comp_norm{i}(j);
      if ~strcmp(sN.status,'done'),
	alldone = 0;
        [x,sN] = som_norm_variable(D(:,i), sN, 'do');
        D(:,i) = x;
        sD.comp_norm{i}(j) = sN;
      end
    end
  end
  if isempty(method),
    if alldone,
      warning('No ''undone'' normalizations found. Data left unchanged.');
    else
      fprintf(1,'Normalizations have been redone.\n');
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% action

% add the new normalizations to the old ones
for i = comps,
  if ~isempty(csNorm{i}),
    [x,sN] = som_norm_variable(D(:,i), csNorm{i}, 'do');
    D(:,i) = x;
    if struct_mode,
      if isempty(sD.comp_norm{i}), sD.comp_norm{i} = sN;
      else sD.comp_norm{i} = [sD.comp_norm{i}, sN]; end
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% output

if struct_mode,
  switch sD.type
   case 'som_map', sD.codebook = D;
   case 'som_data', sD.data = D;
   otherwise, error('Illegal struct.')
  end
else
  sD = D;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: ilarinieminen

repo: SOM-Toolbox


text:
function [tpr,tnr,info] = vl_roc(labels, scores, varargin)
%VL_ROC   ROC curve.
%   [TPR,TNR] = VL_ROC(LABELS, SCORES) computes the Receiver Operating
%   Characteristic (ROC) curve [1]. LABELS is a row vector of ground
%   truth labels, greater than zero for a positive sample and smaller
%   than zero for a negative one. SCORES is a row vector of
%   corresponding sample scores, usually obtained from a
%   classifier. The scores induce a ranking of the samples where
%   larger scores should correspond to positive labels.
%
%   Without output arguments, the function plots the ROC graph of the
%   specified data in the current graphical axis.
%
%   Otherwise, the function returns the true positive and true
%   negative rates TPR and TNR. These are vectors of the same size of
%   LABELS and SCORES and are computed as follows. Samples are ranked
%   by decreasing scores, starting from rank 1. TPR(K) and TNR(K) are
%   the true positive and true negative rates when samples of rank
%   smaller or equal to K-1 are predicted to be positive. So for
%   example TPR(3) is the true positive rate when the two samples with
%   largest score are predicted to be positive. Similarly, TPR(1) is
%   the true positive rate when no samples are predicted to be
%   positive, i.e. the constant 0.
%
%   Setting a label to zero ignores the corresponding sample in the
%   calculations, as if the sample was removed from the data. Setting
%   the score of a sample to -INF causes the function to assume that
%   that sample was never retrieved. If there are samples with -INF
%   score, the ROC curve is incomplete as the maximum recall is less
%   than 1.
%
%   [TPR,TNR,INFO] = VL_ROC(...) returns an additional structure INFO
%   with the following fields:
%
%   info.auc:: Area under the ROC curve (AUC).
%     This is the area under the ROC plot, the parametric curve
%     (FPR(S), TPR(S)). The PLOT option can be used to plot variants
%     of this curve, which affects the calculation of a corresponding
%     AUC.
%
%   info.eer:: Equal error rate (EER).
%     The equal error rate is the value of FPR (or FNR) when the ROC
%     curves intersects the line connecting (0,0) to (1,1).
%
%   info.eerThreshold:: EER threshold.
%     The value of the score for which the EER is attained.
%
%   VL_ROC() accepts the following options:
%
%   Plot:: []
%     Setting this option turns on plotting unconditionally. The
%     following plot variants are supported:
%
%     tntp:: Plot TPR against TNR (standard ROC plot).
%     tptn:: Plot TNR against TPR (recall on the horizontal axis).
%     fptp:: Plot TPR against FPR.
%     fpfn:: Plot FNR against FPR (similar to a DET curve).
%
%     Note that this option will affect the INFO.AUC value computation
%     too.
%
%   NumPositives:: []
%   NumNegatives:: []
%     If either of these parameters is set to a number, the function
%     pretends that LABELS contains the specified number of
%     positive/negative labels. NUMPOSITIVES/NUMNEGATIVES cannot be
%     smaller than the actual number of positive/negative entries in
%     LABELS. The additional positive/negative labels are appended to
%     the end of the sequence as if they had -INF scores (as explained
%     above, the function interprets such samples as `not
%     retrieved'). This feature can be used to evaluate the
%     performance of a large-scale retrieval experiment in which only
%     a subset of highly-scoring results are recorded for efficiency
%     reason.
%
%   Stable:: false
%     If set to true, TPR and TNR are returned in the same order
%     of LABELS and SCORES rather than being sorted by decreasing
%     score.
%
%   About the ROC curve::
%     Consider a classifier that predicts as positive all samples whose
%     score is not smaller than a threshold S. The ROC curve represents
%     the performance of such classifier as the threshold S is
%     changed. Formally, define
%
%       P = overall num. of positive samples,
%       N = overall num. of negative samples,
%
%     and for each threshold S
%
%       TP(S) = num. of samples that are correctly classified as positive,
%       TN(S) = num. of samples that are correctly classified as negative,
%       FP(S) = num. of samples that are incorrectly classified as positive,
%       FN(S) = num. of samples that are incorrectly classified as negative.
%
%     Consider also the rates:
%
%       TPR = TP(S) / P,      FNR = FN(S) / P,
%       TNR = TN(S) / N,      FPR = FP(S) / N,
%
%     and notice that, by definition,
%
%       P = TP(S) + FN(S) ,    N = TN(S) + FP(S),
%       1 = TPR(S) + FNR(S),   1 = TNR(S) + FPR(S).
%
%     The ROC curve is the parametric curve (FPR(S), TPR(S)) obtained
%     as the classifier threshold S is varied in the reals. The TPR is
%     the same as `recall' in a PR curve (see VL_PR()).
%
%     The ROC curve is contained in the square with vertices (0,0) The
%     (average) ROC curve of a random classifier is a line which
%     connects (1,0) and (0,1).
%
%     The ROC curve is independent of the prior probability of the
%     labels (i.e. of P/(P+N) and N/(P+N)).
%
%   REFERENCES:
%   [1] http://en.wikipedia.org/wiki/Receiver_operating_characteristic
%
%   See also: VL_PR(), VL_DET(), VL_HELP().

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

[tp, fp, p, n, perm, varargin] = vl_tpfp(labels, scores, varargin{:}) ;
opts.plot = [] ;
opts.stable = false ;
opts = vl_argparse(opts,varargin) ;

% compute the rates
small = 1e-10 ;
tpr = tp / max(p, small) ;
fpr = fp / max(n, small) ;
fnr = 1 - tpr ;
tnr = 1 - fpr ;

do_plots = ~isempty(opts.plot) || nargout == 0 ;
if isempty(opts.plot), opts.plot = 'fptp' ; end

% --------------------------------------------------------------------
%                                                      Additional info
% --------------------------------------------------------------------

if nargout > 2 || do_plots
  % Area under the curve. Since the curve is a staircase (in the
  % sense that for each sample either tn is decremented by one
  % or tp is incremented by one but the other remains fixed),
  % the integral is particularly simple and exact.

  switch opts.plot
    case 'tntp', info.auc = -sum(tpr .* diff([0 tnr])) ;
    case 'fptp', info.auc = +sum(tpr .* diff([0 fpr])) ;
    case 'tptn', info.auc = +sum(tnr .* diff([0 tpr])) ;
    case 'fpfn', info.auc = +sum(fnr .* diff([0 fpr])) ;
    otherwise
      error('''%s'' is not a valid PLOT type.', opts.plot);
  end

  % Equal error rate. One must find the index S in correspondence of
  % which TNR(S) and TPR(s) cross. Note that TPR(S) is non-decreasing,
  % TNR(S) is non-increasing, and from rank S to rank S+1 only one of
  % the two quantities can change. Hence there are exactly two types
  % of crossing points:
  %
  %  1) TNR(S) = TNR(S+1) = EER and TPR(S) <= EER, TPR(S+1) > EER,
  %  2) TPR(S) = TPR(S+1) = EER and TNR(S) > EER, TNR(S+1) <= EER.
  %
  % Moreover, if the maximum TPR is smaller than 1, then it is
  % possible that neither of the two cases realizes. In the latter
  % case, we return EER=NaN.

  s = max(find(tnr > tpr)) ;
  if s == length(tpr)
    info.eer = NaN ;
    info.eerThreshold = 0 ;
  else
    if tpr(s) == tpr(s+1)
      info.eer = 1 - tpr(s) ;
    else
      info.eer = 1 - tnr(s) ;
    end
    info.eerThreshold = scores(perm(s)) ;
  end
end

% --------------------------------------------------------------------
%                                                                 Plot
% --------------------------------------------------------------------

if do_plots
  cla ; hold on ;
  switch lower(opts.plot)
    case 'tntp'
      hroc = plot(tnr, tpr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [1 0], 'r--', 'linewidth', 2) ;
      spline([0 1], [0 1], 'k--', 'linewidth', 1) ;
      plot(1-info.eer, 1-info.eer, 'k*', 'linewidth', 1) ;
      xlabel('true negative rate') ;
      ylabel('true positive rate (recall)') ;
      loc = 'sw' ;

    case 'fptp'
      hroc = plot(fpr, tpr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [0 1], 'r--', 'linewidth', 2) ;
      spline([1 0], [0 1], 'k--', 'linewidth', 1) ;
      plot(info.eer, 1-info.eer, 'k*', 'linewidth', 1) ;
      xlabel('false positive rate') ;
      ylabel('true positive rate (recall)') ;
      loc = 'se' ;

    case 'tptn'
      hroc = plot(tpr, tnr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [1 0], 'r--', 'linewidth', 2) ;
      spline([0 1], [0 1], 'k--', 'linewidth', 1) ;
      plot(1-info.eer, 1-info.eer, 'k*', 'linewidth', 1) ;
      xlabel('true positive rate (recall)') ;
      ylabel('false positive rate') ;
      loc = 'sw' ;

    case 'fpfn'
      hroc = plot(fpr, fnr, 'b', 'linewidth', 2) ;
      hrand = spline([0 1], [1 0], 'r--', 'linewidth', 2) ;
      spline([0 1], [0 1], 'k--', 'linewidth', 1) ;
      plot(info.eer, info.eer, 'k*', 'linewidth', 1) ;
      xlabel('false positive (false alarm) rate') ;
      ylabel('false negative (miss) rate') ;
      loc = 'ne' ;

    otherwise
      error('''%s'' is not a valid PLOT type.', opts.plot);
  end

  grid on ;
  xlim([0 1]) ;
  ylim([0 1]) ;
  axis square ;
  title(sprintf('ROC (AUC: %.2f%%, EER: %.2f%%)', info.auc * 100, info.eer * 100), ...
        'interpreter', 'none') ;
  legend([hroc hrand], 'ROC', 'ROC rand.', 'location', loc) ;
end

% --------------------------------------------------------------------
%                                                        Stable output
% --------------------------------------------------------------------

if opts.stable
  tpr(1) = [] ;
  tnr(1) = [] ;
  tpr_ = tpr ;
  tnr_ = tnr ;
  tpr = NaN(size(tpr)) ;
  tnr = NaN(size(tnr)) ;
  tpr(perm) = tpr_ ;
  tnr(perm) = tnr_ ;
end

% --------------------------------------------------------------------
function h = spline(x,y,spec,varargin)
% --------------------------------------------------------------------
prop = vl_linespec2prop(spec) ;
h = line(x,y,prop{:},varargin{:}) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function PVal=computeCoherenceSignificance(TH,N)
% PVal=computeCoherenceSignificance(TH,N)
%
% PURPOSE:
% Returns a rough estimate of the significance value of a
% coherence threshold. 'Rough' because it assumes
% gaussian noise and takes no account of the fact that voxels
% are spatially correlated with their neighbours
% and have non-delta temporal auto-correlation functions.
%
% From Bantettini el al, 1993,
% Processing Strategies for Time Course Data Sets
% MRM 30:161-173 (1993) pp 171
%
% ARW 01/14/03 Wrote it
% $Author: wade $
% $Date: 2003/01/14 23:07:47 $

if (~exist('TH','var') | ~exist('N','var'))
    error('You must pass both a coherence threshold "TH" and the number of TRs "N"');
end

% It's just the complimentary error function...
PVal=1-erf(TH*sqrt(N/2))+eps;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function z = r82_to_c8 ( x, y )

%*****************************************************************************80
%
%% R82_TO_C8 makes a complex number from two real values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the real and imaginary parts.
%
%    Output, complex Z, the value of ( A + B*i).
%
  z = x + y * i;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_bug2148

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_connectivitysimulation ft_freqanalysis ft_connectivityanalysis ft_connectivityplot ft_freqdescriptives ft_checkdata

% ft_checkdata (ie fixdimord) always wants to reverts dimord to 'chan'
% although 'labelcmb' is present. The only way to make e.g.
% ft_connectivitplot work atm is to use:
% ft_checkdata(data, 'cmbstyle', 'full');
% This will make dimord chan_chan_freq and convert labelcmb to label.

cfg             = [];
cfg.ntrials     = 500;
cfg.triallength = 1;
cfg.fsample     = 200;
cfg.nsignal     = 4;
cfg.method      = 'ar';
cfg.params(:,:,1) = [.8   0   0 .2;
  0  .9  .5  0;
  .4   0  .5 .3;
  0  .2   0 .7];
cfg.params(:,:,2) = [-.5   0   0  0;
  0 -.8   0  0;
  0   0 -.2  0;
  0   0   0 .1];
cfg.noisecov      = [.3 0  0  0;
  0 1  0  0;
  0 0 .2  0;
  0 0  0  .4];

data            = ft_connectivitysimulation(cfg);
data.label      = {'AFz', 'Cz', 'Fz', 'Oz'};

% freqanalysis
cfgf           = [];
cfgf.method    = 'mtmfft';
cfgf.output    = 'fourier';
cfgf.tapsmofrq = 2;
cfgf.keeptrials= 'yes';
freq           = ft_freqanalysis(cfgf, data);
ft_freqdescriptives([], freq);

% connectivityanalysis
cfgc           = [];
cfgc.channelcmb = {'all' freq.label{1}};
cfgc.method    = 'coh';
c1             = ft_connectivityanalysis(cfgc, freq);
ft_connectivityplot([], c1);

% plotting
cfgp = [];
cfgp.layout    = 'EEG1010.lay';
cfgp.directionality = 'inflow';
cfgp.parameter = 'cohspctrm';
cfgp.refchannel = 'AFz';
ft_topoplotTFR(cfgp, c1);

% do JM's juicy transforms
tmp = ft_checkdata(freq, 'cmbstyle', 'fullfast');
c1             = ft_connectivityanalysis(cfgc, freq);
try
  ft_freqdescriptives([], tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end

ft_topoplotTFR(cfgp, c1);
ft_connectivityplot([], c1);

cfga = [];
cfga.parameter = 'fourierspctrm';
try
  ft_appendfreq(cfga, freq, tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end
vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.unit = 'cm';
vol.cond = 0.5;
vol.type = 'singlesphere';
cfgs.headmodel = vol;
cfgs.frequency = 10;
cfgs.elecfile = 'standard_1020.elc';
ft_sourceanalysis(cfgs, tmp);

tmp = ft_checkdata(tmp, 'cmbstyle', 'sparse');
c1             = ft_connectivityanalysis(cfgc, freq);
try
  ft_freqdescriptives([], tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end
ft_topoplotTFR(cfgp, c1);
ft_connectivityplot([], c1);
try
  ft_appendfreq(cfga, freq, tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end

vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.cond = 0.5;
vol.unit = 'cm';
vol.type = 'singlesphere';
cfgs.headmodel = vol;
cfgs.frequency = 10;
cfgs.elecfile = 'standard_1020.elc';
ft_sourceanalysis(cfgs, tmp);

tmp = ft_checkdata(tmp, 'cmbstyle', 'sparsewithpow');
c1             = ft_connectivityanalysis(cfgc, freq);
ft_freqdescriptives([], tmp);
ft_topoplotTFR(cfgp, c1);
ft_connectivityplot([], c1);
try
  ft_appendfreq(cfga, freq, tmp);
catch err
  fprintf('error: %s @ %s @ line %i\n ', err.message, err.stack(1).name, err.stack(1).line);
end

vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.cond = 0.5;
vol.unit = 'cm';
vol.type = 'singlesphere';
cfgs.headmodel = vol;
cfgs.elecfile = 'standard_1020.elc';
cfgs.frequency = 10;
ft_sourceanalysis(cfgs, tmp);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ stats ] = ml_get_cv_grid_states(test,train)
%ML_GET_CV_GRID_STATS Compute statistics on K-fold Cross Validation (CV).
%
%   input -----------------------------------------------------------------
%
%       o test  : cell array (M x 1), M is the number of parameters for
%                                     which K-fold CV has been evaluated on.
%
%           test{i}.accuracy  : (1 x K), K number of
%           test{i}.precision : (P x K), P is the total number of classes
%           test{i}.recall    : (P x K)
%           test{i}.fmeasure  : (P x K)
%
%
%       o train : cell array (M x 1),
%
%               train{i}: structure is the same as for test
%
%   output ----------------------------------------------------------------
%
%       o stats : struct , returns statistics on train and test sets
%
%           stats.train  : (statistics on the training data set)
%
%               - stats.train.acc : (statistics on the accuracy for the training set)
%
%                   - stats.train.acc.mean
%                   - stats.train.acc.std
%
%       the structure is the same for the test statistics
%
%% Check if grid search was done on > 1 parameters

[P, N] = size(train);

if N > 1
    test = test(:);
    train = train(:);
end

%% Extract Statitistics
M = length(train);
K = length(train{1}.accuracy);

% For Training Data
stats.train.acc.mean = zeros(1,K);
stats.train.acc.std  = zeros(1,K);

stats.train.fmeasure.mean = zeros(1,K);
stats.train.fmeasure.std  = zeros(1,K);

stats.train.fpr.mean = zeros(1,K);
stats.train.fpr.std  = zeros(1,K);

stats.train.tnr.mean = zeros(1,K);
stats.train.tnr.std  = zeros(1,K);

% bla bli
% Model Statistics for SVM
stats.model.totSV.mean   = zeros(1,K);
stats.model.totSV.std    = zeros(1,K);

stats.model.ratioSV.mean = zeros(1,K);
stats.model.ratioSV.std  = zeros(1,K);

stats.model.posSV.mean   = zeros(1,K);
stats.model.posSV.std    = zeros(1,K);

stats.model.negSV.mean   = zeros(1,K);
stats.model.negSV.std    = zeros(1,K);

stats.model.boundSV.mean = zeros(1,K);
stats.model.boundSV.std  = zeros(1,K);

% For Testing Data
if ~isempty(test)
    stats.test.acc.mean = zeros(1,K);
    stats.test.acc.std  = zeros(1,K);

    stats.test.fmeasure.mean = zeros(1,K);
    stats.test.fmeasure.std  = zeros(1,K);

    stats.test.fpr.mean = zeros(1,K);
    stats.test.fpr.std  = zeros(1,K);

    stats.test.tnr.mean = zeros(1,K);
    stats.test.tnr.std  = zeros(1,K);
end

for m=1:M

    % For Testing Data
    if ~isempty(test)
        stats.test.acc.mean(m)      = mean(test{m}.accuracy);
        stats.test.acc.std(m)       = std(test{m}.accuracy);

        stats.test.fmeasure.mean(m) = mean(test{m}.fmeasure);
        stats.test.fmeasure.std(m)  = std(test{m}.fmeasure);

        stats.test.fpr.mean(m)      = mean(test{m}.fpr);
        stats.test.fpr.std(m)       = std(test{m}.fpr);

        stats.test.tnr.mean(m)      = mean(test{m}.tnr);
        stats.test.tnr.std(m)       = std(test{m}.tnr);
    end

    % For Training Data
    stats.train.acc.mean(m)      = mean(train{m}.accuracy);
    stats.train.acc.std(m)       = std(train{m}.accuracy);

    stats.train.fmeasure.mean(m) = mean(train{m}.fmeasure);
    stats.train.fmeasure.std(m)  = std(train{m}.fmeasure);

    stats.train.fpr.mean(m)      = mean(train{m}.fpr);
    stats.train.fpr.std(m)       = std(train{m}.fpr);

    stats.train.tnr.mean(m)      = mean(train{m}.tnr);
    stats.train.tnr.std(m)       = std(train{m}.tnr);

    % Model Statistics for SVM
    stats.model.totSV.mean(m)   = mean(train{m}.totSV);
    stats.model.totSV.std(m)    = std(train{m}.totSV);

    stats.model.ratioSV.mean(m) = mean(train{m}.ratioSV);
    stats.model.ratioSV.std(m)  = std(train{m}.ratioSV);

    stats.model.posSV.mean(m)   = mean(train{m}.posSV);
    stats.model.posSV.std(m)    = std(train{m}.posSV);

    stats.model.negSV.mean(m)   = mean(train{m}.negSV);
    stats.model.negSV.std(m)    = std(train{m}.negSV);

    stats.model.boundSV.mean(m) = mean(train{m}.boundSV);
    stats.model.boundSV.std(m)  = std(train{m}.boundSV);

end

%% Reshape stats if grid search on > 1 parameters
if N > 1
    % For Testing Data
    if ~isempty(test)
        stats.test.acc.mean = reshape(stats.test.acc.mean, P, N);
        stats.test.acc.std  = reshape(stats.test.acc.std, P, N);

        stats.test.fmeasure.mean = reshape(stats.test.fmeasure.mean, P, N);
        stats.test.fmeasure.std  = reshape(stats.test.fmeasure.std, P, N);

        stats.test.fpr.mean = reshape(stats.test.fpr.mean, P, N);
        stats.test.fpr.std  = reshape(stats.test.fpr.std, P, N);

        stats.test.tnr.mean = reshape(stats.test.tnr.mean, P, N);
        stats.test.tnr.std  = reshape(stats.test.tnr.std, P, N);

    end

    % For Training Data
    stats.train.acc.mean = reshape(stats.train.acc.mean, P, N);
    stats.train.acc.std  = reshape(stats.train.acc.std, P, N);

    stats.train.fmeasure.mean = reshape(stats.train.fmeasure.mean, P, N);
    stats.train.fmeasure.std  = reshape(stats.train.fmeasure.std, P, N);

    stats.train.fpr.mean = reshape(stats.train.fpr.mean, P, N);
    stats.train.fpr.std  = reshape(stats.train.fpr.std, P, N);

    stats.train.tnr.mean = reshape(stats.train.tnr.mean, P, N);
    stats.train.tnr.std  = reshape(stats.train.tnr.std, P, N);

    % For model Statistics
    stats.model.totSV.mean   = reshape(stats.model.totSV.mean, P, N);
    stats.model.totSV.std    = reshape(stats.model.totSV.std, P, N);

    stats.model.ratioSV.mean = reshape(stats.model.ratioSV.mean, P, N);
    stats.model.ratioSV.std  = reshape(stats.model.ratioSV.std, P, N);

    stats.model.posSV.mean   = reshape(stats.model.posSV.mean, P, N);
    stats.model.posSV.std    = reshape(stats.model.posSV.std, P, N);

    stats.model.negSV.mean   = reshape(stats.model.negSV.mean, P, N);
    stats.model.negSV.std    = reshape(stats.model.negSV.std, P, N);

    stats.model.boundSV.mean = reshape(stats.model.boundSV.mean, P, N);
    stats.model.boundSV.std  = reshape(stats.model.boundSV.std, P, N);

end

end

author: epfl-lasa


text:
% Author: Junaid Ahmed Ansari @ Robotics Research Center, IIIT Hyderabad
% Email : junaid.ansari@research.iiit.ac.in
function generateCandidateFor2D2Dmatching(scoreMat_3D2D, thres, fpath, fno, detectionsQ, detectionsT)

    fname = sprintf('%s%06d.txt',fpath,fno);
    fid   = fopen(fname,'w+');

    Q_dno = zeros(length(detectionsQ), 1);
    T_dno = zeros(length(detectionsT), 1);

    % save detection numbers into array for speed
    for i = 1:length(detectionsQ)
        Q_dno(i) = detectionsQ{i}.dno;
    end

    for i = 1:length(detectionsT)
        T_dno(i) = detectionsT{i}.dno;
    end

    % find candidate matches based on 3D2D for 2D2D matching
    for i = 1:length(detectionsQ)
        i
        cand_match = [];
        cand_match = T_dno(find(scoreMat_3D2D(i,:) > thres))
        cand_match = [detectionsQ{i}.dno; length(cand_match); cand_match]'
        fprintf(fid,'%d ',cand_match);
        fprintf(fid,'\n');
    end

    fclose(fid);
end

repo: MOTBeyondPixels


text:
% textgui() - make sliding vertical window. This window contain text
%             with optional function calls at each line.
%
% Usage:
%   >> textgui( commandnames, helparray, 'key', 'val' ...);
%
% Inputs:
%   commandnames - name of the commands. Either char array or cell
%                  array of char. All style are 'pushbuttons' exept
%                  for empty commands.
%   helparray    - cell array of commands to execute for each menu
%                  (default is empty)
%
% Optional inputs:
%  'title'        - window title
%  'fontweight'   - font weight. Single value or cell array of value for each line.
%  'fontsize'     - font size. Single value or cell array of value for each line.
%  'fontname'     - font name. Single value or cell array of value for each line.
%  'linesperpage' - number of line per page. Default is 20.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% Example:
%  textgui({ 'function1' 'function2' }, {'hthelp(''function1'')' []});
%  % this function will call a pop_up window with one button on
%  % 'function1' which will call the help of this function.
%

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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

function tmp = textgui( textmenu, helparray, varargin);

XBORDERS   = 0; % pixels
TOPORDINATE = 1.09;
TOPLINES    = 2;

if nargin < 1
	help textgui;
	return;
end;
if nargin <2
	helparray = cell(1, length(	textmenu ));
end;

% optional parameters
% -------------------
if nargin >2
    for i = 1:length(varargin)
        if iscell(varargin{i}), varargin(i) = { varargin(i) }; end;
    end;
    g=struct(varargin{:});
else
    g = [];
end;
try, g.title;    catch, g.title = ''; end;
try, g.fontname; catch, g.fontname = 'courier'; end;
try, g.fontsize; catch, g.fontsize = 12; end;
try, g.fontweight; catch, g.fontweight = 'normal'; end;
try, g.linesperpage; catch, g.linesperpage = 20; end;

if isempty( helparray )
	helparray = cell(1,200);
else
	tmpcell = cell(1,200);
	helparray = { helparray{:} tmpcell{:} };
end;

% number of elements
% ------------------
if iscell(textmenu)	nblines = length(textmenu);
else				nblines = size(textmenu,1);
end;

% generating the main figure
% --------------------------
fig = figure('position', [100 100 800 25*15], 'menubar', 'none', 'numbertitle', 'off', 'name', 'Help', 'color', 'w');
pos = get(gca,'position'); % plot relative to current axes
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)]./100;

% generating cell arrays
% ----------------------
if isstr(g.fontname),     tmp1(1:nblines) = {g.fontname}; g.fontname = tmp1; end;
if isnumeric(g.fontsize), tmp2(1:nblines) = {g.fontsize}; g.fontsize = tmp2; end;
if isstr(g.fontweight),   tmp3(1:nblines) = {g.fontweight}; g.fontweight = tmp3; end;
switch g.fontname{1}
    case 'courrier', CHAR_WIDTH = 11; % pixels
    case 'times', CHAR_WIDTH = 11; % pixels
    otherwise, CHAR_WIDTH = 11;
end;

topordi = TOPORDINATE;
if ~isempty(g.title)
    addlines = size(g.title,1)+1+TOPLINES;
    if nblines+addlines < g.linesperpage
        divider = g.linesperpage;
    else
        divider = nblines+addlines;
    end;
    ordinate      = topordi-topordi*TOPLINES/divider;
    currentheight = topordi/divider;

    for index=1:size(g.title,1)
      ordinate      = topordi-topordi*(TOPLINES+index-1)/divider;
      h = text( -0.1, ordinate, g.title(index,:), 'unit', 'normalized', 'horizontalalignment', 'left', ...
		 'fontname', g.fontname{1}, 'fontsize', g.fontsize{1},'fontweight', fastif(index ==1, 'bold', ...
		  'normal'), 'interpreter', 'none' );
    end;
    %h = uicontrol( gcf, 'unit', 'normalized', 'style', 'text', 'backgroundcolor', get(gcf, 'color'), 'horizontalalignment', 'left', ...
    %		'position', [-0.1 ordinate 1.1 currentheight].*s*100+q, 'string', g.title, ...
    %	       'fontsize', g.fontsize{1}+1, 'fontweight', 'bold', 'fontname', g.fontname{1});
 else
    addlines = TOPLINES;
    if nblines+addlines < g.linesperpage
        divider = g.linesperpage;
    else
        divider = nblines+addlines;
    end;
end;

maxlen = size(g.title,2);
for i=1:nblines
	if iscell(textmenu)	tmptext = textmenu{i};
	else			tmptext = textmenu(i,:);
	end;
	ordinate      = topordi-topordi*(i-1+addlines)/divider;
	currentheight = topordi/divider;

	if isempty(helparray{i})
		h = text( -0.1, ordinate, tmptext, 'unit', 'normalized', 'horizontalalignment', 'left', ...
		 'fontname', g.fontname{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i}, 'interpreter', 'none' );
%		h = uicontrol( gcf, 'unit', 'normalized', 'style', 'text', 'backgroundcolor', get(gcf, 'color'), 'horizontalalignment', 'left', ...
%			'position', [-0.1 ordinate 1.1 currentheight].*s*100+q, 'string', tmptext, 'fontname', g.fontname{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i} );
	else
		h = text( -0.1, ordinate, tmptext, 'unit', 'normalized', 'color', 'b', 'horizontalalignment', 'left', ...
		 'fontname', g.fontname{i}, 'buttondownfcn', helparray{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i}, 'interpreter', 'none' );
%		h = uicontrol( gcf,  'callback', helparray{i}, 'unit','normalized', 'style', 'pushbutton',  'horizontalalignment', 'left', ...
%			'position', [-0.1 ordinate 1.1 currentheight].*s*100+q, 'string', tmptext, 'fontname', g.fontname{i}, 'fontsize', g.fontsize{i},'fontweight', g.fontweight{i});
	end;
	maxlen = max(maxlen, length(tmptext));
end;
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1:2) maxlen*CHAR_WIDTH+2*XBORDERS 400]);
% = findobj('parent', gcf);
%set(h, 'fontname', 'courier');

if nblines+addlines < g.linesperpage
    zooming  = 1;
else
    zooming  = (nblines+addlines)/g.linesperpage;
end;

slider(gcf, 0, 1, 1, zooming, 0);
axis off;

author: PatternRecognition

repo: OpenBMI


text:
function [ a, seed ] = npart_rsf_lex_random ( n, npart, seed )

%*****************************************************************************80
%
%% NPART_RSF_LEX_RANDOM returns a random RSF NPART partition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Input, integer NPART, the number of parts of the partition.
%    1 <= NPART <= N.
%
%    Input/output, integer SEED, a seed for the random number
%    generator.
%
%    Output, integer A(NPART), contains the partition.
%    A(1) through A(NPART) contain the nonzero integers which
%    sum to N.
%
  npartitions = npart_enum ( n, npart );

  [ rank, seed ] = i4_uniform ( 1, npartitions, seed );

  a = npart_rsf_lex_unrank ( rank, n, npart );

  return
end


text:
function do_all(dbName)

%DO_ALL execute all codes in the ttimes package.
% DO_ALL(dbName) runs each of the codes in the ttimes package on the
% database specified by dbName. See the ttimes.dbload function for a
% description.
%
% See also ttimres.dbload

ttimes.map(dbName);
ttimes.depth_section(dbName);
ttimes.tt_curve(dbName);
ttimes.arrival_histogram(dbName);
ttimes.write_lotos(dbName);
ttimes.dbload(dbName)

author: geoscience-community-codes

repo: GISMO


text:
function [fnormhat,t]=instfreq(x,t,L,trace);
%INSTFREQ Instantaneous frequency estimation.
%	[FNORMHAT,T]=INSTFREQ(X,T,L,TRACE) computes the instantaneous
%	frequency of the analytic signal X at time instant(s) T, using the
%	trapezoidal integration rule.
%	The result FNORMHAT lies between 0.0 and 0.5.
%
%	X : Analytic signal to be analyzed.
%	T : Time instants	        (default : 2:length(X)-1).
%	L : If L=1, computes the (normalized) instantaneous frequency
%	    of the signal X defined as angle(X(T+1)*conj(X(T-1)) ;
%	    if L>1, computes a Maximum Likelihood estimation of the
%	    instantaneous frequency of the deterministic part of the signal
%	    blurried in a white gaussian noise.
%	    L must be an integer       	(default : 1).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                (default : 0).
%	FNORMHAT : Output (normalized) instantaneous frequency.
%	T : Time instants.
%
%	Examples :
%	 x=fmsin(70,0.05,0.35,25); [instf,t]=instfreq(x); plot(t,instf)
%	 N=64; SNR=10.0; L=4; t=L+1:N-L; x=fmsin(N,0.05,0.35,40);
%	 sig=sigmerge(x,hilbert(randn(N,1)),SNR);
%	 plotifl(t,[instfreq(sig,t,L),instfreq(x,t)]); grid;
%	 title ('theoretical and estimated instantaneous frequencies');
%
%	See also  KAYTTH, SGRPDLAY.

%	F. Auger, March 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least one parameter required');
end;
[xrow,xcol] = size(x);
if (xcol~=1),
 error('X must have only one column');
end

if (nargin == 1),
 t=2:xrow-1; L=1; trace=0.0;
elseif (nargin == 2),
 L = 1; trace=0.0;
elseif (nargin == 3),
 trace=0.0;
end;

if L<1,
 error('L must be >=1');
end
[trow,tcol] = size(t);
if (trow~=1),
 error('T must have only one row');
end;

if (L==1),
 if any(t==1)|any(t==xrow),
  error('T can not be equal to 1 neither to the last element of X');
 else
  fnormhat=0.5*(angle(-x(t+1).*conj(x(t-1)))+pi)/(2*pi);
 end;
else
 H=kaytth(L);
 if any(t<=L)|any(t+L>xrow),
  error('The relation L<T<=length(X)-L must be satisfied');
 else
  for icol=1:tcol,
   if trace, disprog(icol,tcol,10); end;
   ti = t(icol); tau = 0:L;
   R = x(ti+tau).*conj(x(ti-tau));
   M4 = R(2:L+1).*conj(R(1:L));

   diff=2e-6;
   tetapred = H * (unwrap(angle(-M4))+pi);
   while tetapred<0.0 , tetapred=tetapred+(2*pi); end;
   while tetapred>2*pi, tetapred=tetapred-(2*pi); end;
   iter = 1;
   while (diff > 1e-6)&(iter<50),
    M4bis=M4 .* exp(-j*2.0*tetapred);
    teta = H * (unwrap(angle(M4bis))+2.0*tetapred);
    while teta<0.0 , teta=(2*pi)+teta; end;
    while teta>2*pi, teta=teta-(2*pi); end;
    diff=abs(teta-tetapred);
    tetapred=teta; iter=iter+1;
   end;
   fnormhat(icol,1)=teta/(2*pi);
  end;
 end;
end;

author: HeLiangHIT


text:
function results = tensor_rpca_adal2( data, params )
% Solve
%       min_{X,E} \sum_i ||X_(i)||_* + \lambda1*||E||_1
%       s.t. X + E = T

% data.T
% params
% X, V are cell arrays of tensors.
%
% Algorithm: ADAL

tic;
T = data.T;
N = length( size(params.X0) );
X = cell( 1, N );
U = cell( 1, N );
V = cell( 1, N );
for i = 1:N
    X{i} = params.X0;
    V{i} = params.V0{i};
end
E = params.E0;

lambda = params.lambda;
mu = params.mu1;
verbose = params.verbose;
tnorm = norm(T);

for iter = 1:params.max_iter
    % solve X_i's
    for i = 1:N
        [X{i}, junk, U{i}] = tensor_shrinkage( T+mu*V{i}-E, mu, i );
    end

    % solve E
    Ep = E;
    D = cell( 1, N );
    for i = 1:N
        D{i} = X{i} - T - mu*V{i};
    end
    P = -(1/N)*ten_sum_all(D);
    E = shrinkage_t( P, lambda*mu/N );

    % compute optimality stats
    pres = 0;
    tdiff = cell( 1, N );
    for i = 1:N
        tdiff{i} = X{i} + E - T;
    end
    pres = tensor_array_norm(tdiff);

    Ediff = E - Ep;
    pres = pres / (sqrt(N)*tnorm);
    dres = norm(Ediff) / norm(Ep);
    Y = ten_sum_all(X) / N;
    rel_err = norm(Y-data.X) / norm(data.X);

    % print
    if verbose
    fprintf('Iter: %d,   pinf: %3.2e,   dinf: %3.2e,    rel_err: %3.2e\n', iter, pres, dres, rel_err );
    end

    if max(pres, dres) < params.opt_tol
%     if pres < params.opt_tol
        break;
    end

    % update Lagrange multipliers
    for i = 1:N
        V{i} = V{i} - tdiff{i}/mu;
    end

%     if rem(iter,10) == 0
%         mu = max(mu*0.9, params.mu_min);
%     end
end

results.X = Y;
results.E = E;
results.V = V;
results.T = T;
results.U = U;
results.iter = iter;
results.cpu = toc;
results.mu = mu;
results.lambda = lambda;

end

author: andrewssobral

repo: lrslibrary


text:
function r8_atan_test ( )

%*****************************************************************************80
%
%% R8_ATAN_TEST tests R8_ATAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  ntest = 8;

  xtest(1:8) = [ 1.0, 1.0, 0.0, -1.0, -1.0, -1.0,  0.0,  1.0 ];
  ytest(1:8) = [ 0.0, 1.0, 1.0,  1.0,  0.0, -1.0, -1.0, -1.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_ATAN_TEST\n' );
  fprintf ( 1, '  R8_ATAN computes the arc-tangent given Y and X;\n' );
  fprintf ( 1, '  ATAN2 is the system version of this routine.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X     Y     ATAN2(Y,X)   R8_ATAN(Y,X)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest
    x = xtest(i);
    y = ytest(i);
    fprintf ( 1, '  %12f  %12f  %12f  %12f\n', ...
      x, y, atan2 ( y, x ), r8_atan ( y, x ) );
  end

  return
end


text:
function  res = partialDCT(n,m,J)

res.adjoint = 0;
res.n = n;
res.m = m;
res.J = J(1:m);

% Register this variable as a partialDCT class
res = class(res,'partialDCT');

author: thomaskuestner


text:
function example2 ( )

%*****************************************************************************80
%
%% EXAMPLE2 uses BVP4C to solve the EXAMPLE2 problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    Original MATLAB version by Shampine, Kierzenka, Reichelt.
%    This version by John Burkardt.
%
%  Reference:
%
%    Lawrence Shampine, Jacek Kierzenka, Mark Reichelt,
%    Solving boundary value problems for ordinary differential equations
%    in MATLAB with bvp4c.
%
  global p

  p = 0.00001;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE2:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use BVP4C to solve the following boundary value problem:\n' );
  fprintf ( 1, '  y" + 3 p y / ( p + t^2 )^2 = 0\n' );
  fprintf ( 1, '  y(-0.1) = -0.1 / sqrt ( p + 0.01 )\n' );
  fprintf ( 1, '  y(+0.1) = +0.1 / sqrt ( p + 0.01 )\n' );
  fprintf ( 1, '  Use initial guesses:\n' );
  fprintf ( 1, '  y(x) = 0\n' );
  fprintf ( 1, '  y''(x) = 10\n' );
  fprintf ( 1, '  Use parameter value p = %g\n', p );
%
%  Set SOLINIT, the structure defining the initial guess.
%
  x_init = linspace ( -0.1, 0.1, 11 );
  y_init = [ 0.0, 10.0 ];
  solinit = bvpinit ( x_init, y_init );
%
%  Because this problem is difficult to solve accurately, we must request
%  a more stringent tolerance.  This is done through the optional OPTIONS
%  input.
%
  options = bvpset ( 'RelTol', 1.0E-04 );
%
%  Have BVP4C solve the problem.
%
  sol = bvp4c ( @example2_ode, @example2_bc, solinit, options );
%
%  Use DEVAL to evaluate the solution.
%
  x = linspace ( -0.1, 0.1, 101 );
  sol = deval ( sol, x );

  exact = x ./ sqrt ( p + x.^2 );
%
%  Display a plot of Y and the exact solution.
%
  plot ( x, sol(1,:), 'b*', ...
         x, exact, 'r-', 'Linewidth', 2 );
  xlabel ( '<--- X --->', 'Fontsize', 16 );
  ylabel ( '<--- Y(X) --->', 'Fontsize', 16 );
  title ( 'EXAMPLE2: Exact solution (red), Estimate (blue)', 'Fontsize', 16 )
  grid on
  filename = 'example2.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot file as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE2:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function dydx = example2_ode ( x, y )

%*****************************************************************************80
%
%% EXAMPLE2_ODE evaluates the right hand side of the ODE.
%
%  Discussion:
%
%    We assume that the differential equation has been rewritten as a
%    system of first order equations of the form
%
%      dydx = f(x,y)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which the ODE is to be evaluated.
%
%    Input, real Y(M), the value of the solution at X.
%
%    Output, real DYDX(M), the value of the right hand side given X and Y.
%
  global p

  dydx(1) = y(2);
  dydx(2) = - 3.0 * p * y(1) / ( p + x^2 )^2;

  return
end
function bc = example2_bc ( ya, yb )

%*****************************************************************************80
%
%% EXAMPLE2_BC evaluates the boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real YA(M), YB(M), the solution value at the left and right endpoints.
%
%    Output, real BC(M), the value of the boundary conditions.
%
  global p

  bc(1) = ya(1) + 0.1 / sqrt ( p + 0.01 );
  bc(2) = yb(1) - 0.1 / sqrt ( p + 0.01 );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
% LTFAT - Signal processing tools
%
%  Peter L. Søndergaard, 2007 - 2018.
%
%  General
%    RMS            -  Root Mean Square norm of signal.
%    SETNORM      -  Normalize signal by specified norm.
%    GAINDB         -  Scale input signal
%    CRESTFACTOR    -  Compute the crest factor of a signal.
%    UQUANT         -  Simulate uniform quantization.
%
%  Window functions
%    FIRWIN         -  FIR windows (Hanning,Hamming,Blackman,...).
%    FIRKAISER      -  FIR Kaiser-Bessel window.
%    FIR2LONG       -  Extend FIR window to LONG window.
%    LONG2FIR       -  Cut LONG window to FIR window.
%    FREQWIN        -  Freq responses (Gauss,Gammatone,Butterworth,...)
%
%  Wavelet functions
%    FREQWAVELET    -  Frequency responses of wavelets (Cauchy, Morse)
%
%  Filtering
%    FIRFILTER      -  Construct an FIR filter.
%    BLFILTER       -  Construct a band-limited filter.
%    WARPEDBLFILTER -  Warped, band-limited filter.
%    FREQFILTER     -  Construct a full length frequency side filter.
%    PFILT          -  Apply filter with periodic boundary conditions.
%    MAGRESP        -  Magnitude response plot.
%    TRANSFERFUNCTION - Compute the transfer function of a filter.
%    PGRPDELAY      -  Periodic Group Delay
%
%  Ramping
%    RAMPUP         -  Rising ramp.
%    RAMPDOWN       -  Falling ramp.
%    RAMPSIGNAL     -  Ramp a signal.
%
%  Thresholding methods
%    THRESH         -  Coefficient thresholding.
%    LARGESTR       -  Keep largest ratio of coefficients.
%    LARGESTN       -  Keep N largest coefficients.
%    DYNLIMIT       -  Limit the dynamical range.
%    GROUPTHRESH    -  Group thresholding.
%
%  Image processing
%    RGB2JPEG       -  Convert RGB values to the JPEG colour model
%    JPEG2RGB       -  Convert values from the JPEG colour model to RGB
%
%  Tools for OFDM
%    QAM4           -  Quadrature amplitude modulation, order 4
%    IQAM4          -  Inverse QAM of order 4
%
%  Tonal-transient separation
%    tfjigsawsep      - Tonal-transient-residual separation using the T-F jigsaw puzzle algorithm.
%    plottfjigsawsep  - Plot the separated layers.
%
%  For help, bug reports, suggestions etc. please visit
%  http://github.com/ltfat/ltfat/issues
%

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [varargout] = triangleRayIntersection (V_ori,R,V,F,optStruct)

% function [V_intersect,L_intersect,T] = triangleRayIntersection (V_ori,R,V,F,optStruct)
% ------------------------------------------------------------------------
%
% Ray/triangle intersection using the algorithm proposed by Muller and
% Trumbore (1997)
%
% Based on implementation by: Jarek Tuszynski (jaroslaw.w.tuszynski@saic.com)
%
%
% Algorithm:
%  Function solves
%        |t|
%    M * |u| = (o-v0)
%        |v|
%  for [t; u; v] where M = [-d, v1-v0, v2-v0]. u,v are barycentric coordinates
%  and t - the distance from the ray origin in |d| units
%  ray/triangle intersect if u>=0, v>=0 and u+v<=1
%
% Note:
%  The algorithm is able to solve several types of problems:
%  * many faces / single ray  intersection
%  * one  face  / many   rays intersection
%  * one  face  / one    ray  intersection
%  * many faces / many   rays intersection
%  In order to allow that to happen all imput arrays are expected in Nx3
%  format, where N is number of vertices or rays. In most cases number of
%  vertices is different than number of rays, so one of the imputs will
%  have to be cloned to have the right size. Use "repmat(A,size(B,1),1)".
%
% Input (all arrays in in Nx3 format, where N is number of vertices or rays):
%  * orig : ray's origin
%  * dir  : ray's direction
%  * vert0, vert1, vert2: vertices of the triangle
%  * options: aditional customization options
%    * options.triangle - 'one sided' or 'two sided' (default) - how to treat
%        triangles. In 'one sided' version only intersections in single
%        direction are counted and intersections with back facing
%           tringles are ignored
%    * options.ray - 'ray' (default) or 'segment' - how to treat ray as an
%        infinite line (ray) or as line segment defined by a vector
%    * option.border - controls border handling. If 'normal'(default)
%        border points are included, but can be easily lost due to
%        rounding errors. If option.border='inclusive' borders points are
%        included, with a margin of option.eps. If option.border='exclusive'
%        borders points are excluded, with margin of option.eps.
%    * options.epsilon (default = 1e-5)
%
% Output:
%   * Intersect - boolean array of length N
%   * t   - distance from the ray origin to the intersection point in |dir|
%   * u,v - barycentric coordinates of the intersection point units
%
% Based on:
%  *"Fast, minimum storage ray-triangle intersection". Tomas Muller and
%    Ben Trumbore. Journal of Graphics Tools, 2(1):21--28, 1997.
%    http://www.graphics.cornell.edu/pubs/1997/MT97.pdf
%  * http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/
%  * http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/raytri.c
%
% Change log:
% 2013/12/18 Changed for patch type input and Cartesian coordinate output
% 2019/07/03 Added varargout type output
% 2021/08/18 Listed as depricated in favour of triSurfRayTrace
% ------------------------------------------------------------------------

%%

warning('This function is depricated, it is recommended to update your code to use triSurfRayTrace instead');

%%

%Get triangle vertices
V0=V(F(:,1),:);
V1=V(F(:,2),:);
V2=V(F(:,3),:);

%% verify that inputs are in correct format

if (size(V_ori,1)==3 && size(V_ori,2)~=3)
    V_ori =V_ori';
end

if (size(R,1)==3 && size(R,2)~=3)
    R =R';
end

if (size(V0,1)==3 && size(V0,2)~=3)
    V0=V0';
end

if (size(V1,1)==3 && size(V1,2)~=3)
    V1=V1';
end

if (size(V2,1)==3 && size(V2,2)~=3)
    V2=V2';
end

if (any(size(V_ori)~=size(V0)) || ...
        any(size(V_ori)~=size(V1)) || ...
        any(size(V_ori)~=size(V2)) || ...
        any(size(V_ori)~=size(R  )) )
    error('All input vectors have to be of the same size.');
end

if (size(V_ori,2)~=3)
    error('All input vectors have to be in Nx3 format.');
end

%% Read options from optStruct

%Create defaults
if (nargin<5)
    optStruct.eps      = 1e-5;
    optStruct.triangle = 'two sided';
    optStruct.ray      = 'ray';
    optStruct.border   = 'normal';
end

%% initialize default output
L_intersect = false(size(V_ori,1),1);
T=nan(size(V_ori,1),1);
V_bar=T;
V_intersect=nan(size(V_ori));

%% Find faces parallel to the ray
E1 = V1-V0;          % find vectors for two edges sharing vert0
E2 = V2-V0;
Tvec  = V_ori -V0;          % distance from vert0 to ray origin
Pvec  = cross(R, E2,2);  % begin calculating determinant - also used to calculate U parameter

detMat   = sum(E1.*Pvec,2);   % determinant of the matrix M = dot(E1,Pvec)

logicParalel = (abs(detMat)<optStruct.eps);    % if determinant is near zero then ray lies in the plane of the triangle

if all(logicParalel) % if all parallel than no intersections
    varargout{1}=V_intersect;
    varargout{2}=L_intersect;
    varargout{3}=T;
    return
end

switch optStruct.border
    case 'normal'
        zeroLim=0.0;
    case 'inclusive'
        zeroLim=optStruct.eps;
    case 'exclusive'
        zeroLim=-optStruct.eps;
end

%% Different behavior depending on one or two sided triangles
if strcmpi(optStruct.triangle,'two sided')          % treats triangles as two sided

    detMat(logicParalel) = 1;                       % change to avoid division by zero

    U_bar  = sum(Tvec.*Pvec,2)./detMat;             % calculate U parameter used to test bounds
    L_ok = (~logicParalel & U_bar>=-zeroLim & U_bar<=1.0+zeroLim);% mask which allows performing next 2 operations only when needed

    if ~any(L_ok) % if all ray/plane intersections are outside the triangle than no intersections
        varargout{1}=V_intersect;
        varargout{2}=L_intersect;
        varargout{3}=T;
        return
    end

    Qvec = cross(Tvec(L_ok,:), E1(L_ok,:),2); % prepare to test V parameter
    V_bar(L_ok,:) = sum(R(L_ok,:).*Qvec,2) ./ detMat(L_ok,:);  % calculate V parameter used to test bounds
    L_intersect = (V_bar>=-zeroLim & U_bar+V_bar<=1.0+zeroLim & L_ok);

    T(L_ok,:) = sum(E2(L_ok,:).*Qvec,2)./detMat(L_ok,:);

    if ~(strcmpi(optStruct.ray,'ray'))
        L_intersect = (L_intersect & T>=-zeroLim & T<=1.0+zeroLim);
    end

else % treats triangles as one sided

    U_bar = sum(Tvec.*Pvec,2);                   % calculate U parameter used to test bounds
    L_ok = (detMat>optStruct.eps & U_bar>=0.0 & U_bar<=detMat);        % mask which allows performing next 2 operations only when needed

    if ~any(L_ok) % if all ray/plane intersections are outside the triangle than no intersections
        varargout{1}=V_intersect;
        varargout{2}=L_intersect;
        varargout{3}=T;
        return
    end

    Qvec = cross(Tvec(L_ok,:), E1(L_ok,:),2); % prepare to test V parameter
    V_bar(L_ok,:) = sum(R(L_ok,:).*Qvec,2);        % calculate V parameter used to test bounds
    L_intersect = (detMat>optStruct.eps & U_bar>=-zeroLim & V_bar>=-zeroLim & U_bar+V_bar<=detMat*(1+zeroLim));

    T(L_ok,:)  = sum(E2(L_ok,:).*Qvec,2);
    inv_det = zeros(size(detMat));
    inv_det(L_ok,:) = 1./detMat(L_ok,:);
    T = T.*inv_det;  % calculate t - distance from origin to the intersection in |d| units
    U_bar = U_bar.*inv_det;
    V_bar = V_bar.*inv_det;

    if ~(strcmpi(optStruct.ray,'ray'))
        L_intersect = (L_intersect & T>=-zeroLim & T<=1.0+zeroLim); % intersection between origin and destination
    end
end

V_intersect=((1-U_bar(:,ones(1,3))-V_bar(:,ones(1,3))).*V0) + (U_bar(:,ones(1,3)).*V1) + (V_bar(:,ones(1,3)).*V2);
V_intersect(~L_intersect,:)=NaN;
T(~L_intersect,:)=NaN;

%% Collect output
varargout{1}=V_intersect;
varargout{2}=L_intersect;
varargout{3}=T;

%%
% _*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 data_num = p07_data_num ( )

%*****************************************************************************80
%
%% P07_DATA_NUM returns the dimension of the data vector for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the dimension of the data vector.
%
  data_num = 9;

  return
end


text:
function h = dtiImportFibers(h, filename, xform)
%
% handles = dtiImportFibers(handles, filename)
%
% Imports non-mrDiffusion fiber/path file formats. Currently supports:
%
% 1. MetroTrac (*.pdb) format
% 2. Camino Raw (*.Bfloat) format
%
%
% HISTORY:
% 2006.08.17 RFD: wrote it.
%

fg = mtrImportFibers(filename, xform);
if ~isempty(fg)
    h = dtiAddFG(fg, h);
end
return;


text:
function g = disk_grid_fibonacci ( n, r, c )

%*****************************************************************************80
%
%% DISK_GRID_FIBONACCI computes Fibonacci grid points inside a disk.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Swinbank, James Purser,
%    Fibonacci grids: A novel approach to global modelling,
%    Quarterly Journal of the Royal Meteorological Society,
%    Volume 132, Number 619, July 2006 Part B, pages 1769-1793.
%
%  Parameters:
%
%    Input, integer N, the number of points desired.
%
%    Input, real R, the radius of the disk.
%
%    Input, real C(2), the coordinates of the center of the disk.
%
%    Output, real G(2,N), the grid points.
%
  r0 = r / sqrt ( n - 0.5 );
  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0;

  i = 1 : n;

  gr = r0 * sqrt ( i - 0.5 );
  gt = 2.0 * pi * i / phi;

  g = zeros ( 2, n );

  g(1,1:n) = c(1);
  g(1,1:n) = g(1,1:n) + gr .* cos ( gt );
  g(2,1:n) = c(2);
  g(2,1:n) = g(2,1:n) + gr .* sin ( gt );

  return
end


text:
function [xi,w]=fifthOrderHexagonCubPoints()
%%FIFTHORDERHEXAGONCUBPOINTS Generate fifth-order cubature points for
%               integration over a two-dimensional hexagon with vertices
%               (+/-1,0) and (+/-(1/2),+/-(1/2)*sqrt(3)).
%
%INPUTS: None
%
%OUTPUTS: xi A numDim X numCubaturePoints matrix containing the cubature
%            points (Each "point" is a vector).
%          w A numCubaturePoints X 1 vector of the weights associated with
%            the cubature points.
%
%This function implements Formula H_2 5-1 in [1], pg. 336, 7 points.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

r=sqrt(14/25);
s=sqrt(14/100);
t=sqrt(42/100);

xi=[[0;0],PMCombos([r;0]),PMCombos([s;t])];

V=3*sqrt(3)/2;
A=(258/1008)*V;
B=(125/1008)*V;
w=[A;B*ones(6,1)];
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Ax,Ay,b,K] = convexhullFromSampled(z,fz,xL,xU);

k1 = max((fz(2:end)-fz(1))./(z(2:end)-xL))+1e-12;
k2 = min((fz(2:end)-fz(1))./(z(2:end)-xL))-1e-12;
k3 = min((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))+1e-12;
k4 = max((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))-1e-12;
Ax = [-k1;k2;-k3;k4];
Ay = [1;-1;1;-1];
b =  [k1*(-z(1)) + fz(1);-(k2*(-z(1)) + fz(1));k3*(-z(end)) + fz(end);-(k4*(-z(end)) + fz(end))];
K.f = 0;
K.l = length(b);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [Hx,Hy,Ez,time] = MaxwellHNonCon2D(Hx, Hy, Ez, FinalTime)

% function [Hx,Hy,Ez] = MaxwellHNonCon2D(Hx, Hy, Ez, FinalTime)
% Purpose  : Integrate TM-mode Maxwell's until FinalTime starting with initial conditions Hx,Hy,Ez
%

Globals2D;
time = 0;

% Runge-Kutta residual storage
resHx = zeros(Np,K); resHy = zeros(Np,K); resEz = zeros(Np,K);

% compute time step size
rLGL = JacobiGQ(0,0,N); rmin = abs(rLGL(1)-rLGL(2));
dtscale = dtscale2D; dt = min(dtscale)*rmin*2/3;

% find non-conforming neighbors
neighbors = BuildHNonCon2D(N+1, 1e-6);

% outer time step loop
tstep = 1;
while (time<FinalTime)

  if(time+dt>FinalTime), dt = FinalTime-time; end

   for INTRK = 1:5
      % compute right hand side of TM-mode Maxwell's equations
      [rhsHx, rhsHy, rhsEz] = MaxwellHNonConRHS2D(Hx,Hy,Ez, neighbors);

      % initiate and increment Runge-Kutta residuals
      resHx = rk4a(INTRK)*resHx + dt*rhsHx;
      resHy = rk4a(INTRK)*resHy + dt*rhsHy;
      resEz = rk4a(INTRK)*resEz + dt*rhsEz;

      % update fields
      Hx = Hx+rk4b(INTRK)*resHx;
      Hy = Hy+rk4b(INTRK)*resHy;
      Ez = Ez+rk4b(INTRK)*resEz;
   end;

   % Increment time
   time = time+dt; tstep = tstep+1;

   if(~mod(tstep, 5))
     PlotField2D(N, x, y, Ez);
     hold on; ha = trimesh(EToV, VX, VY); set(ha, 'Color', 'black'); hold off
     axis([-1 1 -1 1 -1 1]); drawnow; pause(.02);
   end

end
return

repo: Aero-matlab


text:
function val=erfcxComplex(z)
%%ERFCXCOMPLEX Evaluate the scaled complementary error function in a manner
%             valid for complex arguments. The complementary error function
%             is one minus the error function (erf). The scaled
%             complementary error function is exp(z^2)*erfc(z). The
%             functions erfcx and erfc that are built into Matlab cannot
%             handle complex arguments.
%
%INPUTS: z A scalar, vector, or matrix of values at which one wishes to
%          evaluate the scaled complementary error function.
%
%OUTPUTS: val A matrix having the same dimensions as z in which the values
%            of the scaled complementary error function are computed for
%            the values in z.
%
%The scaled complementary error function is just Faddeeva(1i*z) based on
%the definition of the Faddeeva function.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

val=Faddeeva(1i*z);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function [g,NR,SI,TI] = regiongrow(f,S,T)
%REGIONGROW Image segmentation using region growing.
%   [G,NR,SI,TI] = REGIONGROW(F,S,T). S can be a matrix of the same size
%   as the input image F, with a 1 at the coordinates of every seed
%   point and 0s elsewhere. S can also be a single seed value. Similarly
%   T can be a matrix of the same size as F, containing a threshold
%   value for each pixel in F. T can also be a scalar, in which case it
%   becomes a global threshold. All values in S and T must be in the
%   range [0,1].
%
%   G is the segmented image resulting from region growing, with each
%   region labeled by a different integer, NR is the number of regions,
%   SI is the final seed image used by the algorithm, and TI is the
%   image consisting of the pixels in F that satisfied the threshold
%   test, but before they were processed for connectivity.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

f = tofloat(f);
% If S is a scalar, obtain the seed image.
if numel(S) == 1
   SI = f == S;
   S1 = S;
else
   % S is a matrix. Eliminate duplicate, connected seed locations to
   % reduce the number of loop executions in the following section of
   % code.
   SI = bwmorph(S,'shrink',Inf);
   S1 = f(SI); % Matrix of seed values.
end

TI = false(size(f));
for K = 1:length(S1)
   seedvalue = S1(K);
   S = abs(f - seedvalue) <= T; % Re-use variable S.
   TI = TI | S;
end

% Use function imreconstruct with SI as the marker image to obtain the
% regions corresponding to each seed in S. Function bwlabel assigns a
% different integer to each connected region.
[g,NR] = bwlabel(imreconstruct(SI,TI));

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function [lp,dlp] = priorTransform(g,dg,ig,prior,x)

% Transformed univariate hyperparameter prior distribution.
% Compute log-likelihood and its derivative or draw a random sample.
% The prior distribution is parameterized as:
%
%  prior(y), where y = g(x),
%
% g is the transformation function (monotonically increasing and invertable),
% dg is its derivative, ig is its inverse, prior is the distribution to be
% transformed and x(1xN) contains query hyperparameters for prior evaluation.
% Note that p(x) = prior(g(x)) is not necessarily normalised w.r.t. x.
%
% For example, to put a nonnegative gamma prior on the lenthscale
% parameter ell of covSEiso, one can use
%   {@priorTransform,@exp,@exp,@log,{@priorGamma,k,t}}
% since ell is represented in the log domain.
%
% For more help on design of priors, try "help priorDistributions".
%
% Copyright (c) by Roman Garnett and Hannes Nickisch, 2014-09-08.
%
% See also PRIORDISTRIBUTIONS.M.

if nargin<4, error('g, dg, ig and prior parameters need to be provided'), end
if nargin<5, lp = ig(feval(prior{:})); return, end      % apply inverse sampling

if ~iscell(prior), prior = {prior}; end
[lp,dlp] = feval(prior{:},g(x)); dlp = dlp.*dg(x);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa