text:
function simNoise = generateNoise(nPoints, sizeVec)
% simNoise = generateNoise(nPoints, sizeVec)
%
% generates Poisson noise with given number of points (nPoints) and given
% image size (sizeVec = [xlim1, xlim2, ylim1, ylim2]).

simNoise(:,1) = sizeVec(1) + rand(nPoints,1)*(sizeVec(2) - sizeVec(1));
simNoise(:,2) = sizeVec(3) + rand(nPoints,1)*(sizeVec(4) - sizeVec(3));

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function p = triangle_circumcenter ( n, t )

%*****************************************************************************80
%
%% TRIANGLE_CIRCUMCENTER computes the circumcenter of a triangle in ND.
%
%  Discussion:
%
%    Three ND points A, B and C lie on a circle.
%
%    The circumcenter P has the formula
%
%      P = ( Area ( PBC ) * A + Area ( APC) * B + Area ( ABP ) * C )
%        / ( Area ( PBC )     + Area ( APC )    + Area ( ABP ) )
%
%    The details of the formula rely on information supplied
%    by Oscar Lanzi III.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real T(N,3), the triangle vertices.
%
%    Output, real P(N,1), the circumcenter of the triangle.
%
  a = r8vec_normsq_affine ( n, t(1:n,2), t(1:n,3) );
  b = r8vec_normsq_affine ( n, t(1:n,3), t(1:n,1) );
  c = r8vec_normsq_affine ( n, t(1:n,1), t(1:n,2) );

  pbc = a * ( - a + b + c );
  apc = b * (   a - b + c );
  abp = c * (   a + b - c );

  p = zeros ( n, 1 );

  p(1:n,1) = ( pbc * t(1:n,1) + apc * t(1:n,2) + abp * t(1:n,3) ) ...
           / ( pbc            + apc            + abp );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = rodman_determinant ( n, alpha )

%*****************************************************************************80
%
%% RODMAN_DETERMINANT returns the determinant of the RODMAN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real ALPHA, the parameter.
%
%    Output, real VALUE, the determinant.
%
  value = ( 1.0 - alpha )^( n - 1 ) * ( 1.0 + alpha * ( n - 1 ) );

  return
end


%% APAP-conssite
% %{
imfolder = 'images\APAP-conssite';
im_n = 5;
imfile = cell(im_n,1);
imfile{1} = [imfolder '\' 'DSC03002.jpg'];
imfile{2} = [imfolder '\' 'DSC03003.jpg'];
imfile{3} = [imfolder '\' 'DSC03004.jpg'];
imfile{4} = [imfolder '\' 'DSC03005.jpg'];
imfile{5} = [imfolder '\' 'DSC03006.jpg'];

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = [1,2; 1,3; 1,4; 2,3; 2,4; 2,5; 3,4; 3,5; 4,5];

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 1, 'equi', 0, imfolder );
%}


text:
function [Y, R, E] = Isomap(D, n_fcn, n_size, options);

% ISOMAP   Computes Isomap embedding using the algorithm of
%             Tenenbaum, de Silva, and Langford (2000).
%
% [Y, R, E] = isomap(D, n_fcn, n_size, options);
%
% Input:
%    D = N x N matrix of distances (where N is the number of data points)
%    n_fcn = neighborhood function ('epsilon' or 'k')
%    n_size = neighborhood size (value for epsilon or k)
%
%    options.dims = (row) vector of embedding dimensionalities to use
%                        (1:10 = default)
%    options.comp = which connected component to embed, if more than one.
%                        (1 = largest (default), 2 = second largest, ...)
%    options.display = plot residual variance and 2-D embedding?
%                        (1 = yes (default), 0 = no)
%    options.overlay = overlay graph on 2-D embedding?
%                        (1 = yes (default), 0 = no)
%    options.verbose = display progress reports?
%                        (1 = yes (default), 0 = no)
%
% Output:
%    Y = Y.coords is a cell array, with coordinates for d-dimensional embeddings
%         in Y.coords{d}.  Y.index contains the indices of the points embedded.
%    R = residual variances for embeddings in Y
%    E = edge matrix for neighborhood graph
%

%    BEGIN COPYRIGHT NOTICE
%
%    Isomap code -- (c) 1998-2000 Josh Tenenbaum
%
%    This code is provided as is, with no guarantees except that
%    bugs are almost surely present.  Published reports of research
%    using this code (or a modified version) should cite the
%    article that describes the algorithm:
%
%      J. B. Tenenbaum, V. de Silva, J. C. Langford (2000).  A global
%      geometric framework for nonlinear dimensionality reduction.
%      Science 290 (5500): 2319-2323, 22 December 2000.
%
%    Comments and bug reports are welcome.  Email to jbt@psych.stanford.edu.
%    I would also appreciate hearing about how you used this code,
%    improvements that you have made to it, or translations into other
%    languages.
%
%    You are free to modify, extend or distribute this code, as long
%    as this copyright notice is included whole and unchanged.
%
%    END COPYRIGHT NOTICE

%%%%% Step 0: Initialization and Parameters %%%%%

N = size(D,1);
if ~(N==size(D,2))
     error('D must be a square matrix');
end;
if n_fcn=='k'
     K = n_size;
     if ~(K==round(K))
         error('Number of neighbors for k method must be an integer');
     end
elseif n_fcn=='epsilon'
     epsilon = n_size;
else
     error('Neighborhood function must be either epsilon or k');
end
if nargin < 3
     error('Too few input arguments');
elseif nargin < 4
     options = struct('dims',1:10,'overlay',1,'comp',1,'display',1,'verbose',1);
end
INF =  1000*max(max(D))*N;  %% effectively infinite distance

if ~isfield(options,'dims')
     options.dims = 1:10;
end
if ~isfield(options,'overlay')
     options.overlay = 1;
end
if ~isfield(options,'comp')
     options.comp = 1;
end
if ~isfield(options,'display')
     options.display = 1;
end
if ~isfield(options,'verbose')
     options.verbose = 1;
end
dims = options.dims;
comp = options.comp;
overlay = options.overlay;
displ = options.display;
verbose = options.verbose;

Y.coords = cell(length(dims),1);
R = zeros(1,length(dims));

%%%%% Step 1: Construct neighborhood graph %%%%%
disp('Constructing neighborhood graph...');

if n_fcn == 'k'
     [tmp, ind] = sort(D);
     for i=1:N
          D(i,ind((2+K):end,i)) = INF;
     end
elseif n_fcn == 'epsilon'
     warning off    %% Next line causes an unnecessary warning, so turn it off
     D =  D./(D<=epsilon);
     D = min(D,INF);
     warning on
end

D = min(D,D');    %% Make sure distance matrix is symmetric

if (overlay == 1)
     E = int8(1-(D==INF));  %%  Edge information for subsequent graph overlay
end

% Finite entries in D now correspond to distances between neighboring points.
% Infinite entries (really, equal to INF) in D now correspond to
%   non-neighoring points.

%%%%% Step 2: Compute shortest paths %%%%%
disp('Computing shortest paths...');

% We use Floyd's algorithm, which produces the best performance in Matlab.
% Dijkstra's algorithm is significantly more efficient for sparse graphs,
% but requires for-loops that are very slow to run in Matlab.  A significantly
% faster implementation of Isomap that calls a MEX file for Dijkstra's
% algorithm can be found in isomap2.m (and the accompanying files
% dijkstra.c and dijkstra.dll).

tic;
for k=1:N
     D = min(D,repmat(D(:,k),[1 N])+repmat(D(k,:),[N 1]));
     if ((verbose == 1) & (rem(k,20) == 0))
          disp([' Iteration: ' num2str(k) '     Estimated time to completion: 'num2str((N-k)*toc/k/60) ' minutes']);
     end
end

%%%%% Remove outliers from graph %%%%%
disp('Checking for outliers...');
n_connect = sum(~(D==INF));        %% number of points each point connects to
[tmp, firsts] = min(D==INF);       %% first point each point connects to
[comps, I, J] = unique(firsts);    %% represent each connected component once
size_comps = n_connect(comps);     %% size of each connected component
[tmp, comp_order] = sort(size_comps);  %% sort connected components by size
comps = comps(comp_order(end:-1:1));
size_comps = size_comps(comp_order(end:-1:1));
n_comps = length(comps);               %% number of connected components
if (comp>n_comps)
     comp=1;                              %% default: use largest component
end
disp(['  Number of connected components in graph: ' num2str(n_comps)]);
disp(['  Embedding component ' num2str(comp) ' with ' num2str(size_comps(comp)) ' points.']);
Y.index = find(firsts==comps(comp));

D = D(Y.index, Y.index);
N = length(Y.index);

%%%%% Step 3: Construct low-dimensional embeddings (Classical MDS) %%%%%
disp('Constructing low-dimensional embeddings (Classical MDS)...');

opt.disp = 0;
[vec, val] = eigs(-.5*(D.^2 - sum(D.^2)'*ones(1,N)/N - ones(N,1)*sum(D.^2)/N + sum(sum(D.^2))/(N^2)), max(dims), 'LR', opt);

h = real(diag(val));
[foo,sorth] = sort(h);  sorth = sorth(end:-1:1);
val = real(diag(val(sorth,sorth)));
vec = vec(:,sorth);

D = reshape(D,N^2,1);
for di = 1:length(dims)
     if (dims(di)<=N)
         Y.coords{di} = real(vec(:,1:dims(di)).*(ones(N,1)*sqrt(val(1:dims(di)))'))';
         r2 = 1-corrcoef(reshape(real(L2_distance(Y.coords{di}, Y.coords{di})),N^2,1),D).^2;
         R(di) = r2(2,1);
         if (verbose == 1)
             disp(['  Isomap on ' num2str(N) ' points with dimensionality ' num2str(dims(di)) '  --> residual variance = ' num2str(R(di))]);
         end
     end
end

clear D;

%%%%%%%%%%%%%%%%%% Graphics %%%%%%%%%%%%%%%%%%

if (displ==1)
     %%%%% Plot fall-off of residual variance with dimensionality %%%%%
     figure;
     hold on
     plot(dims, R, 'bo');
     plot(dims, R, 'b-');
     hold off
     ylabel('Residual variance');
     xlabel('Isomap dimensionality');

     %%%%% Plot two-dimensional configuration %%%%%
     twod = find(dims==2);
     if ~isempty(twod)
         figure;
         hold on;
         plot(Y.coords{twod}(1,:), Y.coords{twod}(2,:), 'ro');
         if (overlay == 1)
             gplot(E(Y.index, Y.index), [Y.coords{twod}(1,:); Y.coords{twod}(2,:)]');
             title('Two-dimensional Isomap embedding (with neighborhood graph).');
         else
             title('Two-dimensional Isomap.');
         end
         hold off;
     end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function YamlStruct = ReadYaml(yaml_file)
% This function reads Yaml file into struct
% Example
% >> yaml_file = 'EnaspolMain.yaml';
% >> YamlStruct = ReadYaml(yaml_file)
%
%  %======================================================================
%{
		Copyright (c) 2011
		This program is a result of a joined cooperation of Energocentrum
		PLUS, s.r.o. and Czech Technical University (CTU) in Prague.
        The program is maintained by Energocentrum PLUS, s.r.o. and
        licensed under the terms of MIT license. Full text of the license
        is included in the program release.

        Author(s):
		Jiri Cigler, Dept. of Control Engineering, CTU Prague
		Jan  Siroky, Energocentrum PLUS s.r.o.

        Implementation and Revisions:

        Auth  Date        Description of change
        ----  ---------   -------------------------------------------------
        jc    01-Mar-11   First implementation
        jc    02-Mar-11   .jar package initialization moved to external fun
        jc    18-Mar-11   Warning added when imported file not found
%}
%======================================================================

InitYaml();

import('org.yaml.snakeyaml.Yaml');

yamlreader = Yaml();
yml = fileread(yaml_file);
jymlobj = yamlreader.load(yml);

Data = Hash2Struct(jymlobj);
work_folder=fileparts(yaml_file);
if isfield(Data,'import')
    yml = '';
    for i=1:numel(Data.import)
        fToImport=[work_folder filesep Data.import{i}];
        try
            fxx=fileread(fToImport);
            yml=sprintf('%s\n%s',yml,  fxx);
        catch
            warning('YAMLMatlab:FileNotFoundException','YAMLMatlab: File %s not found',fToImport);
        end
    end
    jymlobj = yamlreader.load(yml);
    YamlStruct = Hash2Struct(jymlobj);
else
    YamlStruct  =Data;
end
end % end of function

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function kern = whiteKernExpandParam(kern, params)

% WHITEKERNEXPANDPARAM Create kernel structure from WHITE kernel's parameters.
% FORMAT
% DESC returns a white noise kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : whiteKernParamInit, whiteKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

kern.variance = params(1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [present]  = isReactionInSubSystem(model,reactions,subSystem)
% Determine whether a reaction is in a given subSystem.
%
% USAGE:
%
%    [subSystems]  = getModelSubSystems(model)
%
% INPUT:
%    model:                 A COBRA model struct with at least rxns and
%                           subSystems fields
%    reactions:             Either a string identifying a reaction, or a
%                           cell array of strings identifying multiple
%                           reactions, or a double vector identifying the
%                           positions.
%    subSystem:             A String identifying a subsystem.
% OUTPUT:
%    present:               a boolean vector for each provided reaction.
%
% USAGE:
%    %Get all subSystems present in the model.
%    [subSystems]  = getModelSubSystems(model)
%
% .. Author: - Thomas Pfau Nov 2017

if ischar(reactions)
    reactions = {reactions};
end

if iscell(reactions)
    [~,reactions] = ismember(reactions,model.rxns);
end

if isfield(model, 'subSystems')
    present = cellfun(@(x) any(ismember(x,subSystem)),model.subSystems(reactions));
else
    present = false(numel(reactions));
end

author: opencobra

repo: cobratoolbox


text:
function imgOut = ConvertXYZtoIPT(img, inverse)
%
%       imgOut = ConvertXYZtoIPT(img, inverse)
%
%
%        Input:
%           -img: image to convert from XYZ to IPT or from IPT to XYZ.
%           -inverse: takes as values 0 or 1. If it is set to 1 the
%                     transformation from XYZ to IPT is applied, otherwise
%                     the transformation from IPT to XYZ.
%
%        Output:
%           -imgOut: converted image in XYZ or IPT.
%
%     Copyright (C) 2011  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(inverse == 0)
    imgLMS = ConvertXYZtoLMS(img, 0);

    imgOut = ConvertLMStoIPT(imgLMS, 0);

else
    imgLMS = ConvertLMStoIPT(img, 1);

    imgOut = ConvertXYZtoLMS(imgLMS, 1);
end

end

author: banterle


text:
function ft = f13_title ( )

%*****************************************************************************80
%
%% F13_TITLE returns the title for function f13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string FT, the function title.
%
  ft = 'Rescaled Runge function';

  return
end


text:
function c = bitshift(a,b)
	checkinputs(a,b);

	if numel(b)==1
		if b<0
			c = int64rightshift(int64(a),int64(-b));
		else
			c = int64leftshift(int64(a),int64(b));
		end
	else
		c = zeros(size(b),'int64');
		neg = b < 0;
		if numel(a)==1
			c(neg)  = int64rightshift(int64(a),int64(-b(neg)));
			c(~neg) = int64leftshift(int64(a),int64(b(~neg)));
		else
			c(neg)  = int64rightshift(int64(a(neg)),int64(-b(neg)));
			c(~neg) = int64leftshift(int64(a(~neg)),int64(b(~neg)));
		end
	end
end

author: Sable

repo: mcbench-benchmarks


function gradu = singletrack_gradu_f(in1,in2,in3)
%SINGLETRACK_GRADU_F
%    GRADU = SINGLETRACK_GRADU_F(IN1,IN2,IN3)

%    This function was generated by the Symbolic Math Toolbox version 8.3.
%    09-Jan-2020 11:52:34

I_z = in3(2,:);
M = in3(1,:);
T = in2(2,:);
V_vx = in1(4,:);
V_vy = in1(5,:);
c_f = in3(8,:);
delta = in2(1,:);
deltamax = in3(5,:);
l_f = in3(3,:);
maxbrakeWForce = in3(6,:);
maxmotorWForce = in3(7,:);
psi_dot = in1(6,:);
t2 = 1.0./I_z;
t3 = 1.0./M;
t4 = T.*5.0e+1;
t5 = delta.*5.0e+1;
t6 = deltamax.*5.0e+1;
t8 = V_vy.*1i;
t9 = V_vx.*1.0e+2;
t13 = l_f.*psi_dot.*1i;
t7 = -t4;
t10 = -t5;
t11 = -t6;
t12 = tanh(t9);
t24 = V_vx+t8+t13;
t14 = t7+5.0e+1;
t15 = t7-5.0e+1;
t18 = t6+t10;
t21 = t10+t11;
t26 = angle(t24);
t16 = exp(t14);
t17 = exp(t15);
t19 = exp(t18);
t22 = exp(t21);
t20 = t16+1.0;
t23 = t17+1.0;
t25 = t19+1.0;
t29 = t22+1.0;
t27 = 1.0./t20;
t30 = 1.0./t23;
t32 = 1.0./t25;
t34 = 1.0./t29;
t28 = t27.^2;
t31 = t30.^2;
t33 = t32.^2;
t35 = t34.^2;
t36 = t27.*5.0e+1;
t37 = t27-1.0;
t38 = t30.*5.0e+1;
t40 = deltamax.*t32;
t41 = t32-1.0;
t43 = t34-1.0;
t39 = -t36;
t42 = -t38;
t44 = -t40;
t45 = deltamax.*t43;
t47 = t16.*t28.*5.0e+1;
t49 = t17.*t31.*5.0e+1;
t50 = t16.*t28.*2.5e+3;
t53 = t17.*t31.*2.5e+3;
t55 = t6.*t19.*t33;
t56 = deltamax.*t19.*t33.*-5.0e+1;
t57 = t6.*t22.*t35;
t58 = t30.*t37;
t59 = deltamax.*t22.*t35.*-5.0e+1;
t61 = t37.*t38;
t64 = t34.*t41;
t67 = t4.*t16.*t28.*t30;
t69 = t4.*t17.*t31.*t37;
t70 = t5.*t19.*t33.*t34;
t72 = t5.*t22.*t35.*t41;
t46 = -t45;
t48 = -t47;
t51 = -t49;
t52 = -t50;
t54 = -t53;
t60 = T.*t58;
t63 = t4.*t58;
t65 = delta.*t64;
t68 = T.*t30.*t50;
t71 = T.*t37.*t53;
t90 = t56+t59+t64+t70+t72;
t62 = -t60;
t66 = -t65;
t74 = t39+t42+t63+5.0e+1;
t86 = t26+t44+t46+t65;
t87 = t48+t51+t58+t67+t69;
t89 = t52+t54+t61+t68+t71;
t73 = t30+t37+t62;
t75 = exp(t74);
t77 = t40+t45+t66;
t76 = t75+1.0;
t78 = sin(t77);
t81 = cos(t77);
t79 = 1.0./t76;
t80 = t79.^2;
t82 = maxmotorWForce.*t79;
t83 = t79-1.0;
t84 = maxbrakeWForce.*t12.*t83;
t91 = maxmotorWForce.*t75.*t80.*t89;
t92 = maxbrakeWForce.*t12.*t75.*t80.*t89;
t85 = -t84;
t93 = -t92;
t88 = t82+t85;
t94 = t91+t93;
gradu = reshape([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t3.*(c_f.*t78.*t90-c_f.*t81.*t86.*t90+(t73.*t78.*t88.*t90)./2.0),-t3.*((t73.*t94)./2.0+(t87.*t88)./2.0+(t73.*t81.*t94)./2.0+(t81.*t87.*t88)./2.0),0.0,-t3.*(c_f.*t81.*t90+c_f.*t78.*t86.*t90+(t73.*t81.*t88.*t90)./2.0),-t3.*((t73.*t78.*t94)./2.0+(t78.*t87.*t88)./2.0),0.0,-t2.*(c_f.*l_f.*t81.*t90+c_f.*l_f.*t78.*t86.*t90+(l_f.*t73.*t81.*t88.*t90)./2.0),-t2.*((l_f.*t73.*t78.*t94)./2.0+(l_f.*t78.*t87.*t88)./2.0),0.0,0.0,0.0,1.0],[3,7]);


text:
function [ out, component ] = ScoreVideoToComponentVideo( scoreVideo )
%  ScoreVideoToComponentVideo - use component analysis to only retain
%  largest component
%--------------------------------------------------------------------------
%   Params: scoreVideo - the score video
%
%   Returns: out - component video only containing biggest component
%            component - number for max component or 0 if no component was
%            large enough
%--------------------------------------------------------------------------

out = 0;
[L,num] = bwlabeln(scoreVideo);
max = 0;
component = 0;

for i = 1:num
    temp = sum(sum(sum(L==i)));
    if (temp>max)
        component=i;
        max = temp;
    end
end
if component == 0
    out = scoreVideo;
else
    %check if a good portion of biggest component is classified correct.
    %If not then don't show it on screen
    %Good values are 0.007, .006, .0065, .005, .004
    percentage = 0.006 * ( size(scoreVideo,1) * size(scoreVideo,2) * size(scoreVideo,3) );
    %display(max);
    %display(percentage);
    if (max < percentage)
        component = 0;
    else
        out = (L==component);
    end
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function F = spm_Welcome
% Open SPM's welcome splash screen
% FORMAT F = spm_Welcome
% F        - welcome figure handle
%__________________________________________________________________________
% Copyright (C) 2014-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_Welcome.m 7571 2019-04-10 15:23:16Z guillaume $

%-Open startup window, set window defaults
%--------------------------------------------------------------------------
RectW = spm('WinSize','W',1); % 500x280 figure
Rect0 = spm('WinSize','0',1);
Pos   = [Rect0(1)+(Rect0(3)-RectW(3))/2,...
         Rect0(2)+(Rect0(4)-RectW(4))/2,...
         RectW(3),...
         RectW(4)];
WS    = [1 1 1 1]; % spm('WinScale')

PF = spm_platform('fonts');
PF.helvetica = 'helvetica';

F = figure('IntegerHandle','off',...
    'Name',spm('Version'),...
    'NumberTitle','off',...
    'Tag','Welcome',...
    'Units','pixels',...
    'Position',Pos,...
    'Resize','on',...
    'Pointer','arrow',...
    'Color',[1 1 1]*.8,...
    'MenuBar','none',...
    'DefaultUicontrolFontName',PF.helvetica,...
    'Visible','off');

%-Text
%--------------------------------------------------------------------------
hA = axes('Parent',F,...
    'Units','pixels',...
    'Position',[10 10 80 260].*WS,...
    'Visible','Off');

text(0.5,0.5,'SPM',...
    'Parent',hA,...
    'FontName',PF.times,'FontSize',72,...
    'FontAngle','Italic','FontWeight','Bold',...
    'Rotation',90,...
    'VerticalAlignment','Middle',...
    'HorizontalAlignment','Center',...
    'Color',[1 1 1]*.6);

uicontrol(F,'Style','Text',...
    'String','Statistical Parametric Mapping',...
    'Position',[100 245 390 030].*WS,...
    'FontName',PF.helvetica,'FontSize',18,'FontAngle','Italic',...
    'FontWeight','Bold',...
    'ForegroundColor',[1 1 1]*.5,...
    'BackgroundColor',[1 1 1]*.8);

uipanel('Parent',F,'Title','',...
    'Units','pixels',...
    'Position',[100 130 390 110].*WS,...
    'BackgroundColor',[1 1 1]*.75);

uipanel('Parent',F,'Title','',...
    'Units','pixels',...
    'Position',[100 030 390 085].*WS,...
    'BackgroundColor',[1 1 1]*.75);

uicontrol(F,'Style','Text',...
    'String',spm('Ver'),...
    'Position',[105 200 380 030].*WS,...
    'FontName',PF.helvetica,'FontSize',22,'FontWeight','Bold',...
    'ForegroundColor',[0 0.5 0],...
    'BackgroundColor',[1 1 1]*.75);

uicontrol(F,'Style','Text',...
    'Position',[105 175 380 020].*WS,...
    'String','developed by members and collaborators of the',...
    'FontName',PF.helvetica,'FontSize',10,'FontAngle','Italic',...
    'BackgroundColor',[1 1 1]*.75);

uicontrol(F,'Style','Text',...
    'Position',[105 157 380 020].*WS,...
    'String','Wellcome Centre for Human Neuroimaging',...
    'FontName',PF.helvetica,'FontSize',12,'FontWeight','Bold',...
    'BackgroundColor',[1 1 1]*.75);

uicontrol(F,'Style','Text',...
    'Position',[105 137 380 020].*WS,...
    'String','Institute of Neurology, University College London',...
    'FontName',PF.helvetica,'FontSize',11,...
    'BackgroundColor',[1 1 1]*.75);

uicontrol(F,'Style','Text',...
    'String',['Copyright (c) 1991,1994-' datestr(now,'yyyy')],...
    'Position',[100 005 390 020].*WS,...
    'ForegroundColor',[1 1 1]*.5,...
    'BackgroundColor',[1 1 1]*.8,...
    'FontName',PF.helvetica,'FontSize',11,...
    'HorizontalAlignment','center',...
    'BackgroundColor',[1 1 1]*.8);

%-Objects with Callbacks - PET, M/EEG, fMRI, About SPM, SPMweb, Quit
%--------------------------------------------------------------------------
set(F,'DefaultUicontrolFontSize',12,'DefaultUicontrolInterruptible','on');

uicontrol(F,'Style','pushbutton',...
    'String','PET & VBM',...
    'Position',[135 075 100 030].*WS,...
    'CallBack','delete(gcbf);spm(''Clean'');spm(''PET'')',...
    'FontName',PF.helvetica,...
    'ForegroundColor',[0.5 0 0]);

uicontrol(F,'Style','pushbutton',...
    'String','M/EEG',...
    'Position',[245 075 100 030].*WS,...
    'CallBack','delete(gcbf);spm(''Clean'');spm(''EEG'')',...
    'FontName',PF.helvetica,...
    'ForegroundColor',[0.5 0 0]);

uicontrol(F,'Style','pushbutton',...
    'String','fMRI',...
    'Position',[355 075 100 030].*WS,...
    'CallBack','delete(gcbf);spm(''Clean'');spm(''FMRI'')',...
    'FontName',PF.helvetica,...
    'ForegroundColor',[0.5 0 0]);

uicontrol(F,'Style','pushbutton',...
    'String','About SPM',...
    'Position',[135 040 100 030].*WS,...
    'CallBack','spm_help;',...
    'FontName',PF.helvetica,...
    'ForegroundColor',[0 0.5 0]);

uicontrol(F,'Style','pushbutton',...
    'String','SPMweb',...
    'Position',[245 040 100 030].*WS,...
    'CallBack','web(''https://www.fil.ion.ucl.ac.uk/spm/'')',...
    'FontName',PF.helvetica,...
    'ForegroundColor',[0 0.5 0]);

uicontrol(F,'Style','pushbutton',...
    'String','Quit',...
    'FontName',PF.helvetica,...
    'Position',[355 040 100 030].*WS,...
    'CallBack','delete(gcbf)',...
    'Interruptible','off',...
    'ForegroundColor',[1 0 0]);

if ~ismac
    set(findobj(F,'Style','pushbutton'),'BackgroundColor',[1 1 1]*0.75);
end

%-
%--------------------------------------------------------------------------
set(findobj(F),'Units','normalized');
set(findobj(F,'FontUnits','points'),'FontUnits','normalized');
set(F,'Visible','on');
spm_figure('Focus',F);

author: spm


text:
function [nxx,nyy] =  zoom_size(nx,ny,factor)
nxx = nx*factor + 0.5;
nyy = ny*factor +0.5;
end


function [boundingbox] = bbreg(boundingbox,reg)
	%calibrate bouding boxes
	if size(reg,2)==1
		reg=reshape(reg,[size(reg,3) size(reg,4)])';
	end
	w=[boundingbox(:,3)-boundingbox(:,1)]+1;
	h=[boundingbox(:,4)-boundingbox(:,2)]+1;
	boundingbox(:,1:4)=[boundingbox(:,1)+reg(:,1).*w boundingbox(:,2)+reg(:,2).*h boundingbox(:,3)+reg(:,3).*w boundingbox(:,4)+reg(:,4).*h];
end


text:
function [x,f,exitCode]=convexQuadProg(G,a,C,b,numEqConst,epsVal,maxIter)
%%CONVEXQUADPROG Perform quadratic programming on a convex problem with
%                equality and inequality constraints. Specifically, this
%                algorithm solves the optimization problem:
%                minimize_x a'*x+(1/2)*x'*G*x
%                such that C(:,1:numEqConst)'*x=b(1:numEqConst)
%                     and C(:,(numEqConst+1):end)'*x>=b((numEqConst+1):end)
%                A dual active set algorithm for strictly convex problems
%                is used. Strict convexity means that G is positive
%                definite. The algorithm is robust to poorly conditioned
%                matrices. Thus, to handle the semidefinite case,
%                eigenvalue thresholding is applied prior to running the
%                algorithm.
%
%INPUTS: G An nXn real, positive (semi)definite symmetric matrix (n>0).
%        a An nX1 real vector.
%        C An nXm real matrix (m>=0). The first numEqConst constraints (the
%          equality constraints) must be linearly independent or the
%          algorithm will identify the problem as infeasible.
%        b An mX1 vector.
% numEqConst The number of constraints in C that are equality constraints
%          (numEqConst<=m). If this parameter is omitted or an empty matrix
%          is passed, the default numEqConst=0 is used.
%   epsVal A small value that is used to determine whether there is no step
%          in the primal space such that a constraint becomes feasible,
%          meaning that a step must be taken in the dual space and a
%          constraint dropped. This value should be >0 due to finite
%          precision errors. If this parameter is omitted or an empty
%          matrix is passed, the default value of eps() is used.
%  maxIter The maximum number of iterations to allow. The default value if
%          this parameter is omitted or an empty matrix is passed is 200.
%
%OUTPUTS: x The nX1 solution to the optimization problem or an empty matrix
%           if the problem is infeasible.
%         f The value of the objective function at x, or an empty ,matrix
%           if the problem is infeasible.
%  exitCode A value indicating how the algorithm terminated. Possible
%           values are:
%           0 The algorithm found the optimal value x.
%           1 The maximum number of iterations elapsed.
%          -1 The optimization problem is infeasible.
%
%The algorithm of [1] is implemented. However, the algorithmic description
%in [1] omits a few minor details regarding deleting elements from uPlus
%and updating s removing a constraint. The changes are commented in the
%code. Also, the algorithm of [1] is only described for inequality
%constraints. The change to allow equality constraints is performed by
%forcing the first numEqConst columns of C to be chosen to be added to the
%active set and then making sure that they are never removed from the
%active set. If an attempt is made to remove them from the active set, then
%the corresponding dual variables are set to zero rather than removed.
%
%Another change from the algorithm of [1] deals with the qr formulation of
%certain terms. Though the qr formulation of the H and NStar matrices is
%used, the single column deletions and insertions are not performed. Though
%slower than the algorithm in the paper, one does not need to worry about
%cumulative finite precision errors due to many incremental updates.
%
%Also, algorithm [1] did not include an epsVal term to try to deal with
%finite precision effects limiting the ability to determine whether no
%primal step can activate a certain constraint.
%
%EXAMPLE 1:
%This is the example worked out step-by-step in [1]:
% G=[4, -2;-2, 4];
% b=[0;0;2];
% a=[6;0];
% C=[1,0,1;
%    0,1,1];
% [x,f,exitCode]=convexQuadProg(G,a,C,b)
%The optimal solution can be seen to be x=[0.5;1.5] and f=6.5;
%
%EXAMPLE 2:
%In [1], the algorithm was tested on many random problems. Below is an
%implementation of the method described in the paper for n=81 dimensions
%and m=3*n=243 constraints. The first qStar=m/3=81 constraints end up being
%enforced.
% n=81;
% m=3*n;
% qStar=m/3;
%
% G=rand(n,n);
% G=G+G';
% %G is not symmetric and all elements are in the range of -1 to 1
% SCur=sum(G(1,2:end));
% G(1,1)=1+SCur+rand(1);
% for curRow=2:n
%    SPrev=SCur;
%    SCur=sum(G(curRow,[1:(curRow-1),(curRow+1):end]));
%    G(curRow,curRow)=G(curRow-1,curRow-1)+SCur+SPrev+rand();
% end
%
% xStar=-0.5+rand(n,1);
% C=-1+2*rand(n,m);
% C=bsxfun(@times,C,1./sqrt(sum(C.*C,1)));
% s=zeros(m,1);
% s((qStar+1):end)=rand(m-qStar,1);
% b=-(s-C'*xStar);%Sign changed from description in paper to be feasible.
% u=30*rand(m,1);
% u((qStar+1):end)=0;
% a=C*u-G*xStar;
% [x,f,exitCode]=convexQuadProg(G,a,C,b)
%
%EXAMPLE 3:
%This is an example of a problem having both inequality and equality
%constraints:
% G=[4, -2;-2, 4];
% b=[2;1.6];
% a=[6;0];
% C=[1,0;
%    1,1];
% numEqConst=1;%One equality constraint.
% [x,f,exitCode]=convexQuadProg(G,a,C,b,numEqConst)
%The optimal solution is x=[0.4;1.6] and f=6.56.
%
%REFERENCES:
%[1] D.Goldfarb and A. Idnani, "A numerically stable dual method for
%    solving strictly convex quadratic programs," Mathematical Programming,
%    vol. 27, no. 1, pp. 1-33, Sep. 1983.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(numEqConst))
   numEqConst=0;
end

if(numEqConst>size(C,1))
   error('numEqConst is larger than the number of constraints in C')
end

if(nargin<6||isempty(epsVal))
    epsVal=eps();
end

if(nargin<7||isempty(maxIter))
    maxIter=200;
end

%This portion is to compensate for the possibility that G is positive
%semidefinite.
[V,D]=eig(G);
d=diag(D);
dNew=max(d,2*max(size(G))*eps(max(d)));
if(~all(dNew==d))
    %Only change G if some eigenvalues were too small.
    D=diag(dNew);
    G=V*D/V;
    G=(G+G')/2;%Ensure symmetry.
end

%If there are no constraints, the return the unconstrained solution.
if(isempty(b))
    x=-G\a;
    f=(1/2)*a'*x;
    exitCode=0;
    return
end

%Step 0: Find the unconstrained minimum and initialize parameters.
H=inv(G);
x=-G\a;
f=(1/2)*a'*x;
%Values are added to and deleted from A. Though this is inefficient, it is
%questionable whether any particular data structure could make it notably
%more efficient than what Matlab does without reprogramming this function
%in C.
A=[];
q=0;
u=[];
NStar=[];
L=chol(G,'lower');

skipStep1=false;
for curIter=1:maxIter
    %Step 1: Choose a violated constant, if any. We shall choose the most
    %violated constant.
    if(skipStep1==false)
        s=C'*x-b;
        %Force the equality constraints to be added first.
        if(curIter<=numEqConst)
            p=curIter;
        else
            [~,p]=min(s);

            if(any(p==A)||s(p)>=0)
                %If no constraints are violated, the current solution is
                %feasible and optimal.
                exitCode=0;
                return;
            end
        end

        uPlus=[u;0];%Dual variables with one for constraint p added.
        nPlus=C(:,p);
    end

    %Step 2: Check for feasibility and determine a new S-pair
    %Step 2a) Determine the step direction.
    z=H*nPlus;

    %r is empty when no constraints have been added.
    if(q==0)
        r=[];
    else
        r=NStar*nPlus;
    end

    %Step 2b) Determine the step length
    %Step 2bi) Find t1, the maximum step without violating dual
    %          feasibility.
    t1=Inf;
    l=[];
    for i=1:q
        %The added requirement that uPlus(i)>0 arises due to the equality
        %constraints. Otherwise, all of the dual variables enforced would
        %be greater than zero due to the nature of the dual problem.
        if(r(i)>0&&uPlus(i)>0)
            t1Cur=uPlus(i)/r(i);
            if(t1Cur<t1)
                t1=t1Cur;
                l=i;
            end
        end
    end

    %Step 2bii: Find the minimum step in the primal space such that the
    %pth constraint becomes feasible.
    zProd=z'*nPlus;

    %This comparison is not described in [1], but appears to be necessary
    %due to finite precision errors in some instances.
    if(zProd<epsVal)
        t2=Inf;
    else
        t2=-s(p)/zProd;
    end
    t=min(t1,t2);

    %Step 2c: Determine a new S-pair and take the step.
    %Step 2ci: If no step exists in the primal or dual spaces.
    if(~isfinite(t)||(t==0&&~all(uPlus==0)))
        %The problem is infeasible. The t=0 condition comes from adding
        %equality constraints to inequality constraints.
        x=[];
        f=[];
        exitCode=-1;
        return;
    end

    if(~isfinite(t2))
        %Step 2cii: If a step exists only in the dual space.
        uPlus=uPlus+t*[-r;1];
        if(l>numEqConst)%Do not delete equality constraints
            %Delete constraint A(l).
            A(l)=[];
            %Delete the dual variables associated with the constraint.
            %This deletion is not mentioned in the steps given in the paper,
            %but is necessary.
            uPlus(l)=[];
            q=q-1;

            %x did not change, so s(p) does not need to be recomputed.

            [NStar,H]=updateNStarH(q,A,C,L);
        else
            uPlus(l)=0;
        end
        skipStep1=true;
    else
        %Step 2ciii: Take a step in primal and dual space.
        x=x+t*z;
        f=f+t*z'*nPlus*((1/2)*t+uPlus(q+1));
        uPlus=uPlus+t*[-r;1];

        if(t==t2)%If a full step is taken, add constraint p
            u=uPlus;
            q=q+1;
            A(q)=p;

            [NStar,H]=updateNStarH(q,A,C,L);
            skipStep1=false;
        else%If a partial step is taken, then drop constraint k=A(l)
            if(l>numEqConst)%Do not delete equality constraints
                A(l)=[];
                %Delete the dual variables associated with the constraint.
                %This deletion is not mentioned in the steps given in the
                %paper, but is necessary.
                uPlus(l)=[];
                q=q-1;

                [NStar,H]=updateNStarH(q,A,C,L);
            else
                uPlus(l)=0;
            end
            %s must be recomputed for a partial step. This is not mentioned
            %in the paper.
            s=C'*x-b;
            skipStep1=true;
        end
    end
end

%Maximum number of iterations exceeded.
exitCode=1;
end

function [NStar,H]=updateNStarH(q,A,C,L)
%UPDATENSTARH This function is based on the discussion in Section 4 of [1].
%             It uses a qr decomposition, but does not use any of the qr
%             updating methods.
    N=C(:,A);

    opts.LT=true;
    opts.UT=false;
    %Solve B=L\N;
    B=linsolve(L,N,opts);
    [Q,R]=qr(B);
    R=R(1:q,1:q);

    opts.LT=false;
    opts.UT=true;
    %Solve J=(L')\Q;
    J=linsolve(L',Q,opts);
    J1=J(:,1:q);
    J2=J(:,(q+1):end);

    H=J2*J2';
    NStar=R\(J1');
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:
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http ://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Bilateral Guided Upsampling fits an affine bilateral grid gamma for an
%   operator f: input_image --> output_image.
%
% function [gamma, A, b, lambda_spatial, intensity_options] =
%   bguFit(input_image, edge_image, ...
%   output_image, output_weight, grid_size, lambda_spatial, intensity_options)
%
% Inputs:
%
% input_image is a double tensor: the input to the operator f.
%   Dimensions: height x width x num_input_channels
% edge_image is a double matrix: the "edges" we try to respect. For
%   example, the input luminance (try rgb2luminance or rgb2gray).
%   Dimensions: height x width (x 1)
% output_image is a double tensor: the output f(input_image).
%   Dimensions: height x width x num_output_channels
%
% [Optional] weight_image is a double tensor: whether f is defined at each
%   pixel. Weights need only be non-negative (>= 0) and can be any double.
%   Dimensions: height x width x num_output_channels
%   (same size as output_image)
%
% [Optional] grid_size is:
%   [height width depth num_output_channels num_input_channels]
%   If not specified or empty, defaults to:
%   [round(input_height / 16), round(input_width / 16), 8, ...
%    output_channels, input_channels + 1]
%
% [Optional] lambda_spatial controls spatial smoothness, and defaults to 1.
%   lambda_spatial must be positive.
%
% [Optional] intensity_options is a struct:
%   .type:
%     The type of constraint you want on the grid in the intensity
%     direction. Must be one of:
%       'none': no constraints
%       'first': constrain the first derivative dg/dz
%       'second': constrain the second derivative d2g/dz2
%     Default: 'second'
%   .value:
%     The value towards which the first or second derivative should be.
%     Note that this value is in *absolute* output units and is independent
%     of the grid spacing. I.e., if you want your derivatives to be close
%     to 1 in each bin, then set type to 'first' and value to 1, not
%     1/grid_size(3).
%       type = 'none': ignored.
%       type = 'first' or 'second': any double.
%     Default: 0
%   .lambda:
%     The strength of the constraint relative to the other terms.
%     Default: 4e-6 for first derivative, 4e-7 for second derivative.
%
% Outputs:
%
% A, b are the optional output sparse matrix and right hand side vector
% used to solve for gamma. gamma = reshape(A \ b, grid_size).
%
% lambda_spatial and intensity_options are echoed back as optional outputs
% so the test harness can ask for default parameters then retrieve them.
function [gamma, A, b, lambda_spatial, intensity_options] = ...
    bguFit(input_image, edge_image, ...
    output_image, output_weight, grid_size, lambda_spatial, intensity_options)

DEFAULT_LAMBDA_SPATIAL = 1;

DEFAULT_FIRST_DERIVATIVE_LAMBDA_Z = 4e-6; % about 0.01 for default bin sizes.
%DEFAULT_FIRST_DERIVATIVE_LAMBDA_Z = 4e-7; % about 0.001 for default bin sizes.

DEFAULT_SECOND_DERIVATIVE_LAMBDA_Z = 4e-7; % about 0.01 for default bin sizes.
%DEFAULT_SECOND_DERIVATIVE_LAMBDA_Z = 4e-8; % about 0.001 for default bin sizes.

if ~isa(input_image, 'double')
    error('input_image must be double.');
end

if ~isa(output_image, 'double')
    error('model_image must be double.');
end

if ~exist('edge_image', 'var') || isempty(edge_image) || ...
        ~ismatrix(edge_image) || ~isa(edge_image, 'double')
    error('edge_image must be a double matrix (one channel).');
end

if isempty(output_weight)
    output_weight = ones(size(output_image));
end

if ~isa(output_weight, 'double')
    error('weight must be double matrix');
end

if ndims(input_image) < 2
    error('input_image must be at least two-dimensional.');
end

if ndims(output_image) < 2
    error('output_image must be at least two-dimensional.');
end

if ~isequal(size(input_image, 1), size(output_image, 1)) || ...
    ~isequal(size(input_image, 2), size(output_image, 2))
    error('input_image and output_image must have the same width and height.');
end

if ~exist('grid_size', 'var') || isempty(grid_size)
    grid_size = getDefaultAffineGridSize(input_image, output_image);
end

if ~exist('lambda_spatial', 'var') || isempty(lambda_spatial)
    lambda_spatial = DEFAULT_LAMBDA_SPATIAL;
end

if lambda_spatial <= 0
    error('lambda_spatial must be positive.');
end

% If you pass in nothing, default to second derivative.
if ~exist('intensity_options', 'var') || isempty(intensity_options)
    intensity_options.type = 'second';
    intensity_options.value = 0;
    intensity_options.lambda = DEFAULT_SECOND_DERIVATIVE_LAMBDA_Z;
    intensity_options.enforce_monotonic = false;
end

% If you ask for a constraint but are missing some of the parameters.
if strcmp(intensity_options.type, 'first')
    if ~isfield(intensity_options, 'lambda')
        intensity_options.lambda = DEFAULT_FIRST_DERIVATIVE_LAMBDA_Z;
    end

    if ~isfield(intensity_options, 'value')
        intensity_options.value = 0;
    end

    if ~isfield(intensity_options, 'enforce_monotonic')
        intensity_options.enforce_monotonic = false;
    end
elseif strcmp(intensity_options.type, 'second')
    if ~isfield(intensity_options, 'lambda')

        intensity_options.lambda = DEFAULT_SECOND_DERIVATIVE_LAMBDA_Z;
    end

    if ~isfield(intensity_options, 'value')
        intensity_options.value = 0;
    end

    if ~isfield(intensity_options, 'enforce_monotonic')
        intensity_options.enforce_monotonic = false;
    end
else
    if ~isfield(intensity_options, 'enforce_monotonic')
        intensity_options.enforce_monotonic = false;
    end
end

input_height = size(input_image, 1);
input_width = size(input_image, 2);
grid_height = grid_size(1);
grid_width = grid_size(2);
grid_depth = grid_size(3);
affine_output_size = grid_size(4);
affine_input_size = grid_size(5);

% Size of each grid cell in pixels (# pixels per bin).
bin_size_x = input_width / grid_width;
bin_size_y = input_height / grid_height;
bin_size_z = 1 / grid_depth;

num_deriv_y_rows = (grid_height - 1) * grid_width * grid_depth ...
    * affine_output_size * affine_input_size;
num_deriv_x_rows = grid_height * (grid_width - 1) * grid_depth ...
    * affine_output_size * affine_input_size;

% Set up data term Ax = b.
%
% x is the bilateral grid. It is vectorized as a column vector of size n,
% where n = grid_height * grid_width * grid_depth * ...
%   affine_output_size * affine_input_size.
%
% The right hand side b is the *output* image, vectorized as output_image(:).
% I.e., it's a m x 1 column vector:
% [ red0 red1 ... redN green0 green1 ... greenN blue0 blue1 ... blueN ]'
%
% A is an m x n (sparse) matrix.
% It slices into the bilateral grid with linear interpolation at edge_image
% to get an affine model, then applies it.

% TODO: vectorize
% Build slice matrices for each (i,j) entry of the affine model.
weight_matrices = cell(affine_output_size, affine_input_size);
slice_matrices = cell(affine_output_size, affine_input_size);
for j = 1:affine_input_size
    for i = 1:affine_output_size
        %fprintf('Building weight and slice matrices, i = %d, j = %d\n', i, j);
        [weight_matrices{i, j}, slice_matrices{i, j}] = ...
            buildAffineSliceMatrix(input_image, edge_image, grid_size, i, j);
    end
end

% Concat them together.
slice_matrix = sparse(0, 0);
weight_matrix = sparse(0, 0);
for j = 1:affine_input_size
    for i = 1:affine_output_size
        %fprintf('Concatenating affine slice matrices, i = %d, j = %d\n', i, j);
        slice_matrix = [slice_matrix; slice_matrices{i, j}];
        weight_matrix = blkdiag(weight_matrix, weight_matrices{i, j});
    end
end

%fprintf('Building apply affine model matrix\n');
apply_affine_model_matrix = buildApplyAffineModelMatrix(...
    input_image, affine_output_size);

% Complete slicing matrix.
%fprintf('Building full slice matrix\n');
sqrt_w = sqrt(output_weight(:)); % weighted least squares takes sqrt(w).
W_data = spdiags(sqrt_w, 0, numel(output_weight), numel(output_weight));
A_data = W_data * apply_affine_model_matrix * weight_matrix * slice_matrix;
b_data = output_image(:) .* sqrt_w;

% ----- Add deriv_y constraints -----
%fprintf('Building d/dy matrix\n');
A_deriv_y = (bin_size_x * bin_size_z / bin_size_y) * lambda_spatial * ...
    buildDerivYMatrix(grid_size);
b_deriv_y = zeros(num_deriv_y_rows, 1);

% ----- Add deriv_x constraints -----
%fprintf('Building d/dx matrix\n');
A_deriv_x = (bin_size_y * bin_size_z / bin_size_x) * lambda_spatial * ...
    buildDerivXMatrix(grid_size);
b_deriv_x = zeros(num_deriv_x_rows, 1);

% ----- Add intensity constraints -----
%fprintf('Building d/dz matrix\n');
if strcmp(intensity_options.type, 'first')
    A_intensity = (bin_size_x * bin_size_y / bin_size_z) * ...
        intensity_options.lambda * buildDerivZMatrix(grid_size);
    value = intensity_options.lambda * intensity_options.value;
    m = size(A_intensity, 1);
    b_intensity = value * ones(m, 1);
elseif strcmp(intensity_options.type, 'second')
    A_intensity = (bin_size_x * bin_size_y) / (bin_size_z * bin_size_z) * ...
        intensity_options.lambda * ...
        buildSecondDerivZMatrix(grid_size);
    value = intensity_options.lambda * intensity_options.value;
    m = size(A_intensity, 1);
    b_intensity = value * ones(m, 1);
end

% ----- Assemble final A matrix -----
%fprintf('Assembling final sparse system\n');
if ~strcmp(intensity_options.type, 'none')
    A = [A_data; A_deriv_y; A_deriv_x; A_intensity];
    b = [b_data; b_deriv_y; b_deriv_x; b_intensity];
else
    A = [A_data; A_deriv_y; A_deriv_x];
    b = [b_data; b_deriv_y; b_deriv_x];
end

% ----- Solve -----
%fprintf('Solving system\n');
gamma = A \ b;
gamma = reshape(gamma, grid_size);

author: mahmoudnafifi


text:
function [W] = kW2W(kW)
% Convert power from kilowatts to watts.
% Chad A. Greene 2012
W = kW*1000 ;


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% 	Problem 4- Step responses

z1=[];
p1=[-1 -.2+10j  -.2-10j];
k1=10;

H=zpk(z1,p1,k1);

t=0:.1:20
y1=step(H,t)

z2=-3
p2=[-1 -.2+10j  -.2-10j];
k2=10;

[num,den]=zp2tf(z2,p2,k2);

y2=step(num,den,t)

plot(t,y1,t,y2, ':')
legend('Step1', 'Step2')


text:
classdef ConicalSensorState < AbstractSensorState
    %ConicalSensorState Summary of this class goes here
    %   Detailed explanation goes here

    properties
        sensor ConicalSensor

        activeTf(1,1) logical = true;
        steeringModel(1,1) AbstractSensorSteeringModel
        angle(1,1) double {mustBeGreaterThan(angle,0)} = deg2rad(10);
        range(1,1) double {mustBeGreaterThan(range,0)} = 1;
    end

    methods
        function obj = ConicalSensorState(sensor, activeTf, steeringModel, angle, range)
            arguments
                sensor(1,1) ConicalSensor
                activeTf(1,1) logical
                steeringModel(1,1) AbstractSensorSteeringModel
                angle(1,1) double {mustBeGreaterThanOrEqual(angle,0)}
                range(1,1) double {mustBeGreaterThan(range, 0)}
            end

            obj.sensor = sensor;
            obj.activeTf = activeTf;
            obj.steeringModel = steeringModel;
            obj.angle = angle;
            obj.range = range;
        end

        function sensor = getSensor(obj)
            sensor = obj.sensor;
        end

        function activeTf = getSensorActiveState(obj)
            activeTf = obj.activeTf;
        end

        function setSensorActiveState(obj, activeTf)
            obj.activeTf = activeTf;
        end

        function steeringModel = getSensorSteeringMode(obj)
            steeringModel = obj.steeringModel;
        end

        function setSteeringModel(obj, steeringModel)
            obj.steeringModel = steeringModel;
        end

        function angle = getSensorAngle(obj)
            angle = obj.angle;
        end

        function setSensorAngle(obj, angle)
            obj.angle = angle;
        end

        function range = getSensorMaxRange(obj)
            range = obj.range;
        end

        function setSensorMaxRange(obj, range)
            obj.range = range;
        end

        function angle = getMaxAngle(obj)
            angle = obj.getSensorAngle();
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%
for i=1:length(lambda_vec)
	lambda=lambda_vec(i);
	[theta]=trainLinearReg(X, y, lambda);
	lambda=0;	% Set to zero when calculating error since already been trained
	% We are calc J as error not as cost. So, lambda should not be included when
	% calculating error for thetas that have been trained, else it will be biased
	% Refer to 2.1 of Exercise 5, error is computed without lambda
	[e_train]=linearRegCostFunction(X, y, theta, lambda);
	[e_val]=linearRegCostFunction(Xval, yval, theta, lambda);	% J over all CV set for new set of theta

% Accumulating error from i=1:m
	if (i==1)
		error_train=e_train;
		error_val=e_val;
	else
		error_train=[error_train; e_train];
		error_val=[error_val; e_val];
	end
end

% Above methods are also accepted when submitted, following method are also accepted

%for i=1:length(lambda_vec)
%lambda=lambda_vec(i);
%[theta] = trainLinearReg(X, y, lambda);

% For training set
%m=length(y);
%thetas=theta(2:end,1);
%h=X*theta;
%lambda=0;
%error_train(i)=(1/(2.*m))*sum((h-y).^2)+(lambda/(2.*m)).*sum(thetas.^2);

% For CV set
%n=length(yval);
%h2=Xval*theta;
%lambda=0;
%error_val(i)=(1/(2.*n))*sum((h2-yval).^2)+(lambda/(2.*n)).*sum(thetas.^2);
%end

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

end

author: yhyap

repo: machine-learning-coursera


text:
function initializeDateListBox( dates, hdls )
%INITIALIZEDATELISTBOX Initialize date list box:
%   dateListBox: a list of dates in format yyyymmdd for the user to select
%   generateDataButton: a button to generate data
%   forward/backward buttons: buttons to move forward/backward in time
%   Detailed explanation goes here

% Create list box for choosing date to display eddies
date_list = arrayfun(@num2str, dates, 'UniformOutput', false);

dateListBox = uicontrol(hdls.fig, ...
    'Style', 'listbox', ...
    'Units', 'normalized', ...
    'String', date_list, ...
    'Position', [0.875 .6 .1 .4], ...
    'Callback', {@dateListBoxCallback, hdls});

% Forward button
uicontrol(hdls.fig, ...
    'Style', 'pushbutton', ...
    'Units', 'normalized', ...
    'String', 'Forward', ...
    'Callback', {@forward_button_Callback, dateListBox, hdls}, ...
    'position', [0.925 0.33 .05 .05]);

% backward button
uicontrol(hdls.fig, ...
    'Style', 'pushbutton', ...
    'Units', 'normalized', ...
    'String', 'Backward', ...
    'Callback', {@backward_button_Callback, dateListBox, hdls}, ...
    'position', [0.875 0.33 .05 .05]);

hdls.eddyInfoText = uicontrol(hdls.fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'String', 'Eddy Information', ...
    'position', [0.875 0 0.1 .3]);
end

function dateListBoxCallback(hObject, ~, hdls)

old_pointer = start_busy_pointer(hdls);

hdls.currentDateIndex = get(hObject, 'Value');
hdls.generate_data();

end_busy_pointer(hdls, old_pointer);

end

function forward_button_Callback(~, ~, dateListBox, hdls)

old_pointer = start_busy_pointer(hdls);

curr_date_index = get(dateListBox, 'Value');
if(curr_date_index == length(get(dateListBox, 'String')))
    return;
end
set(dateListBox, 'Value', curr_date_index + 1);
hdls.currentDateIndex = hdls.currentDateIndex + 1;
hdls.generate_data();

end_busy_pointer(hdls, old_pointer);

end

function backward_button_Callback(~, ~, dateListBox, hdls)

old_pointer = start_busy_pointer(hdls);

if hdls.currentDateIndex == 1
    return;
end
hdls.currentDateIndex = hdls.currentDateIndex - 1;
set(dateListBox, 'Value', hdls.currentDateIndex);
hdls.generate_data();

end_busy_pointer(hdls, old_pointer);

end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
input_image_path = 'image1.jpg';

im = double(imread(input_image_path));

[masked_image] = separate(im,1);

imshow(uint8(masked_image.*im))


text:
% clc, clear, close all
%
% %% compile
% % mex kdtree_build.cpp
% % mex kdtree_k_nearest_neighbors.cpp
% % disp('necessary files compiled.');
%
% %% test1: small test for evaluation and visualizatoin
% rand('seed',1);
% p = rand( 30,  2 ); % input data
% q = [ .5, .5 ]'; k = 5; % query data (buggy)
%
% % execute query by kd-tree and by linear scan
% tree = kdtree_build( p );
% idxs1 = kdtree_k_nearest_neighbors( tree,q,k );
% dist = zeros( size(p,1), 1 );
% for i=1:size(p,1)
%     dist( i ) = sqrt( sum( ( p(i,:)-q' ).^2 ) );
% end
% [IGNORE, idxs2] = sort( dist,'ascend' );
% idxs2 = idxs2(1:k); % keep only first k entries
%
% %%% visualize indexes returned by query
% disp( ['(fullsearch)   kNN indexes: ', sprintf('%d ', idxs2 ) ] );
% disp( ['(kdtreesearch) kNN indexes: ', sprintf('%d ', idxs1 ) ] );
%
% %%% visualize
% xlim( [0 1] );
% ylim( [0 1] );
% hold on; axis equal; axis off;
% plot( p(:,1), p(:,2), '.b');
% plot(q(1), q(2),'.r');
% plot(p(idxs2,1), p(idxs2,2),'or');
% plot(p(idxs1,1), p(idxs1,2),'+y');
% legend('database', 'query', '(matlab) result', '(kdtree) result');
%
% %%% convert input data to C++ format (for testing)
% % for i=1:size(p,1)
% %    disp(sprintf('A[%d][0] = %.2f; A[%d][1] = %.2f;', i-1, p(i,1), i-1, p(i,2) ));
% % end
%
% return

%% test2: extensive random performance and correctness test
% clc, clear, close all;
% NUMTESTS = 10;
% NUMQUERS = 100;
% N = 10000;
% counter = 1;
% kdtree_times = zeros( NUMTESTS*NUMQUERS,1 );
% matlab_times = zeros( NUMTESTS*NUMQUERS,1 );
%
% for j=1:NUMTESTS
%     disp(sprintf('executing test #%d',j));
%     p = rand( 10000,   2 ); % input data
%     q = rand( NUMQUERS,2 ); % query data
%     k = ceil(N/10*rand( 1,1 ));   % query size (number of kNN to extract)
%
%     %% execute query by kd-tree
%     tree = kdtree_build( p );
%
%     for m=1:NUMQUERS
%         % kdtree based
%         tic
%         idxs_kdtree = kdtree_k_nearest_neighbors( tree,q(m,:)',k );
%         kdtree_times( counter ) = toc();
%
%         % matlab based
%         tic
%         dist = zeros( size(p,1), 1 );
%         for i=1:size(p,1)
%             dist( i ) = sqrt( sum( ( p(i,:)-q(m,:) ).^2 ) );
%         end
%         [IGNORE, idxs2] = sort( dist,'ascend' );
%         idxs_matlab = idxs2(1:k);
%         matlab_times( counter ) = toc();
%
%         if( ~all( idxs_kdtree==idxs_matlab) )
%             error('kdtree gave incorrect results');
%         end
%
%         counter = counter + 1;
%     end
% end
%
% subplot( 211 ), plot( kdtree_times ); ylim( [0,0.05] );
% subplot( 212 ), plot( matlab_times ); ylim( [0,0.05] );
% disp( sprintf('kdtree average time (ms): %d ', mean(kdtree_times)/1000 ) );
% disp( sprintf('matlab average time (ms): %d ', mean(matlab_times)/1000 ) );
% return;

% test 3: check if distances returned are correct!
clc, clear, close all
rand('twister',1);
p = rand( 1000,  2 ); % input data
q = [ .5, .5 ];
k = 100;

% execute query by kd-tree and by linear scan
tree = kdtree_build( p );
[idxs1, dists] = kdtree_k_nearest_neighbors( tree,q,k );
dist = zeros( size(p,1), 1 );
for i=1:size(p,1)
    dist( i ) = sqrt( sum( ( p(i,:)-q ).^2 ) );
end
[IGNORE, idxs2] = sort( dist,'ascend' );
idxs2 = idxs2(1:k); % keep only first k entries
disp( sprintf('search is correct?: %d', all(idxs1==idxs2)));

% compare distances against computed
matdist = zeros( length(idxs1),1 );
for i=1:length(idxs1)
    matdist(i) = sqrt( sum((q-p(idxs1(i),:)).^2) );
end
figure(2);
subplot 211, plot( matdist );
subplot 212, plot( dists );

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function nd = ndim(c)

%tstoolbox/@core/ndim
%   Syntax:
%     * ndim(c)
%
%   Input Arguments:
%     * c - core object
%
%   return number of dimensions, a scalar value has 0 dimensions
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

nd = length(c.dlens);
if nd == 1
	if c.dlens(1) == 1
		nd = 0;
	end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function y=idft(p1,p2,p3,p4)

% y=idft(x,...);
%
% Exactly the same as the ifft function except that the output of this function
% is scaled so that sigpower(x) is the same as sigpower(dft(x));

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% 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/>.

error(nargchk(1,4,nargin));

scale=sqrt(length(p1));
switch nargin
  case 1
    y=ifft(p1)*scale;
  case 2
    y=ifft(p1,p2)*scale;
  case 3
    y=ifft(p1,p2,p3)*scale;
  case 4
    y=ifft(p1,p2,p3,p4)*scale;
  otherwise
    error('Check code...');
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function new_hist = merge_histograms(sp, parts)
% Merges the histograms of several superpixels by weighting the individual
% histograms by the superpixel sizes. This gives the actual histogram of
% the combined superpixel.

partsize = zeros(1, length(parts));
for i = 1:length(parts)
    partsize(i) = sp{parts(i)}.size;
end

for fn = 1:length(sp{1}.hist)
    new_hist{fn} = zeros(length(sp{1}.hist{fn}),1);
    cc = 0;
    for i = 1:length(parts)
        cur_hist = sp{parts(i)}.hist{fn};
        if sum(cur_hist) > 0 % skip empty histograms
            new_hist{fn} = new_hist{fn} + partsize(i)*cur_hist;
            cc = cc + partsize(i);
        end
    end
    if cc > 0
        new_hist{fn} = new_hist{fn}/cc; % normalization
    else
        % If there are only empty histograms, the output histogram remains as a vector of zeros
    end
end

author: Cloud-CV

repo: object-proposals


text:
function M_hat = run_mc(params)
  M = params.M;
  Idx = params.Idx;
  warning('off','all');
  [numr,numc] = size(M);
  rank = 2;
  Known = find(Idx);
  data = M(Known);
  [U,S,V] = svp(Known,data,numr,numc,rank);
  M_hat = U*diag(S)*V';
  warning('on','all');
end

author: andrewssobral

repo: lrslibrary


text:
function laith()
% LAITH Graphical Simulation of an elevator system
% laith() creates a GUI with the number of floors as input. Granting this number
% creates a building with this number of floors and an elevator that moves up
% and down in it according to the will of animated people (modeled as balls).
%
% To run, type laith in the command window or press F5 from this editor window.
%
% Laith Alkurdi is the name of a friend who worked with me on the code.
%
% Husam Aldahiyat
% Feb. 2009
% numandina@gmail.com
%
%
%% Creating Figure and Uicontrols
% create figure
figure('units','normalized','position',[.1 .1 .8 .8],'color',[1 1 1],...
	'menubar','none','numbertitle','off','name','Elevator')

% create axes
axes('position',[.35 .1 .65 .9])

% create pushbutton (Start)
pb=uicontrol('style','pushbutton','units','normalized','position',...
	[.1 .5 .04 .025],'backgroundcolor',[1 1 1],'string','Start',....
	'fontweight','bold','callback',@dogo);

% create edit box
ed2=uicontrol('style','edit','units','normalized','position',[.1 .525 .04 .025],...
	'backgroundcolor',[1 1 1],'fontweight','bold');

% create text (Floors)
uicontrol('style','text','units','normalized','position',[.1 .55 .04 .025],...
	'backgroundcolor',[1 1 1],'string','Floors','fontweight','bold')

% creates debug text
tee=uicontrol('style','text','foregroundcolor',[1 0 0],'units','normalized',...
	'position',[.025 .6 .25 .35],'backgroundcolor',[0 0 0],'max',2,'fonts',15,...
	'fontname','courier');

% this is so the the floor numbers on the z-axis will have dotted lines
% showing their connection to the elevator
grid on

% initially, we can't see the axes
axis off

% this is just so we can move the camera around during play
rotate3d

%% Drawing the Elevator
% when we press the pushbutton we go here
	function dogo(varargin)

		% if we press the 'Start' button, change its string to 'Stop'
		if strcmp(get(pb,'string'),'Start')
			set(pb,'string','Stop')
		else

			% if we press the 'Stop' button, change its string to 'Start' and
			set(pb,'string','Start')

			% clear the axes
			cla

			% and halt operation
			return
		end

		% make axes appear
		axis on

		% clear debugging text
		set(tee,'string','')

		% get number of floors and turn it from string to number
		flrz=str2double(get(ed2,'string'));

		% check input
		if flrz < 2
			errordlg('Number of floors needs to be at least 2','Stupid!')
			set(pb,'string','Start')
			return
		end

		% this is compensated later by adding 1 to flrz
		% it's stupid but I like to do things like this
		flrz=flrz-1;

		% the following (a) and (b) matrices are very important. they create our
		% 'cube' shape.
		% the first column of (a) determines the min and max values for the (x),
		% the second for (y) and third for (z). So initially we have a unit cube
		% because the min and max for (x),(y) and (z) are (0) and (1).
		a= [0 0 0;
			1 0 0;
			1 1 0;
			0 1 0;
			0 0 1;
			1 0 1;
			1 1 1;
			0 1 1];

		% this (b) matrix specifies the faces of said cube. The cube has six faces
		% and eight vertices (points). the first row of (b) specifies which points
		% of (a) lie on the [b] face.
		% For example, the first face is (1,2,6,5), which are the rows of (a).
		% This mean the points [0,0,0],[1,0,0],[1,0,1] and [0,0,1] lie on the
		% first face, and so on.
		b= [1 2 6 5;
			2 3 7 6;
			3 4 8 7;
			4 1 5 8;
			1 2 3 4;
			5 6 7 8];

		% we want to make the elevator shaft first.
		ac=a;

		% take the max for (z) and multiply it by the number of floors. I have
		% decided that each floor is (0.2) high.
		ac(5:8,3)=.2*(flrz+1);

		% make the 3D patch that has height of the number of floors. this will act
		% as the elevator shaft, or building which our elevator will move
		% through in.
		% the facealpha makes the shaft kind of transparent so you would see the
		% people inside
		patch('vertices',ac,'faces',b,'edgecolor','k','facecolor',[.3 .3 .3],...
			'facevertexalphadata',0.5,'facealpha','flat');

		% next we make the elevator itself
		c=a;

		% for this, just make the (z) column of (a) from (0) to (0.2),
		% this means our elevator is placed initially between (0) and (0.2)
		% in the (z) axis
		c(:,3)=c(:,3).*.2;

		% create our elevator!
		elv=patch('vertices',c,'faces',b,'edgecolor','w','facecolor',[0 0 1],...
			'facevertexalphadata',0.3,'facealpha','flat');

		% next we will make the two elevator doors
		c1=a;

		% the (x) of these doors is zero (flat!)
		c1(:,1)=c1(:,1).*0;

		% the (z) is the same as the elevator
		c1(:,3)=[0;0;0;0;.2;.2;.2;.2];

		% for the first door, the (y) goes from (0) to (0.5)
		c1(:,2)=c1(:,2).*.5;

		% create door number one!
		door1 = patch('vertices',c1,'faces',b,'edgecolor','k','facecolor',...
			[0 0 .8],'facevertexalphadata',1,'facealpha','flat','linewidth',2);

		% same thing for our second door
		c2=a;

		% flat (x)
		c2(:,1)=c2(:,1).*0;

		% (z) from (0) to (0.2)
		c2(:,3)=[0;0;0;0;.2;.2;.2;.2];

		% (y) from (0.5) to (1)
		c2(:,2)=c2(:,2).*.5+.5;

		% create second door!
		door2 = patch('vertices',c2,'faces',b,'edgecolor','k','facecolor',[0 0 1],'facevertexalphadata',1,'facealpha','flat',...
			'linewidth',2);

		% next we want to make flat surfaces for each ground as a reference

		% for each floor go through the loop
		for jk=1:flrz+1

			pk=a;

			% make the (z) equal a number multiple of (0.2)
			% since the min and max of (z) are the same, the resulting surface
			% is flat
			pk(:,3)=repmat(.2*(jk-1),8,1);

			% create this floor
			patch('vertices',pk,'faces',b,'edgecolor','k','facecolor','k','facevertexalphadata',0.9,'facealpha','flat');

		end

		% by now we'll have a transparent elevator shaft with floors and an
		% elevator inside. we'll next want to cover one of the shaft sides with
		% a kind of curtain patch right above and below the elevator.
		% to do this, we will create two patches, one right above the elevator
		% and extends to the roof, while the other will have a min of the ground
		% and a max of the lower limits of our elevators.
		% these two patches will have the cyan face colour

		% flat (x)
		ar=[0 0 0;
			0 0 0;
			0 1 0;
			0 1 0;
			0 0 1;
			0 0 1;
			0 1 1;
			0 1 1];

		% max (z) depends on number of floors
		ar(5:8,3)=.2*(flrz+1);

		% for our first patch, min (z) stays at (0), while max (z) changes with
		% the min height of the elevator, so initially, max (z) is (0)

		af1=ar;
		af1(5:8,3)=repmat(0,4,1);
		fp1=patch('vertices',af1,'faces',b,'edgecolor','k','facecolor','c',...
			'facevertexalphadata',1,'facealpha','flat');

		% for our second patch, max (z) stays at the roof, while min (z) changes
		% with elevator max height. initially, min (z) is (0.2)
		af2=ar;
		af2(1:4,3)=repmat(0.2,4,1);
		fp2=patch('vertices',af2,'faces',b,'edgecolor','k','facecolor','c',...
			'facevertexalphadata',1,'facealpha','flat');

		%  next we'll mark the (z) axis with labels indicating the floor names
		FN=cell(flrz+1,1);
		for kkp=1:flrz+1
			FN{kkp}=['Floor ',num2str(kkp)];
		end

		% by now we have a cell (FN) containing strings for the floors numbers,
		% so we then set the (z) tick label as this string, and the (z) ticks
		% as multiples of (0.2)
		set(gca,'zticklabel',FN,'ztick',.1:.2:(flrz+1)*.2)

		% hide (x) and (y) ticks
		set(gca,'xtick',[]);
		set(gca,'ytick',[]);

		% the following are saved points that will be used later (kind of like
		% lookup table)

		% first are the co-ordinates for the points -inside- the elevator
		in={[.25 .1];[.35 .1];[.5 .1];[.65 .1];[.8 .1];[.925 .1]};

		% these are the co-ordinates for the points outside the elevator after
		% getting out from it
		out={[-.2 .925];[-.2 .75];[-.2 .575];[-.2 .425];[-.2 .275];[-.2 .125]};

		% points outside the elevator waiting to get inside it
		toin={[-.2 .1];[-.2 .25];[-.2 .4];[-.2 .55];[-.2 .7];[-.2 .85]};

		% run function (do1)
		do1

		function do1()

			% create the ground (lobby)

			% the groud axis limits
			v(1)=-.25;
			v(2)=1;
			v(3)=0;
			v(4)=1;

			% (z) is (0) and flat, (y) is same as elevator shaft, while (x) min
			% is slightly lower (so the ground would appear sticking out)
			aaa=[v(1) v(3) 0;
				v(2) v(3) 0;
				v(2) v(4) 0;
				v(1) v(4) 0;
				v(1) v(3) 0;
				v(2) v(3) 0;
				v(2) v(4) 0;
				v(1) v(4) 0];
			b=[1 2 6 5;2 3 7 6;3 4 8 7;4 1 5 8;1 2 3 4;5 6 7 8];

			% create the ground!
			patch('vertices',aaa,'faces',b,'edgecolor',[1 1 1],'facecolor',...
				[.3 .3 .3],'facevertexalphadata',0.5,'facealpha','flat');

			% set proper view (3D)
			view(3)

			% current floor is (0)
			cf=0;

			% dummy variable, means nothing
			dumnop=4;

			% choose number of balls randomly from (1) to (6)
			num=floor(6*rand)+1;

			% create these balls (people)
			createballs(num)

			% counter to be used later
			countz=0;

			% choose heading floor randomly between (2) and max number of floors
			tofloor=floor(rand*(flrz))+1;
%
% [pa1,pa2,aa,aa1,aa2,aa3,ha]=pigeon(.05);
% 			function movp(varargin)
%
% 				for kkl=1:length(ha)
% 					set(ha(kkl),'zdata',get(ha(kkl),'zdata')+.5)
% 					set(ha(kkl),'ydata',get(ha(kkl),'ydata')+.5)
% 					set(ha(kkl),'xdata',get(ha(kkl),'xdata')-.5)
% 				end
% 				ppI=get(pa1,'vertices');
% 				ppII=get(pa2,'vertices');
% 				ppI(:,3)=ppI(:,3)+.5;
% 				ppI(:,2)=ppI(:,2)+.5;
% 				ppI(:,1)=ppI(:,1)-.5;
% 				ppII(:,3)=ppII(:,3)+.5;
% 				ppII(:,2)=ppII(:,2)+.5;
% 				ppII(:,1)=ppII(:,1)-.5;
%
% 				aa2(:,1)=aa2(:,1)-.5;
% 				aa2(:,2)=aa2(:,2)+.5;
% 				aa2(:,3)=aa2(:,3)+.5;
% 				aa1(:,1)=aa1(:,1)-.5;
% 				aa1(:,2)=aa1(:,2)+.5;
% 				aa1(:,3)=aa1(:,3)+.5;
% 				aa3(:,1)=aa3(:,1)-.5;
% 				aa3(:,2)=aa3(:,2)+.5;
% 				aa3(:,3)=aa3(:,3)+.5;
% 				aa(:,1)=aa(:,1)-.5;
% 				aa(:,2)=aa(:,2)+.5;
% 				aa(:,3)=aa(:,3)+.5;
%
% 				aa
% 				aa1
% 				aa2
% 				aa3
%
% 				get(pa1,'faces')
% 				for kk=1:2
% 					set(pa1,'vertices',aa2)
% 					set(pa2,'vertices',aa3)
% 					pause(.1)
% 					set(pa1,'vertices',aa)
% 					set(pa2,'vertices',aa1)
% 					pause(.1)
% 				end
%
%
% 			end
%% Main Loop
			% while 'Stop' pushbutton has not yet been pressed
			while strcmp(get(pb,'string'),'Stop')

				% increment counter
				countz=countz+1;

				% choose speed for opening the door
				dspeed=.1;

				% open door
				opend(dspeed)

				% move balls inside
				for j=1:num

					% move each ball from its location (toin) to its heading
					% point (in)
					hs(j)=movb(rz(j),toin{j},in{j},rand*10+15,j,hs(j)); %#ok
				end

				% choose speed for closing the door
				dspeed=.1;

				% close door
				closed(dspeed)

				% choose speed of going up
				upspeed=.035;

				% update debugging text
				set(tee,'string',{['Number of Stops: ',num2str(countz)];...
					sprintf('Current Stop: Floor %d',tofloor+1);...
					['Number of Dudes: ',num2str(num)]})
% 				movp
				% have elevator go up from (cf) to (tofloor) with speed
				% (upspeed) and ball handles (hs)
				cf=goup(upspeed,tofloor,cf,hs);

				% choose speed of door opening again
				dspeed=.075;

				% open door
				opend(dspeed)

				% next we'll want to move each ball outside the elevator
				% also, we want to get their handles

				% preallocation
				todel=zeros(1,num*2);

				% for each ball go once around the loop
				for kmpf=1:num

					% move ball outside its location to (out) and draw virtual
					% ground (patch) and get the patche's handle
					[hs2,pn]=movb(rz(num),dumnop,out{kmpf},20,kmpf,hs(kmpf));

					% store handles for each ball and ground here
					todel(kmpf*2-1:kmpf*2)=[hs2,pn];
				end

				% now it's time to make the balls and patches gradually
				% disappear

				% initial transperancy is (0.5)
				for k=.5:-.075:0

					% gradually lower the transperancy of everything (every
					% handle) in (todel)
					set(todel,'facealpha',k,'edgealpha',k)
					pause(.05)

				end

				% finally after everything is completely transparent, delete it
				delete(todel)

				% choose speed for closing the door
				speedc=.075;

				% close doors
				closed(speedc)

				% choose a new number of dudes from (1) to (6)
				num=floor(6*rand)+1;

				% choose a new floor number from (2) to number of floors
				tofloor=floor(rand*(flrz))+1;

				% update debugger text
				set(tee,'string',{['Number of Stops: ',num2str(countz)];...
					sprintf('Current Stop: Floor %d',tofloor+1);...
					['Number of Dudes: ',num2str(num)]})

				% create the balls (people)
				createballs(num)

				% choose speed to returning the elevator to the ground
				speedret=.035;

				% return elevator to the ground
				cf=retn(speedret,cf);

			end

%% Create Balls
			function createballs(num)

				% choose random radii for each ball (fat, tall, young, etc...)
				% from (0.02) to (0.05)
				rz=rand(num,1).*.03+.02;

				% next we create the balls by moving them from (toin) to (toin)
				hs=zeros(num,1);
				for n2=1:num
					hs(n2)=movb(rz(n2),toin{n2},toin{n2},10,n2);
				end
			end

		end

%% Raising the Elevator
		function cf=goup(inc_speed,floor,cf,hs)

			% counter for use later
			cc=1;
			while true

				% condition for stopping the elevator: next elevtor position
				% (indicated by current position + speed) is higher than
				% required floor, in which case, set the position of the
				% elevator as the next floor minimum coordinate
				if min(c(:,3))+inc_speed>floor*.2

					% min (z) of elevator is same as floor height
					c(1:4,3)=repmat(floor*.2,4,1);

					% max (z) of elevator is same as floor ceiling (next floor)
					c(5:8,3)=repmat(floor*.2+.2,4,1);

					% move elevator to this position
					set(elv,'vertices',c)

					% do the same thing for the elevator doors
					c1(1:4,3)=repmat(floor*.2,4,1);
					c1(5:8,3)=repmat(floor*.2+.2,4,1);
					c2(1:4,3)=repmat(floor*.2,4,1);
					c2(5:8,3)=repmat(floor*.2+.2,4,1);
					set(door1,'vertices',c1)
					set(door2,'vertices',c2)

					% also move the curtains in the same way
					af1=get(fp1,'vertices');
					af1(5:8,3)=repmat(floor*.2,4,1);
					af2=get(fp2,'vertices');
					af2(1:4,3)=repmat(.2+floor*.2,4,1);
					set(fp1,'vertices',af1)
					set(fp2,'vertices',af2)

					% and move the balls as well
					for k=1:length(hs)
						ppv=get(hs(k),'zdata');
						set(hs(k),'zdata',ppv+floor*.2-min(min(ppv)))
					end

					% update plot
					drawnow()

					% leave loop
					break
				end

				% increment counter
				cc=cc+1;

				% increase height of elevator (z) by its speed
				c(:,3)=c(:,3) + inc_speed;
				set(elv,'vertices',c)

				% same with doorsw
				c1(:,3)=c1(:,3) + inc_speed;
				c2(:,3)=c2(:,3) + inc_speed;
				set(door1,'vertices',c1)
				set(door2,'vertices',c2)

				% change the height of curtains in the same way
				af1=get(fp1,'vertices');
				af1(5:8,3)=repmat(0+cc*inc_speed+cf*.2,4,1);
				af2=get(fp2,'vertices');
				af2(1:4,3)=repmat(0.2+cc*inc_speed+cf*.2,4,1);
				set(fp1,'vertices',af1)
				set(fp2,'vertices',af2)

				% increase the height of the balls as well
				for k=1:length(hs)
					set(hs(k),'zdata',get(hs(k),'zdata')+inc_speed)
				end

				% update plot
				drawnow()

			end

			% update floor number (output)
			cf=floor;
		end

%% Returning the Elevator
		function cf=retn(inc_speed,cf)

			% counter
			cc=1;

			% destination floor
			floor=0;

			% the following is the same as raising the elevator, with the
			% replacement of a few addition with subtractions
			while true

				if min(c(:,3))-inc_speed<floor*.2

					c(1:4,3)=repmat(floor*.2,4,1);
					c(5:8,3)=repmat(floor*.2+.2,4,1);
					set(elv,'vertices',c)

					c1(1:4,3)=repmat(floor*.2,4,1);
					c1(5:8,3)=repmat(floor*.2+.2,4,1);
					c2(1:4,3)=repmat(floor*.2,4,1);
					c2(5:8,3)=repmat(floor*.2+.2,4,1);
					set(door1,'vertices',c1)
					set(door2,'vertices',c2)

					af1=get(fp1,'vertices');
					af1(5:8,3)=repmat(floor*.2,4,1);
					af2=get(fp2,'vertices');
					af2(1:4,3)=repmat(.2+floor*.2,4,1);
					set(fp1,'vertices',af1)
					set(fp2,'vertices',af2)
					drawnow()
					break
				end
				cc=cc+1;
				c(:,3)=c(:,3) - inc_speed;
				set(elv,'vertices',c)

				c1(:,3)=c1(:,3) - inc_speed;
				c2(:,3)=c2(:,3) - inc_speed;
				set(door1,'vertices',c1)
				set(door2,'vertices',c2)

				af1=get(fp1,'vertices');
				af1(5:8,3)=repmat(0-cc*inc_speed+cf*.2,4,1);
				af2=get(fp2,'vertices');
				af2(1:4,3)=repmat(0.2-cc*inc_speed+cf*.2,4,1);
				set(fp1,'vertices',af1)
				set(fp2,'vertices',af2)
				drawnow()

			end
			cf=floor;
		end

%% Opening Elevator Doors
		function opend(dspeed)

			% get location information for both doors
			var1=get(door1,'vertices');

			% for use in the while loop that lies ahead
			y1=var1(:,2);
			var2=get(door2,'vertices');

			% loop until next location of first door is lower than its fully
			% opened location (y) = (0), in which case get out of the loop and
			% manually set both door locations
			while max(y1)-dspeed>0;

				% move both doors relative to the speed, in opposite directions
				var1([3,4,7,8],2)=var1([3,4,7,8],2)-dspeed;
				var2([1,2,5,6],2)=var2([1,2,5,6],2)+dspeed;
				set(door1,'vertices',var1)
				set(door2,'vertices',var2)

				% for the loop condition
				y1=var1(:,2);

				% update plot
				drawnow()
			end

			% set both doors to fully open status
			var1(:,2)=zeros(8,1);
			var2(:,2)=ones(8,1);
			set(door1,'vertices',var1)
			set(door2,'vertices',var2)

			% update plot
			drawnow()

		end

%% Closing Elevator Doors
		function closed(dspeed)

			% very similar to opening the elevator doors, with subtraction
			% operation replacing the addition, as well as having the fully
			% closed conditions being one door having a maximum (y) of (0.5)
			% while the other has a minimum (y) of (0.5).
			var1=get(door1,'vertices');
			y1=var1(:,2);
			var2=get(door2,'vertices');
			while max(y1)+dspeed<0.5;
				var1([3,4,7,8],2)=var1([3,4,7,8],2)+dspeed;
				var2([1,2,5,6],2)=var2([1,2,5,6],2)-dspeed;
				y1=var1(:,2);
				set(door1,'vertices',var1)
				set(door2,'vertices',var2)
				drawnow()
			end
			var1([3,4,7,8],2)=repmat(.5,4,1);
			var2([1,2,5,6],2)=repmat(.5,4,1);
			set(door1,'vertices',var1)
			set(door2,'vertices',var2)
			drawnow()

		end

%% Moving the Balls
		function [h,pn]=movb(r,A,B,V,nm,varargin)

			% r: Radius of ball
			%
			% A: Starting position of centre of ball
			%
			% B: Ending position of centre of ball
			%
			% V: Speed of rotation
			%
			% nm: Number denoting the colour
			%
			% varargin: Could contain handle of previously created ball

			% so previous plots would not disappear
			hold on

			% if no previous ball is created (i.e. we want to create instead
			% of move)
			if ~numel(varargin)
				% create unit ball
				[a,b,cC2]=sphere(25);

				% adjust size of ball by multiplying by the radius
				% also, move ball to starting position by adding
				% and put ball on the ground (0z) by adding r to z
				h=surf(a.*r+A(1),b.*r+A(2),cC2.*r+r);

				% make colors cooler
				set(h,'facecolor','interp','edgecolor','interp')

				% get ball colour data
				s2=get(h,'cdata');

				% take the ball top and change its colour to something
				% now each ball has a cute little colourful hat!
				s2(end-3:end,:)=repmat(.15*nm,4,size(s2,2));
				set(h,'cdata',s2)

				% adjust axis for better viewing
 				axis equal

			end

			% for future reference
			v=axis;

			% for while loop initial status
			ctr=[1e9 1e9];

			% move ball
			if numel(varargin)

				% get handle of ball
				h=varargin{1};

				% create a patch that spans the axis edges (acts as ground)
				% (z) is flat, and depends on the ball's min (z)
				pft=min(min(get(h,'zdata')));
				aaa=[v(1) v(3) pft;
					v(2) v(3) pft;
					v(2) v(4) pft;
					v(1) v(4) pft;
					v(1) v(3) pft;
					v(2) v(3) pft;
					v(2) v(4) pft;
					v(1) v(4) pft];
				b=[1 2 6 5;2 3 7 6;3 4 8 7;4 1 5 8;1 2 3 4;5 6 7 8];

				% plot ground!
				pn=patch('vertices',aaa,'faces',b,'edgecolor',[1 1 1],'facecolor',[.3 .3 .3],...
					'facevertexalphadata',0.5,'facealpha','flat');
			end

			% get min (z) of ball (i.e. get the current floor)
			llpp=get(h,'zdata');
			cf=min(min(llpp));

			% while centre of the ball is at least .1 away from
			% the destination points, keep moving it. notice here that the
			% condition for stopping the ball movement doesn't depend on the
			% next location, but rather on the current location. this leads to
			% the balls to end up in slightly different locations that their
			% destinations, due to the randomization of the ball radius coupled
			% with their relatively high speeds
			while mean(abs(ctr-B))>1e-1

				% get centre of ball
				ctr=[(max(max(get(h,'xdata')))+min(min(get(h,'xdata'))))/2,...
					(max(max(get(h,'ydata')))+min(min(get(h,'ydata'))))/2];

				% vector from current point to destination
				rot=B-ctr;

				% get orthogonal vector (which the ball will spin around)
				yu=rot(2);
				rot(2)=rot(1);
				rot(1)=-yu;

				% rotate around ball about said vector
				% origin is centre of ball
				rotate(h,[rot,0],V,[ctr,r+cf])

				% move ball according to its speed
				% the sign commands are for the direction of movement
				% the sphere moves in X and Y the same distance as the one it has rolled
				% in other words it moves V/360
				% but it doesn't move the entirety of this roll since it is divided between
				% X and Y rolling, we find out how much is shared between X and Y using atan
				set(h,'xdata',get(h,'xdata')...
					-(sign(rot(1))-~sign(rot(1)))...
					*V/360*...
					atan(rot(2)/-rot(1)))
				set(h,'ydata',get(h,'ydata')+(sign(rot(2))+sign(rot(2)))*V/360*atan(rot(1)/-rot(2)))

				% so that axis wouldn't change while animating
				axis(v)

				% update plot
				drawnow()
			end
		end
	end
% EOF
end


text:
function prob_test0744 ( )

%*****************************************************************************80
%
%% TEST0744 tests FRECHET_CDF, FRECHET_CDF_INV and FRECHET_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
  seed = 1213456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0744\n' );
  fprintf ( 1, '  For the Frechet PDF:\n' );
  fprintf ( 1, '  FRECHET_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  FRECHET_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  FRECHET_PDF evaluates the PDF;\n' );

  alpha = 3.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter ALPHA =         %f\n', alpha );

  fprintf ( 1, '\n' );
  fprintf ( 1,'       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = frechet_sample ( alpha, seed );

    pdf = frechet_pdf ( x, alpha );

    cdf = frechet_cdf ( x, alpha );

    x2 = frechet_cdf_inv ( cdf, alpha );

    fprintf ( 1, '  %12f  %12f  %12f  %12f\n',x, pdf, cdf, x2 );

  end

  return
end

author: johannesgerer


text:
function dist = mrFloodFillVolume(calc, sagSize, numSlices, dataRange,startPt)
%
% mrFloodFillVolume
%
%	dist = mrFloodFillVolume(calc, sagSize, numSlices, dataRange,startPt)
%
%	Calculate minimum distances to pts from starting point while remaining
%	in the region specified by the volume calc.  Uses a
%	flood fill algorithm.
%

dist = zeros(1,length(calc));

active = [startPt,1]';
calc(mr3d21d(active(1:3),sagSize,dataRange)) = 0;
thepixel = active';
count = 0;
sum(calc ~= 0)

while ~isempty(active)
	[themin,mindex] = min(active(4,:));
	if (mod(count,100) == 0)
		disp(count)
	end
	count = count+1;
	thepixel = active(:,mindex);
	dist(mr3d21d(thepixel(1:3),sagSize,dataRange)) = thepixel(4);
	active(:,mindex) = [];
	for i = -1:1:1
	   for j = -1:1:1
		for k = -1:1:1
  		    nupixel = [thepixel(1)+i,thepixel(2)+j,thepixel(3)+k];
	  	    if((nupixel(3) >= dataRange(1)) & (nupixel(3) <= dataRange(2)))
			foo = mr3d21d(nupixel,sagSize,dataRange);
			if(calc(foo))
			    if(count == 1)
				foo
			    end
			    calc(foo) = 0;
			    a = i*.9375; b = j*.9375; c = k*.7;  % Convert to mm
			    tmp = thepixel(4)+ sqrt(a*a+b*b+c*c);
			    newpixel = [nupixel,tmp]';
			    active = [active,newpixel];
			end
		    end
		end
	    end
	end
end

author: vistalab


text:
% std_interp() - interpolate, if needed, a list of named data channels
%                for all datasets included in a STUDY. Currently assumes
%                that all channels have uniform locations across datasets.
%
% Usage:      >> [STUDY ALLEEG] = std_interp(STUDY, ALLEEG, chans, method);
%
% Inputs:
%     STUDY    - EEGLAB STUDY structure
%     ALLEEG   - EEGLAB vector containing all EEG dataset structures in the STUDY.
%     chans    - [Cell array] cell array of channel names (labels) to interpolate
%                into the data if they are missing from one of the datasets.
%     method   - [string] griddata() method to use for interpolation.
%                See  >> help eeg_interp() {default:'spherical'}
%
% Important limitation:
%     This function currently presuposes that all datasets have the same channel
%     locations (with various channels from a standard set possibly missing).
%     If this is not the case, the interpolation will not be performed.
%
% Output:
%     STUDY    - study structure.
%     ALLEEG   - updated datasets.
%
% Author: Arnaud Delorme, CERCO, CNRS, August 2006-
%
% See also: eeg_interp()

% Copyright (C) Arnaud Delorme, CERCO, 2006, 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

% Deprecated:
%              - [chanlocs structure] channel location structure containing
%                a full channel structure (missing channels in the current
%                dataset are interpolated).

function [STUDY, ALLEEG] = std_interp(STUDY, ALLEEG, chans, method);

if nargin < 2
    help std_interp;
    return;
end;
if nargin < 3
    chans = [];
end;
if nargin < 4
    method = 'spherical';
end;

% union of all channel structures
% -------------------------------
alllocs = eeg_mergelocs(ALLEEG.chanlocs);

% check electrode names to interpolate
% ------------------------------------
if iscell(chans)
    alllabs = lower({ alllocs.labels });
    for index = 1:length(chans)
        tmpind = strmatch(lower(chans{index}), alllabs, 'exact');
        if isempty(tmpind)
            error( sprintf('Channel named ''%s'' not found in any dataset', chans{index}));
        end;
    end;
end;

% read all STUDY datasets and interpolate electrodes
% ---------------------------------------------------
for index = 1:length(STUDY.datasetinfo)
   tmpind  = STUDY.datasetinfo(index).index;
   tmplocs = ALLEEG(tmpind).chanlocs;

   % build electrode location structure for interpolation
   % ----------------------------------------------------
   [tmp tmp2 id1] = intersect({tmplocs.labels}, {alllocs.labels});
   if isempty(chans)
       interplocs = alllocs;
   elseif iscell(chans)
       [tmp tmp2 id2] = intersect( chans, {alllocs.labels});
       interplocs = alllocs(union(id1, id2));
   else
       interplocs = chans;
   end;

   if length(interplocs) ~= length(tmplocs)

       % search for position of electrode in backup structure
       % ----------------------------------------------
       extrachans = [];
       if isfield(ALLEEG(tmpind).chaninfo, 'nodatchans')
            if isfield(ALLEEG(tmpind).chaninfo.nodatchans, 'labels')
                extrachans = ALLEEG(tmpind).chaninfo.nodatchans;
            end;
       end;
       tmplabels = { tmplocs.labels };
       for i=1:length(interplocs)
           ind = strmatch( interplocs(i).labels, tmplabels, 'exact');
           if ~isempty(ind)
               interplocs(i) = tmplocs(ind); % this is necessary for polhemus
           elseif ~isempty(extrachans)
               ind = strmatch( interplocs(i).labels, { extrachans.labels }, 'exact');
               if ~isempty(ind)
                    fprintf('Found position of %s in chaninfo structure\n', interplocs(i).labels);
                    interplocs(i) = extrachans(ind);
               end;
           end;
       end;

        % perform interpolation
        % ---------------------
        EEG = eeg_retrieve(ALLEEG, index);
        EEG = eeg_checkset(EEG);
        EEG = eeg_interp(EEG, interplocs, method);
        EEG.saved = 'no';
        EEG = pop_saveset(EEG, 'savemode', 'resave');

        % update dataset in EEGLAB
        % ------------------------
        if isstr(ALLEEG(tmpind).data)
            tmpdata = ALLEEG(tmpind).data;
            [ ALLEEG EEG ] = eeg_store(ALLEEG, EEG, tmpind);
            ALLEEG(tmpind).data  = tmpdata;
            ALLEEG(tmpind).saved = 'yes';
            clear EEG;
        else
            [ ALLEEG EEG ] = eeg_store(ALLEEG, EEG, tmpind);
            ALLEEG(tmpind).saved = 'yes';
        end;
    else
        fprintf('No need for interpolation for dataset %d\n', tmpind);
    end;
end;

function checkchans(STUDY, ALLEEG)

    % Check to see if all the channels have the same coordinates
    % (check only the theta field)
    % ----------------------------------------------------------
    for index = 1:length(STUDY.datasetinfo)
       tmpind  = STUDY.datasetinfo(index).index;
       tmplocs = ALLEEG(tmpind).chanlocs;
       [tmp id1 id2] = intersect({tmplocs.labels}, {alllocs.labels});
       for ind = 1:length(id1)
           if tmplocs(id1(ind)).theta ~= alllocs(id2(ind)).theta

               % find datasets with different coordinates
               % ----------------------------------------
               for ind2 = 1:length(STUDY.datasetinfo)
                   tmplocs2 = ALLEEG(ind2).chanlocs;
                   tmpmatch = strmatch(alllocs(id2(ind)).labels, { tmplocs2.labels }, 'exact');
                   if ~isempty(tmpmatch)
                       if alllocs(id2(ind)).theta == tmplocs2(tmpmatch).theta
                           datind = ind2;
                           break;
                       end;
                   end;
               end;

               error(sprintf( [ 'Dataset %d and %d do not have the same channel location\n' ...
                   'for electrode ''%s''' ], datind, tmpind, tmplocs(id1(ind)).labels));
           end;
       end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
function test_ft_spike_isi()

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_spike_isi

clear
nTrials = 100;
data = [];
time       = linspace(0,1,1000);
data.trial(1:nTrials) = {zeros(1,length(time))};
data.time(1:nTrials) = {time};
for iTrial = 1:nTrials
  for iUnit = 1:3
      data.trial{iTrial}(iUnit,:) = double(rand(1,1000)<0.05);
  end
end
data.fsample = 1000;
data.hdr = [];
data.cfg.trl = [];
data.label{1} = 'chan1';
data.label{end+1} = 'chan2';
data.label{end+1} = 'chan3';

% show that the psth works also with the poisson format
spike = ft_checkdata(data,'datatype', 'spike');
for iUnit  = 1:3
  spike.time{iUnit} =   spike.time{iUnit} + 0.001*rand(1,length(spike.time{iUnit}));
end
%%
% now test the isi
cfgIsi = [];
cfgIsi.keeptrials = 'yes';
isih = ft_spike_isi(cfgIsi,spike);
%%
figure
cfgIsi = [];
cfgIsi.spikechannel = 1;
cfgIsi.plotfit = 'no'
cfgIsi.latency = [0 0.2];
h = ft_spike_plot_isi(cfgIsi,isih)
%%
cfgIsi = [];
cfgIsi.keeptrials = 'yes';
isih = ft_spike_isi(cfgIsi,spike);
%%
cfgIsi = [];
cfgIsi.spikechannel = 1;
isihS = ft_spike_isi(cfgIsi,spike);
%%
cfgIsi = [];
cfgIsi.spikechannel = 'all';
isihS = ft_spike_isi(cfgIsi,spike);
%%
cfgRet = [];
cfgRet.window = 'gausswin'
cfgRet.interpolate = 5
cfgRet.scattersize = 1;
cfgRet.density = 'yes'
cfgRet.spikechannel = 1;
figure
H = ft_spike_plot_isireturn(cfgRet,isihS);
%%
figure
H = ft_spike_plot_isireturn([],isihS);
%%

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cout=comp_irdgtiii(cin,a,M)
%COMP_IRDGTIII  Compute inverse real DGT type III.
%
%   This is a computational routine. Do not call it
%   directly.

%   AUTHOR : Peter L. Søndergaard

N=size(cin,1)/M;
W=size(cin,2);
L=N*a;

cin=reshape(cin,M,N,W);

Mhalf=floor(M/2);

cout=zeros(M,N,W,assert_classname(cin));

for m=0:Mhalf-1
  cout(m+1,:,:)=1/sqrt(2)*(cin(2*m+1,:,:)-i*cin(2*m+2,:,:));
  cout(M-m,:,:)=1/sqrt(2)*(cin(2*m+1,:,:)+i*cin(2*m+2,:,:));
end;

if mod(M,2)==1
  cout((M+1)/2,:,:)=cin(M,:,:);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat