text:
function varargout = process_generate_head( varargin )
% PROCESS_GENERATE_HEAD: Generate s head surface from an MRI.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2013

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Generate head surface';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import anatomy'};
    sProcess.Index       = 20;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/LabelFreeSurfer#The_head_surface_looks_bad';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = 'NewSubject';
    % Option: Number of vertices
    sProcess.options.nvertices.Comment = 'Number of vertices [integer]: ';
    sProcess.options.nvertices.Type    = 'value';
    sProcess.options.nvertices.Value   = {10000, '', 0};
    % Option: Erode factor
    sProcess.options.erodefactor.Comment = 'Erode factor [0,1,2,3]: ';
    sProcess.options.erodefactor.Type    = 'value';
    sProcess.options.erodefactor.Value   = {0, '', 0};
    % Option: Fill factor
    sProcess.options.fillfactor.Comment = 'Fill holes factor [0,1,2,3]: ';
    sProcess.options.fillfactor.Type    = 'value';
    sProcess.options.fillfactor.Value   = {2, '', 0};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return;
    end
    % Number of vertices
    nVertices   = sProcess.options.nvertices.Value{1};
    erodeFactor = sProcess.options.erodefactor.Value{1};
    fillFactor  = sProcess.options.fillfactor.Value{1};
    if isempty(nVertices) || isempty(erodeFactor) || isempty(fillFactor)
        bst_report('Error', sProcess, [], 'Invalid values.');
        return;
    end

    % ===== GET SUBJECT =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    if isempty(iSubject)
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" does not exist.']);
        return
    end
    % Check if a MRI is available for the subject
    if isempty(sSubject.Anatomy)
        bst_report('Error', sProcess, [], ['No MRI available for subject "' SubjectName '".']);
        return
    end

    % ===== GENERATE HEAD =====
    % Generate head surface
    HeadFile = tess_isohead(iSubject, nVertices, erodeFactor, fillFactor);
    % Error handling
    if isempty(HeadFile)
        bst_report('Error', sProcess, [], 'An error occurred in tess_isohead function.');
        return;
    end

    OutputFiles = {'import'};
end

author: brainstorm-tools


text:
BASE_DIR = fullfile(mfilename('fullpath'),'../');
addpath(fullfile(BASE_DIR, '../'));
global_variables;

% viewpoint annotation path
%path_ann_view = '../Annotations';
path_ann_view = fullfile(g_pascal3d_root_folder, 'Annotations');

% read ids of validation images
addpath(fullfile(g_pascal3d_root_folder, 'VDPM'));
addpath(fullfile(g_pascal3d_root_folder, 'PASCAL/VOCdevkit/VOCcode'));
VOCinit;
%pascal_init;
ids = textread(sprintf(VOCopts.imgsetpath, 'val'), '%s');
M = numel(ids);

for i = 1:M
    % read ground truth bounding box
    rec = PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    voc12val_records{i} = rec;
end

save(fullfile(BASE_DIR, 'voc12val_records.mat'),'voc12val_records');

author: ShapeNet

repo: RenderForCNN


text:
function m_full = getFullMask(m,frame_d,delay,kernel)
%GETFULLMASK Convert frame rate mask to a sample-by-sample mask
%
%   M_FULL = IOSR.BSS.GETFULLMASK(M,FRAME_D) expands the time-frequency
%   mask M, which has one unit for each frequency channel and frame of
%   length FRAME_D (in samples), to a sample-by-sample mask. The mask M is
%   a time-frequency mask, with one column for each frequency channel, and
%   one row for each time frame. The resulting mask will have dimensions
%   [FRAME_D*size(M,1) size(M,2)].
%
%   M_FULL = IOSR.BSS.GETFULLMASK(M,FRAME_D,DELAY) removes a delay from
%   each frequency channel in the mask. DELAY is a vector, with the same
%   number of elements as M has columns, containing a delay (in samples)
%   that is removed from the corresponding frequency channel. The mask is
%   subsequently zero-padded.
%
%   M_FULL = IOSR.BSS.GETFULLMASK(M,FRAME_D,DELAY,KERNEL) allows smoothing
%   to be applied to the full mask. By default, the full mask contains
%   rectangular transitions at unit boundaries. Specifying KERNEL allows
%   the transitions to be smoother, by convolving the full mask with a
%   two-dimensional kernel (dimensions [frequency time]). The central part
%   of the convolution is returned, so the centre of the KERNEL should be
%   in the centre of the matrix. The kernel is normalised in order to
%   ensure zero gain at DC.
%
%   See also IOSR.BSS.RESYNTHESISE.

%   Copyright 2016 University of Surrey.

    if nargin < 2
        error('iosr:getFullMask:nargin','Not enough input arguments')
    end

    numchans = size(m,2);
    frameCount = size(m,1);

    if nargin < 3
        delay = zeros(1,numchans);
    end
    if nargin < 4
        kernel = 1;
    end

    % Create the sample-by-sample mask
    m_full = zeros(frameCount*frame_d,numchans);
    for i = 1:numchans
        for j = 1:frameCount
            m_full(((j-1)*frame_d+1):((j-1)*frame_d+1)+frame_d-1,i) = m(j,i);
        end
        m_full(:,i) = [m_full(delay(i)+1:end,i); zeros(delay(i),1)];
    end

    % convolve with kernel
    kernel = kernel./sum(abs(kernel(:)));
    m_full = conv2(m_full,kernel,'same');

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function T=findDartboardSubarrays(aperFracs,subarraysPerRing,angShifts,xyPoints)
%%FINDDARTBOARDSUBARRAYS When designing subarrays for a circular aperture,
%          it is noted in [1] that a dartboard pattern tends to work pretty
%          well. It is also suggested that a subarray pattern returned by
%          findBaylissSubarrays be used as a rule of thumb for choosing how
%          many rings and the number of subarrays per ring in the
%          dartboard. When given the positions of elements in a
%
%INPUTS: aperFracs A numRingsX1 or 1XnumRings vector sorted in increasing
%                  order where aperFrac(i) is the fraction of the aperture
%                  size where the ring specified ends. aperFrac(end) must
%                  equal 1 (100% of the aperture size). The aperture width
%                  is taken to be the maximum distance of a point in
%                  xyPoints from the origin.
% subarraysPerRing A numRingsX1 or 1XnumRings vector indicating the number
%                  of subarrays in each ring. These must be positive
%                  integers.
%        angShifts A numRingsX1 or 1XnumRings vector giving angular offsets
%                  in radians by which the first subarray in each ring
%                  begins. This is here to keep the rings from all lining
%                  up, which should improve how the grating lobes of the
%                  subarrays line up. If an empty matrix is passed, then
%                  angular shifts of zero are used.
%         xyPoints A 2XnumPoints set of numPoints points in the aperture
%                  plane corresponding to element positions. The points
%                  should be shifted such that the center of the aperture
%                  is the origin. The width of the aperture is the maximum
%                  distance of any point from the center.
%
%OUTPUTS: T A numSubarraysXnumPoints boolean matrix (a matrix such that
%           T(i,:) is a set of boolean values indicating which elements are
%           in each subarray). The subarrays do not overlap.
%
%EXAMPLE:
%Here, we lay out the elements for a circular array and then we break them
%into a dartboard configuration. We then display the elements with
%different symbols and colors for the subarrays.
% %First, we create a circular array. The element locations are given in
% %terms of the wavelength lambda, so lambda will not appear in the
% %equations for the sum beam.
% xyVals=getShaped2DLattice([49;49],'circular');
%
% aperFracs=[1/5;2/5;3/5;4/5;1];
% subarraysPerRing=[4;4;8;8;8];
% degShifts=[0;2*pi/(2*4);0;2*pi/(2*8);0];
% T=findDartboardSubarrays(aperFracs,subarraysPerRing,degShifts,xyVals);
%
% figure()
% clf
% hold on
% axis square
% numSubarrays=size(T,1);
% dispOpts={'.b','.g','.m','.c','.y','.k','.b','.r',...
%     'xm','xc','xy','xk','xb','xr','xb','xg',...
%     'ob','og','om','oc','oy','ok','ob','or'...
%     '+m','+c','+m','+y','+k','+r','+b','+g'...
%     'sb','sg','sm','sc','sm','sy','sk','sr'};
% for curSubarray=1:numSubarrays
%     points=xyVals(:,T(curSubarray,:)==1);
%     optVal=mod(curSubarray,length(dispOpts))+1;
%     scatter(points(1,:),points(2,:),dispOpts{optVal},'linewidth',2);
% end
% h1=xlabel('x');
% h2=ylabel('y');
% title('Coloring Represents Subarrays')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] U. Nickel, "Subarray configurations for digital beamforming with
%    low sidelobes, adaptive interference suppression and superresolution,"
%    1995, FFM-Report (available from Fraunhofer FKIE, Wachtberg, Germany).
%
%September 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numRings=length(aperFracs);
numSubarrays=sum(subarraysPerRing);
numEls=size(xyPoints,2);

if(aperFracs(numRings)~=1)
   error('The last element in  aperFracs must be one');
end

if(isempty(angShifts))
    angShifts=zeros(numRings,1);
end

%The maximum squared distance from the origin to a point is taken to be the
%radius of the aperture squared.
a=sqrt(max(sum(xyPoints.^2,1)));

temp=cumsum(subarraysPerRing(:));
numSubsInPriorRings=[0;temp(1:(numRings-1))];

T=false(numSubarrays,numEls);
for curEl=1:numEls
    %The fractional distance from the outer edge of the aperture where the
    %element is located.
    fracCur=norm(xyPoints(:,curEl))/a;

    %Determine which ring the element is in.
    [~, ringIdx]=binSearch(aperFracs,fracCur,2);

    %Each ring can have multiple subarrays. The subarrays are uniformly
    %spaced in angle, but the starting angle is not taken to be zero but
    %rather degShifts(ringIdx) radians.
    angle=atan2(xyPoints(2,curEl),xyPoints(1,curEl));

    subArrayWidth=2*pi/subarraysPerRing(ringIdx);

    subArrayInRing=fix(wrapRange(angle+angShifts(ringIdx),0,2*pi)/subArrayWidth)+1;
    subArrayIdx=numSubsInPriorRings(ringIdx)+subArrayInRing;
    T(subArrayIdx,curEl)=1;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ t_stop, y_stop ] = p40_stop ( neqn )

%*****************************************************************************80
%
%% P40_STOP returns the stopping point for problem p40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  e = p40_param ( 'GET', 'E', [] );
  t_stop = 1.0;
  y_stop(1) = - 2.0 * sqrt ( e ) ...
    * exp ( ( 1.0 - t_stop * t_stop ) / ( 2.0 * e ) ) / ...
    ( ...
      2.0 * sqrt ( e ) + ...
      ( ...
        erf ( 1.0 / sqrt ( 2.0 * e ) )  + ...
        erf ( t_stop  / sqrt ( 2.0 + e ) ) ...
      ) ...
      * exp ( 0.5 / e ) * sqrt ( 2.0 * pi ) ...
    );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = isTransform3d(trans, varargin)
%ISTRANSFORM3D Check if input is a affine transformation matrix.
%
%   A = isTransform3d(TRANS) where TRANS should be a transformation matrix.
%   The function accepts transformations given using the following formats:
%   [a b c]   ,   [a b c j] , or  [a b c j]
%   [d e f]       [d e f k]       [d e f k]
%   [g h i]       [g h i l]       [g h i l]
%                                 [0 0 0 1]
%
%   If the transformation matrix should only contain rotation and
%   translation without reflection, scaling, shearing, ... set 'rotation'
%   to true. Default is false.
%
%   Example
%     rot = ...
%         createRotationOx(rand*2*pi)*...
%         createRotationOy(rand*2*pi)*...
%         createRotationOx(rand*2*pi);
%     trans = rot*createTranslation3d(rand(1,3));
%     isTransform3d(trans, 'rot', true)
%
%   See also
%   composeTransforms3d, createBasisTransform3d, recenterTransform3d,
%   transformPoint3d

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2018-07-08
% Copyright 2018-2022

narginchk(1,5)

p = inputParser;
logParValidFunc = @(x) (islogical(x) || isequal(x,1) || isequal(x,0));
addParameter(p,'rotation', 0, logParValidFunc);
valTol = @(x) validateattributes(x,{'numeric'},{'scalar', '>=',eps(class(trans)), '<=',1});
addParameter(p,'tolerance', 1e-8, valTol);
parse(p,varargin{:});
rotation = p.Results.rotation;
tolerance = p.Results.tolerance;

% eventually add null translation
if size(trans, 2) == 3
    trans = [trans zeros(size(trans, 1), 1)];
elseif size(trans, 2) < 3 || size(trans, 2) > 4
    a=false;
    return
end

% eventually add normalization
if size(trans, 1) == 3
    trans = [trans;0 0 0 1];
elseif size(trans, 1) < 3 || size(trans, 1) > 4
    a=false;
    return
end

a=true;

% NaN is invalid
if any(isnan(trans(:)))
    a=false;
    return
end

% Infinity is invalid
if any(isinf(trans(:)))
    a=false;
    return
end

% trans(4,4) has to be a one
if ~isequal(1, trans(4,4))
    a = false;
    return
end

% trans(4,1:3) have to be zeros
if ~isequal(zeros(1,3), trans(4,1:3))
    a = false;
    return
end

if rotation
    % transpose(trans(1:3,1:3)) * trans(1:3,1:3) has to be eye(3)
    if any(abs(eye(3) - (trans(1:3,1:3)'*trans(1:3,1:3))) > tolerance)
        a = false;
        return;
    end

    % determinant of trans(1:3) has to be one
    if abs(1-det(trans)) > tolerance
        a = false;
        return
    end
end

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.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.

function CstActivateLocalWCS(mws, SetNormal, SetOrigin, SetUVector, Activate)

%SetNormal = [0 0 1]
%SetOrigin = [0 0 0]
%SetUVector = [1 0 0]
%Activate = 1 or 0;

if Activate == 1
WCS = invoke(mws,'WCS');
invoke(WCS,'ActivateWCS','local');
invoke(WCS,'SetNormal',num2str(SetNormal(1)),num2str(SetNormal(2)),num2str(SetNormal(3)));
invoke(WCS,'SetOrigin',num2str(SetOrigin(1)),num2str(SetOrigin(2)),num2str(SetOrigin(3)));
invoke(WCS,'SetUVector',num2str(SetUVector(1)),num2str(SetUVector(2)),num2str(SetUVector(3)));
end
if Activate == 0
  WCS = invoke(mws,'WCS');
  invoke(WCS,'ActivateWCS','Global');
end

end

repo: CST-MATLAB-API


text:
function L = laplacian( F )
%LAPLACIAN Vector Laplacian of a CHEBFUN2V.
%   LAPLACIAN(F) returns a CHEBFUN2V representing the vector Laplacian of F.
%
% See also CHEBFUN2V/LAP.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty( F ) )
    L = chebfun2v;
    return
end

% laplacian = f_xx + f_yy
L = diff(F, 2, 1) + diff(F, 2, 2);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%JTRAJ Compute a joint space trajectory between two configurations
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, M) is a joint space trajectory Q (MxN) where the joint
% coordinates vary from Q0 (1xN) to QF (1xN).  A quintic (5th order) polynomial is used
% with default zero boundary conditions for velocity and acceleration.
% Time is assumed to vary from 0 to 1 in M steps.  Joint velocity and
% acceleration can be optionally returned as QD (MxN) and QDD (MxN) respectively.
% The trajectory Q, QD and QDD are MxN matrices, with one row per time step,
% and one column per joint.
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, M, QD0, QDF) as above but also specifies initial
% and final joint velocity for the trajectory.
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, T) as above but the trajectory length is defined
% by the length of the time vector T (Mx1).
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, T, QD0, QDF) as above but specifies initial and
% final joint velocity for the trajectory and a time vector.
%
% See also QPLOT, CTRAJ, SerialLink.jtraj.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function [qt,qdt,qddt] = jtraj(q0, q1, tv, qd0, qd1)
    if length(tv) > 1
        tscal = max(tv);
        t = tv(:)/tscal;
    else
        tscal = 1;
        t = (0:(tv-1))'/(tv-1); % normalized time from 0 -> 1
    end

    q0 = q0(:);
    q1 = q1(:);

    if nargin == 3
        qd0 = zeros(size(q0));
        qd1 = qd0;
    elseif nargin == 5
        qd0 = qd0(:);
        qd1 = qd1(:);
    else
        error('incorrect number of arguments')
    end

    % compute the polynomial coefficients
    A = 6*(q1 - q0) - 3*(qd1+qd0)*tscal;
    B = -15*(q1 - q0) + (8*qd0 + 7*qd1)*tscal;
    C = 10*(q1 - q0) - (6*qd0 + 4*qd1)*tscal;
    E = qd0*tscal; % as the t vector has been normalized
    F = q0;

    tt = [t.^5 t.^4 t.^3 t.^2 t ones(size(t))];
    c = [A B C zeros(size(A)) E F]';

    qt = tt*c;

    % compute optional velocity
    if nargout >= 2
        c = [ zeros(size(A)) 5*A 4*B 3*C  zeros(size(A)) E ]';
        qdt = tt*c/tscal;
    end

    % compute optional acceleration
    if nargout == 3
        c = [ zeros(size(A))  zeros(size(A)) 20*A 12*B 6*C  zeros(size(A))]';
        qddt = tt*c/tscal^2;
    end

author: Allopart

repo: rbpf-gmapping


text:
% Digital Video Stabilization and Rolling Shutter Correction using Gyroscopes
% Copyright (C) 2011 Alexandre Karpenko
%
% 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
% 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/>.

function X = track_video_markers(video_file)

% read corresponding movie
xyloObj = mmreader(['data/' video_file '.mov']);
display(xyloObj);

num_frames = xyloObj.NumberOfFrames;
vid_height = xyloObj.Height;
vid_width = xyloObj.Width;

% initialize markers for first frame
figure(1); imshow(read(xyloObj,1));

x=[];
while isempty(x)
    [x, y] = ginput;
    display([x y]);
end

rect_size = 160;

X = cell(numel(x), 1);

% Read one frame at a time.
for f = 1:num_frames
    frame = read(xyloObj, f);

    for i = 1:numel(x)
        min_x = max(round(x(i) - rect_size/2), 1);
        max_x = min(round(x(i) + rect_size/2), vid_width);

        min_y = max(round(y(i) - rect_size/2), 1);
        max_y = min(round(y(i) + rect_size/2), vid_height);

        patch = double(rgb2gray(frame(min_y:max_y, min_x:max_x, :)));
        min_patch = min(patch(:));
        patch = (patch - min_patch) / (max(patch(:)) - min_patch);

        figure(2); subplot(numel(x), 1, i);
        imshow(patch < 0.2);

        [py, px] = find(patch < 0.2);
        x(i) = min_x + mean(px);
        y(i) = min_y + mean(py);
        X{i} = [X{i}; x(i) y(i)];
    end

    figure(1);
    imshow(frame); hold on;
    colors = {'rx', 'gx', 'bx', 'cx'};
    for i = 1:numel(x)
        plot(x(i), y(i), colors{mod(i-1,4)+1});
    end
    hold off;
end

author: alex-golts

repo: Video-Stabilization


text:
function [event] = read_nexstim_event(filename)

% Use as
%   [event] = read_nexstim_event(filename)

% Written by Vladimir Litvak based on the function nxeGetTriggers
% provided by Nexstim
%
% Copyright (C) 2007, Vladimir Litvak
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_nexstim_event.m 945 2010-04-21 17:41:20Z roboos $

% trigLine - either 1(GATE), 2(TRIG1) or 3(TRIG2)
% trigEdge - either 'rising' or 'falling'

fid=fopen(filename,'r','l');

numChannels = 64;
blockSamples = 14500;
trigThreshold = 2000;
trigChannels = [1 2 3]; % Is fixed for the present Nexstim format.

fseek(fid,0,'eof');
numBytes = ftell(fid);
numSamples = (numBytes/2)/numChannels;
numBlocks = ceil(numSamples/blockSamples);

fseek(fid,0,'bof');
trigPos=[];

event=[];
for i = 1:numBlocks
    blockPos=(i-1)*blockSamples;
    data = fread(fid,[numChannels blockSamples+1],'int16');

    for trigLine=1:length(trigChannels);
        trigPosRising = find(diff(data(trigChannels(trigLine),:))>trigThreshold)+blockPos;
        trigPosFalling = find(diff(data(trigChannels(trigLine),:))<-trigThreshold)+blockPos;

        if ~isempty(trigPosRising)
            for i=1:length(trigPosRising)
                event(end+1).type     = 'rising';
                event(end  ).sample   =  trigPosRising(i);
                event(end  ).value    = trigLine;
                event(end  ).offset   = [];
                event(end  ).duration = [];
            end
        end

        if ~isempty(trigPosFalling)
            for i=1:length(trigPosFalling)
                event(end+1).type     = 'falling';
                event(end  ).sample   =  trigPosFalling(i);
                event(end  ).value    = trigLine;
                event(end  ).offset   = [];
                event(end  ).duration = [];
            end
        end

    end
    fseek(fid,-(2*numChannels),'cof');
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s10.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function s10(x)
% Shekel10 function
function f = s10(x)
a = [4, 1, 8, 6, 3, 2, 5, 8, 6, 7;
     4, 1, 8, 6, 7, 9, 5, 1, 2, 3.6;
     4, 1, 8, 6, 3, 2, 3, 8, 6, 7;
     4, 1, 8, 6, 7, 9, 3, 1, 2, 3.6];
c = [ 0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3, 0.7, 0.5, 0.5];
if size(x,1) == 1
 x = x';
end
for i=1:10
 b = (x - a(:,i)).^2;
 d(i) = sum(b);
end
f = -sum((c+d).^(-1));

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
%% SOBOL_TEST02 tests BIT_HI1_BASE_2.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'SOBOL_TEST02\n' );
fprintf ( 1, '  BIT_HI1_BASE_2 returns the location of the high 1 bit.\n' );

fprintf ( 1, '\n' );
fprintf ( 1, '     I     BIT_HI1_BASE_2\n' );
fprintf ( 1, '\n' );

seed = get_seed ( 0 );

for ( test = 1 : 10 )

  [ i, seed ] = i_random ( 0, 100, seed );

  j = bit_hi1_base_2 ( i );

  fprintf ( 1, '%6d %6d %6d\n', i, j );

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [pos, tri] = mesh_tetrahedron

% MESH_TETRAHEDRON returns the vertices and triangles of a tetrahedron.
%
% Use as
%   [pos, tri] = mesh_tetrahedron;
%
% See also MESH_ICOSAHEDRON, MESH_OCTAHEDRON, MESH_SPHERE

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

v1 = [ 0, 0, 1 ];
v2 = [  sqrt(8/9),          0, -1/3 ];
v3 = [ -sqrt(2/9),  sqrt(2/3), -1/3 ];
v4 = [ -sqrt(2/9), -sqrt(2/3), -1/3 ];

pos = [
  v1
  v2
  v3
  v4
  ];

tri = [
  1 2 3
  1 3 4
  1 4 2
  2 4 3
  ];

  % scale all vertices to the unit sphere
  pos = pos ./ repmat(sqrt(sum(pos.^2,2)), 1,3);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [DAT, StudyTable] = struct2canlab_dataset(data_struct)
% Take an ad hoc data structure data_struct and return a Table format
% object (StudyTable) and a canlab_dataset object DAT.
%
% - field names become variable names
% - Column vectors of numbers or text becomes Subj_Level data
% - Cell vectors containing column vectors of trial-level numeric data become Event_level data
% - Number of elements in vector stored in each field should match num. observations (subjects) for Subj_level array data
% - Trial-level data should contain same-length vectors for all variables for a given subject
%
% See also struct2table and table2canlab_dataset

N = fieldnames(data_struct);

% Empty structures to sort into
CELLS = struct();
COLUMNS = struct();

for i = 1:length(N)

    myvar = data_struct.(N{i});

    % Try to enforce columnar orientation
    if ~iscolumn(myvar), myvar = myvar'; end

    % If still not a column, skip this (matrix?)
    if ~iscolumn(myvar), continue; end

    if iscell(myvar)

        if all(cellfun(@length, myvar) == 1) && ~any(cellfun(@ischar, myvar))
            % all single values, really numbers
            myvar = cat(1, myvar{:});
            COLUMNS.(N{i}) = myvar;
        else
            CELLS.(N{i}) = myvar;
        end

    else
        % Text or event-level data
        COLUMNS.(N{i}) = myvar;

    end

end

%% Get sizes and remove odd-sized

CELLS = remove_odd_sized(CELLS);

COLUMNS = remove_odd_sized(COLUMNS);

% ------------------------------------------------------------------------

%% Make table

StudyTable = struct2table(COLUMNS);

%% Initialize canlab_dataset
% Deal with subject-level data in matrix first

DAT = table2canlab_dataset(StudyTable);

%% Add text data

DAT = add_text_data_from_CELLS(DAT, CELLS);

%% Add event-level data from cells

CELLS = remove_text_data_from_CELLS_struct(DAT, CELLS);

N = fieldnames(CELLS);
nobs = length(CELLS.(N{1}));

% Build array of event-level data
% Deal data from each of nobs into a matrix of vars for each obs
XX = cell(1, nobs);

for i = 1:length(N)

    myvar = CELLS.(N{i});
    XXnames(1, i) = N(i);

    for j = 1:nobs

        datatoadd =  myvar{j};

        if ~iscolumn(datatoadd), datatoadd = datatoadd'; end
        if ~iscolumn(datatoadd), continue; end

        if ~isempty(XX{j}) && size(XX{j}, 1) ~= size(datatoadd, 1)
            continue
        end

        XX{j}(:, i) = datatoadd;

    end

end

DAT.Event_Level.names = XXnames;
DAT.Event_Level.data = XX;
DAT.Event_Level.type = repmat({'numeric'}, 1, length(XXnames));

end % function

function CELLS = remove_odd_sized(CELLS)

N = fieldnames(CELLS);

for i = 1:length(N)

    cellsize(i) = size(CELLS.(N{i}), 1);

end

cellmode = mode(cellsize);
cellmode = cellmode(1);

wh_omit = cellsize ~= cellmode;

for i = 1:length(N)

    if wh_omit(i)

        CELLS = rmfield(CELLS, N{i});

    end

end

end

function DAT = add_text_data_from_CELLS(DAT, CELLS)

nobs = 0;
if ~isempty(DAT.Subj_Level.data), nobs = size(DAT.Subj_Level.data, 1); end

N = fieldnames(CELLS);

for i = 1:length(N)

    cellsize = size(CELLS.(N{i}), 1);

    % Skip if Subj_level nobs is defined and this does not match
    if nobs && cellsize ~= nobs
        continue
    end

    istext = all(cellfun(@ischar, CELLS.(N{i})));

    if istext

        nextvar = max(size(DAT.Subj_Level.data, 2), size(DAT.Subj_Level.textdata, 2));

        DAT.Subj_Level.textdata(:, nextvar) = CELLS.(N{i});
        DAT.Subj_Level.names(nextvar) = N(i);
        DAT.Subj_Level.type(nextvar) = {'text'};

    end

end

end

function CELLS = remove_text_data_from_CELLS_struct(DAT, CELLS)

nobs = 0;
if ~isempty(DAT.Subj_Level.data), nobs = size(DAT.Subj_Level.data, 1); end

N = fieldnames(CELLS);

for i = 1:length(N)

    cellsize = size(CELLS.(N{i}), 1);

    % Remove if Subj_level nobs is defined and this does not match
    if nobs && cellsize ~= nobs

        CELLS = rmfield(CELLS, N{i});

    end

    istext = all(cellfun(@ischar, CELLS.(N{i})));

    if istext

        CELLS = rmfield(CELLS, N{i});

    end

end

end


text:
function make_srrc_lut
    OS_RATE = 8;

    f = firrcos(2*OS_RATE,.25,.25,1,'rolloff','sqrt');
    f = f/sum(abs(f)); % make sure no matter what we don't go beyond 1

    fid = fopen('SRRC.m','w+');
    fprintf(fid,'function y = SRRC\n');
    fprintf(fid,'%%#codegen\n');
    fprintf(fid,'y = [\n');
    fprintf(fid,'%13.12f\n',f);
    fprintf(fid,'];\n');
    fclose(fid);
end

author: Sable

repo: mcbench-benchmarks


text:
function plotPDF(SO3F,h,varargin)
% plot pole density function
%
% Syntax
%   plotPDF(odf,[h1,..,hN])
%   plotPDF(odf,{[h11,h12],h2,hN],'superposition',{[c11,c12],c2,cN})
%   plotPDF(odf,pf.h,'superposition',pf.c)
%
% Input
%  odf - @SO3FunHarmonic
%  h   - @Miller crystallographic directions
%  c   - structure coefficients
%
% Options
%  resolution    - resolution of the plots
%  superposition - plot superposed pole figures
%
% Flags
%  noTitle   - suppress the Miller indices at the top
%  antipodal - include <VectorsAxes.html antipodal symmetry>
%  complete  - plot entire (hemi)--sphere
%
% See also
% S2Grid/plot annotate savefigure Plotting Annotations_demo ColorCoding_demo PlotTypes_demo
% SphericalProjection_demo

if numel(SO3F)>1
  warning(['You try to plot an multivariate function. Plot the desired components ' ...
    'manually. In the following the first component is plotted.'])
end

plotPDF@SO3Fun(SO3F.subSet(1),h,varargin{:});

end

author: mtex-toolbox

repo: mtex

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


text:
function [xi,yi] = snakeinterp(x,y,dmax,dmin)
%----------This function is used for interpolating the snake adaptively

%   dmax: the maximum distance between two snake points
%   dmin: the maximum distance between two snake points
%   d(i,i+1)>dmax, then a new point is added between i and i+1
%   d(i,i+1)<dmin, then either i or i+1 is removed
%

x = x(:); y = y(:);

N = length(x);

d = abs(x([2:N 1])- x(:)) + abs(y([2:N 1])- y(:));

% remove the points which distance to neighbor points is shorter than dmin
IDX = (d<dmin);

idx = find(IDX==0);
x = x(idx);
y = y(idx);

N = length(x);
d = abs(x([2:N 1])- x(:)) + abs(y([2:N 1])- y(:));

IDX = (d>dmax);

z = snakeindex(IDX);

p = 1:N+1;

xi = interp1(p,[x;x(1)],z');
yi = interp1(p,[y;y(1)],z');

N = length(xi);
d = abs(xi([2:N 1])- xi(:)) + abs(yi([2:N 1])- yi(:));

while (max(d)>dmax),

    IDX = (d>dmax);
    z = snakeindex(IDX);

    p = 1:N+1;

    xi = interp1(p,[xi;xi(1)],z');
    yi = interp1(p,[yi;yi(1)],z');

    N = length(xi);
    d = abs(xi([2:N 1])- xi(:)) + abs(yi([2:N 1])- yi(:));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


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

function z = cell2nosyms(x)

  for i=1:length(x)
    if isa(x{i}, 'sym')
      x{i} = double(x{i});
    end
  end

  z = x;

end

repo: octsympy


text:
classdef TestImwrite
    %TestImwrite

    properties (Constant)
        im = cv.imread(fullfile(mexopencv.root(),'test','fruits.jpg'), 'ReduceScale',2);
    end

    methods (Static)
        function test_write
            frmts = TestImwrite.getFormats();
            for i=1:numel(frmts)
                filename = [tempname() frmts(i).ext];
                cObj = onCleanup(@() TestImwrite.deleteFile(filename));
                try
                    if strcmp(frmts(i).ext, '.exr')
                        img = single(TestImwrite.im) / 255;
                    else
                        img = TestImwrite.im;
                    end
                    if any(strcmp(frmts(i).ext, {'.pbm', '.pgm'}))
                        img = cv.cvtColor(img, 'RGB2GRAY');
                    end
                    cv.imwrite(filename, img, frmts(i).opts{:});
                    assert(exist(filename,'file')==2, ...
                        'Failed to write %s', frmts(i).name);
                catch ME
                    %TODO: some codecs are not available on all platforms
                    %error('mexopencv:testskip', 'codecs %s (%s)', ...
                    %    frmts(i).name, frmts(i).ext);
                    continue;
                end
            end
        end

        function test_write_multi
            if true
                fname = which('mri.mat');
                if isempty(fname)
                    error('mexopencv:testskip', 'missing data');
                end
                mri = load(fname);
                imgs = cell(size(mri.D,4), 1);
                for i=1:size(mri.D,4)
                    imgs{i} = ind2gray(mri.D(:,:,1,i), mri.map);
                end
            else
                fnames = dir(fullfile(toolboxdir('images'), 'imdata', 'AT3_*.tif'));
                if isempty(fnames)
                    error('mexopencv:testskip', 'missing data');
                end
                imgs = cell(numel(fnames), 1);
                for i=1:numel(fnames)
                    imgs{i} = imread(fullfile(fnames(i).folder, fnames(i).name));
                end
            end

            filename = [tempname() '.tiff'];
            cObj = onCleanup(@() TestImwrite.deleteFile(filename));

            cv.imwrite(filename, imgs);
            assert(exist(filename,'file')==2, 'Failed to write TIFF');
        end

        function test_verify_lossless_png
            filename = [tempname() '.png'];
            cObj = onCleanup(@() TestImwrite.deleteFile(filename));

            % write image to PNG file (lossless)
            cv.imwrite(filename, TestImwrite.im);
            assert(exist(filename,'file')==2, 'Failed to write PNG');

            % read back image from disk
            img = imread(filename);

            % compare
            assert(isequal(img, TestImwrite.im), 'Images are not equal');
        end

        function test_flip_channels
            filename1 = [tempname() '.png'];
            filename2 = [tempname() '.png'];
            cObj1 = onCleanup(@() TestImwrite.deleteFile(filename1));
            cObj2 = onCleanup(@() TestImwrite.deleteFile(filename2));

            cv.imwrite(filename1, TestImwrite.im, 'FlipChannels',false);
            cv.imwrite(filename2, TestImwrite.im, 'FlipChannels',true);
            im1 = imread(filename1);
            im2 = imread(filename2);
            assert(isequal(im1, flipdim(im2,3)), 'FlipChannels option failed');
        end

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

        function test_error_unrecognized_extension
            %TODO: crashes Octave
            if true
                error('mexopencv:testskip', 'todo');
            end

            try
                cv.imwrite('image.foobar', TestImwrite.im);
                throw('UnitTest:Fail');
            catch e
                %TODO: C++ exception thrown by OpenCV
                %assert(strcmp(e.identifier,'MATLAB:unexpectedCPPexception'));
            end
        end
    end

    %% helper functions
    methods (Static)
        function deleteFile(fname)
            if exist(fname, 'file') == 2
                delete(fname);
            end
        end

        function frmts = getFormats()
            frmts = repmat(struct('name','', 'ext','', 'opts',{{}}), 11, 1);
            frmts(1).name = 'JPEG';
            frmts(1).ext = '.jpg';
            frmts(1).opts = {'JpegQuality',90, 'JpegProgressive',false, ...
                'JpegOptimize',false, 'JpegResetInterval',0, ...
                'JpegLumaQuality',0, 'JpegChromaQuality',0};
            frmts(2).name = 'PNG';
            frmts(2).ext = '.png';
            frmts(2).opts = {'PngCompression',9, 'PngStrategy','RLE'};
            frmts(3).name = 'PPM';
            frmts(3).ext = '.ppm';
            frmts(3).opts = {'PxmBinary',false};
            frmts(4).name = 'TIFF';
            frmts(4).ext = '.tif';
            frmts(5).name = 'BMP';
            frmts(5).ext = '.bmp';
            frmts(6).name = 'WebP';
            frmts(6).ext = '.webp';
            frmts(7).name = 'Sun Raster';
            frmts(7).ext = '.ras';
            frmts(8).name = 'JPEG-2000';
            frmts(8).ext = '.jp2';
            frmts(9).name = 'OpenEXR';
            frmts(9).ext = '.exr';
            frmts(9).opts = {'ExrType','Float'};
            frmts(10).name = 'Radiance HDR';
            frmts(10).ext = '.hdr';
            frmts(11).name = 'PAM';
            frmts(11).ext = '.pam';
            frmts(11).opts = {'PamTupleType','RGB'};
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [ a, ipvt, rcond ] = csico ( a, lda, n )

%*****************************************************************************80
%
%% CSICO factors a complex symmetric matrix.
%
%  Discussion:
%
%    The factorization is done by symmetric pivoting.
%
%    The routine also estimates the condition of the matrix.
%
%    If RCOND is not needed, CSIFA is slightly faster.
%
%    To solve A*X = B, follow CSICO by CSISL.
%
%    To compute inverse(A)*C, follow CSICO by CSISL.
%
%    To compute inverse(A), follow CSICO by CSIDI.
%
%    To compute determinant(A), follow CSICO by CSIDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex A(LDA,N), the symmetric matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N); a block diagonal matrix and the multipliers which
%    were used to obtain it.  The factorization can be written A = U*D*U'
%    where U is a product of permutation and unit upper triangular matrices,
%    U' is the transpose of U, and D is block diagonal with 1 by 1 and
%    2 by 2 blocks.  Only the diagonal and upper triangle are used.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, real RCOND, an estimate of the reciprocal condition of
%    the matrix.  For the system A*X = B, relative perturbations in A and B
%    of size EPSILON may cause relative perturbations in X of size
%    (EPSILON/RCOND).  If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%  Local Parameters:
%
%    Local, complex Z(N), a work vector whose contents are usually
%    unimportant.  If A is close to a singular matrix, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%

%
%  Find norm of A using only upper half.
%
  for j = 1 : n
    z(j) = scasum ( j, a(1:j,j), 1 );
    for i = 1 : j-1
      z(i) = real ( z(i) ) + cabs1 ( a(i,j) );
    end
  end

  anorm = 0.0;
  for j = 1 : n
    anorm = max ( anorm, real ( z(j) ) );
  end
%
%  Factor.
%
  [ a, ipvt, info ] = csifa ( a, lda, n );
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where U*D*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve U*D*W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  k = n;

  while ( 0 < k )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    kp = abs ( ipvt(k) );
    kps = k + 1 - ks;

    if ( kp ~= kps )
      t      = z(kps);
      z(kps) = z(kp);
      z(kp)  = t;
    end

    if ( cabs1 ( z(k) ) ~= 0.0 )
      ek = csign1 ( ek, z(k) );
    end

    z(k) = z(k) + ek;
    z(1:k-ks) = z(1:k-ks) + z(k) * transpose ( a(1:k-ks,k) );

    if ( ks ~= 1 )

      if ( cabs1 ( z(k-1) ) ~= 0.0 )
        ek = csign1 ( ek, z(k-1) );
      end

      z(k-1) = z(k-1) + ek;
      z(1:k-ks) = z(1:k-ks) + z(k-1) * a(1:k-ks,k-1);

    end

    if ( ks ~= 2 )

      if ( cabs1 ( a(k,k) ) < cabs1 ( z(k) ) )
        s = cabs1 ( a(k,k) ) / cabs1 ( z(k) );
        z(1:n) = z(1:n) * s;
        ek = s * ek;
      end

      if ( cabs1 ( a(k,k) ) ~= 0.0 )
        z(k) = z(k) / a(k,k);
      else
        z(k) = 1.0;
      end

    else

      ak = a(k,k) / a(k-1,k);
      akm1 = a(k-1,k-1) / a(k-1,k);
      bk = z(k) / a(k-1,k);
      bkm1 = z(k-1) / a(k-1,k);
      denom = ak * akm1 - 1.0;
      z(k) = ( akm1 * bk - bkm1 ) / denom;
      z(k-1) = ( ak * bkm1 - bk ) / denom;

    end

    k = k - ks;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
%
%  Solve U' * Y = W.
%
  k = 1;

  while ( k <= n )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= 1 )

      z(k) = z(k) + z(1:k-1) * transpose ( a(1:k-1,k) );

      if ( ks == 2 )
        z(k+1) = z(k+1) + z(1:k-1) * transpose ( a(1:k-1,k+1) );
      end

      kp = abs ( ipvt(k) );

      if ( kp ~= k )
        t     = z(k);
        z(k)  = z(kp);
        z(kp) = t;
      end

    end

    k = k + ks;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = 1.0;
%
%  Solve U*D*V = Y.
%
  k = n;

  while ( 0 < k )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= ks )

      kp = abs ( ipvt(k) );
      kps = k + 1 - ks;

      if ( kp ~= kps )
        t      = z(kps);
        z(kps) = z(kp);
        z(kp)  = t;
      end

      z(1:k-ks) = z(1:k-ks) + z(k) * a(1:k-ks,k);

      if ( ks == 2 )
        z(1:k-ks) = z(1:k-ks) + z(k-1) * transpose ( a(1:k-ks,k-1) );
      end

    end

    if ( ks ~= 2 )

      if ( cabs1 ( a(k,k) ) < cabs1 ( z(k) ) )
        s = cabs1 ( a(k,k) ) / cabs1 ( z(k) );
        z(1:n) = z(1:n) * s;
        ynorm = s * ynorm;
      end

      if ( cabs1 ( a(k,k) ) ~= 0.0 )
        z(k) = z(k) / a(k,k);
      else
        z(k) = 1.0;
      end

    else

      ak = a(k,k) / a(k-1,k);
      akm1 = a(k-1,k-1) / a(k-1,k);
      bk = z(k) / a(k-1,k);
      bkm1 = z(k-1) / a(k-1,k);
      denom = ak * akm1 - 1.0;
      z(k) = ( akm1 * bk - bkm1 ) / denom;
      z(k-1) = ( ak * bkm1 - bk ) / denom;

    end

    k = k - ks;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;
%
%  Solve U' * Z = V.
%
  k = 1;

  while ( k <= n )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= 1 )

      z(k) = z(k) + z(1:k-1) * transpose ( a(1:k-1,k) );

      if ( ks == 2 )
        z(k+1) = z(k+1) + z(1:k-1) * transpose ( a(1:k-1,k+1) );
      end

      kp = abs ( ipvt(k) );

      if ( kp ~= k )
        t     = z(k);
        z(k)  = z(kp);
        z(kp) = t;
      end

    end

    k = k + ks;

  end
%
%  Make ZNORM = 1.
%
  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end


text:
function [varargout]=mergeVertices(varargin)

% function [F,V,ind1,ind2]=mergeVertices(F,V,numDigitsMerge)

%% Parse input

switch nargin
    case 2
        F=varargin{1};
        V=varargin{2};
        numDigitsMerge=[];
    case 3
        F=varargin{1};
        V=varargin{2};
        numDigitsMerge=varargin{3};
    otherwise
        error('Wrong number of input arguments');
end

cellMode=isa(F,'cell');

if isempty(numDigitsMerge)
    D=patchEdgeLengths(F,V);
    if cellMode
        mean_D=mean(cellfun(@mean,D)); %Take mean across cell entries
    else
        mean_D=mean(D); %Mean across array
    end
    numDigitsMerge=6-numOrder(mean_D); %base number of digits on mean
end

%% Merge nodes

[~,indKeep,indFix]=unique(pround(V,numDigitsMerge),'rows');
V=V(indKeep,:);

%% Fix indices in face array
if isa(F,'cell')
    for q=1:1:numel(F)
        F{q}=fixFaces(F{q},indFix);
    end
else
    F=fixFaces(F,indFix);
end

%% Collect output

varargout{1}=F;
varargout{2}=V;
varargout{3}=indKeep;
varargout{4}=indFix;

end

%% Fix indices in face array
function F=fixFaces(F,indFix)
    if size(F,1)==1
        F=indFix(F)'; %Fix indices in F
    else
        F=indFix(F); %Fix indices in F
    end
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% ######################################################################################################################################################################################
% We are here talking about spatio-temporal detections, i.e. a set of ground-truth bounding boxes that
%  I will denote by g_t, with t between t_g^b and t_g^e (beginning and end time of the ground-truth)
% versus a detection which is also a set of bounding boxes, denoted by d_t, with t between t_d^e et t_d^e.
%
% a) temporal iou =  T_i / T_u
%  this is the intersection over union between the timing of the the tubes,
% ie mathematically T_i / T_u with
% the intersection T_i = max(0,   max(t_g^b,t_d^b)-min(t_d^e,t_g^e) )
% and the union T_u = min(t_g^b,t_d^b)-max(t_d^e,t_g^e)
%
% b) for each t between max(tgb,tdb)-min(tde,tge), we compute the IoU between g_t and d_t, and average them
%
% Multiplying (a) and (b) is the same as computed the average of the spatial iou over all frames in T_u of the two tubes, with a spatial iou of 0 for frames where only one box exists.
% c) as this is standard in detection problem, if there are multiple detections for the same groundtruth detection, the first one is counted as positive and the other ones as negatives
% ######################################################################################################################################################################################
%{
gt_fnr = 1xn doube
gt_bb = nx4 doubld - [x y w h]
dt_fnr = 1xm double
dt_bb = mx4 double - [x y w h]
%}
% -------------------------------------------------------------------------
function st_iou = compute_spatio_temporal_iou(gt_fnr, gt_bb, dt_fnr, dt_bb)
% -------------------------------------------------------------------------

% time gt begin

tgb = gt_fnr(1);
% time gt end
tge = gt_fnr(end);
%time dt begin
tdb = dt_fnr(1);
tde = dt_fnr(end);
% temporal intersection
T_i = double(max(0, min(tge,tde)-max(tgb,tdb)));

if T_i>0
    T_i = T_i +1;
    % temporal union
    T_u = double(max(tge,tde) - min(tgb,tdb)+1);
    %temporal IoU
    T_iou = T_i/T_u;
    % intersect frame numbers
    int_fnr = max(tgb,tdb):min(tge,tde);

    % find the ind of the intersected frames in the detected frames
    [~,int_find_dt] = ismember(int_fnr, dt_fnr);
    [~,int_find_gt] = ismember(int_fnr, gt_fnr);

    assert(length(int_find_dt)==length(int_find_gt));

    iou = zeros(length(int_find_dt),1);
    for i=1:length(int_find_dt)
        if int_find_gt(i)<1
%             fprintf('error ')
            pf = pf;
        else
            pf = i;
        end

        gt_bound = gt_bb(int_find_gt(pf),:);
        dt_bound = dt_bb(int_find_dt(pf),:)+1;

        % gt_bound = [gt_bound(:,1:2) gt_bound(:,3:4)-gt_bound(:,1:2)];
        % dt_bound = [dt_bound(:,1:2) dt_bound(:,3:4)-dt_bound(:,1:2)];
        iou(i) = inters_union(double(gt_bound),double(dt_bound));
    end
    % finalspatio-temporal IoU threshold
    st_iou = T_iou*mean(iou);
else
    st_iou =0;
end
% % iou_thresh = 0.2,...,0.6 % 'Learing to track paper' takes 0.2 for UCF101 and 0.5 for JHMDB
% if delta >= iou_thresh
%     % consider this tube as valid detection
% end

end

% -------------------------------------------------------------------------
function iou = inters_union(bounds1,bounds2)
% -------------------------------------------------------------------------

inters = rectint(bounds1,bounds2);
ar1 = bounds1(:,3).*bounds1(:,4);
ar2 = bounds2(:,3).*bounds2(:,4);
union = bsxfun(@plus,ar1,ar2')-inters;

iou = inters./(union+eps);

end

author: gurkirt


text:
function [C,R] = incircle(x,y)
% incircle: compute the maximal in-circle of the polygonal convex hull of a set of points in the plane
%
% [C,R] = incircle(x,y)
%    Called with a pair of vectors of the same length,
%    incircle computes the convex hull in 2-dimensions,
%    then finds the maximal radius in-circle that lies
%    entirely within the polygon of the convex hull.
%
% [C,R] = incircle(xy)
%    Called with a single nx2 array, incircle treats
%    this as a list of points in 2 dimensions. Each row
%    of xy is treated as a single point. The convex hull
%    is computed and the maximal in-circle then computed
%    for the convex hull.
%
% [C,R] = incircle(xy,edges)
%    Called with two arguments, the first of which is an
%    nx2 array of data points, and the second a list of
%    edges of the convex hull as would be generated by
%    convhull (or convhulln), the maximal in-circle is
%    generated for the hull as provided. No test is made
%    to assure the hull is truly convex.
%
% All input forms return a row vector of length 2 (C)
% which defines the center of the in-circle, and the
% radius (R) of the in-circle.
%
% Example:
%
% [C,R] = incircle(rand(1000,1),rand(1000,1))
%
% C =
%      0.50191      0.49999
%
% R =
%      0.49788
%
% Example:
%
% x = randn(100,1) + 1;
% y = randn(100,1) - 2.5;
% H = convhull(x,y);
% [C,R] = incircle([x,y],H)
%
% C =
%      1.0212      -2.3458
%
% R =
%      2.1329
%
% t = linspace(0,2*pi,200);
% xc = R*cos(t) + C(1);
% yc = R*sin(t) + C(2);
% plot(x(H),y(H),'-r',x,y,'ko',xc,yc,'-b')
% axis equal
% xlabel X
% ylabel Y
% title 'Convex hull, with polygon in-circle'
% grid on
%
%
% See also: minboundcicle, convhull
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 6/14/09

if (nargin < 1) || (nargin > 2)
  error('INCIRCLE:improperarguments', ...
    'incircle requires exactly 1 or 2 arguments')
elseif (nargin == 2) && isvector(x) && isvector(y) && (numel(x) == numel(y))
  % a pair of vectors
  xy = [x(:),y(:)];
  % compute the hull. I prefer convhulln.
  edges = convhulln(xy);
elseif (nargin == 1) && (size(x,2) == 2)
  % a single list of points as rows of x
  xy = x;
  edges = convhulln(xy);
elseif (nargin == 2) && (size(x,2) == 2) && (size(y,2) == 2)
  % y must be a list of edges from convhulln
  xy = x;
  edges = y;
elseif (nargin == 2) && (size(x,2) == 2) && isvector(y) && (y(1) == y(end))
  % y must be a list of edges from convhull
  xy = x;
  I = y(:);
  % in case the first point was not wrapped in the polygon
  if I(1) ~= I(end)
    I = [I;I(1)];
  end
  edges = [I(1:(end-1)),I(2:end)];
else
  % none of the forms I allow seem to fit.
  % give up and throw an error.
  error('INCIRCLE:invaliddata', ...
    'x and y do not seem to fit into any of the allowed forms for incircle input')
end
ne = size(edges,1);

% the end points of each edge are...
A = xy(edges(:,1),:);
B = xy(edges(:,2),:);

% the normal vector to each edge
N = (B - A)*[0 1;-1 0];

% normalize to unit length
L = sqrt(sum(N.^2,2));
N = N./[L,L];

% a central point inside the hull itself
C0 = mean(A,1);

% ensure the normals point inwards. While
% I can use the signed area for a hull as
% generated by convhull (suggestion by Bruno)
% this test is also efficient, and it deals
% with the case of a hull provided from some
% unknown and less trusted source.
k = sum(N.*bsxfun(@minus,C0,A),2) < 0;
N(k,:) = -N(k,:);

% formulate the linear programming problem.
% given a point C inside the hull, the distance
% from C to the edge that contains A(i,:) is
% given by (dot(N(i,:),C - A(i,:))). If this
% number is negative for any edge of the hull,
% then C is outside of the hull.
%
% Now think of it as a set of slack variables,
% one for each edge of the hull. Given the
% unknown point C that defines the center of
% the in-circle,
%
%  dot(N(i,:),C-A(i,:)) - S(i) == 0
%
% Thus the vector S is of length ne, where ne is
% the number of edges in the convex hull. Every
% element of S must be positive.
%
% Create one more slack variable, a scalar R.
% R is the minimum value that S attains for a
% given point C.
%
%   R >= 0
%   S(i) >= R
%
% The objective for our linear programming problem
% is simple. It is just -R. When we find the
% solution to the LP problem that satisfies the
% equality constraints above between C and S, the
% bound constraint on T, and the inequality
% constraints between S and R, the solution yields
% the maximal radius in-circle that fits inside
% the convex hull polygon.
%
% The unknowns for the LP are a list of ne+3
% variables. The first two variables represent
% the center coordinates of the circle. X(3) = R
% is the circle radius. The remainder of the
% variables are the slack variabls S(i).

% equality constraints defined by the dot products
% with the normal vectors.
Aeq = [N,zeros(ne,1),-eye(ne)];
beq = sum(N.*A,2);

% lower bounds only for the slack variables
LB = [-inf; -inf; 0; zeros(ne,1)];
% there are no upper bounds
UB = [];

% inequalities defined by the slack variable
% constraints
A = [zeros(ne,2),ones(ne,1),-eye(ne)];
b = zeros(ne,1);

% the objective is just -R
f = zeros(ne+3,1);
f(3) = -1;

% just turn off the output message
options = optimset('linprog');
options.Display = 'off';

% linprog does all of the hard work.
result = linprog(f,A,b,Aeq,beq,LB,UB,[],options);

% unpack the circle parameters
C = result(1:2)';
R = result(3);


text:
function [net, info] = cnn_mnist(varargin)
% CNN_MNIST  Demonstrated MatConNet on MNIST

run(fullfile(fileparts(mfilename('fullpath')),...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

opts.expDir = fullfile('data','mnist-baseline') ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.dataDir = 'data/mnist' ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.useBatchNorm = false ;
opts.networkType = 'simplenn' ;
opts.train = struct() ;
[opts, varargin] = vl_argparse(opts, varargin) ;

% --------------------------------------------------------------------
%                                                         Prepare data
% --------------------------------------------------------------------

net = cnn_mnist_init('useBatchNorm', opts.useBatchNorm, ...
                     'networkType', opts.networkType) ;

if exist(opts.imdbPath, 'file')
  imdb = load(opts.imdbPath) ;
else
  imdb = getMnistImdb(opts) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

net.meta.classes.name = arrayfun(@(x)sprintf('%d',x),1:10,'UniformOutput',false) ;

% --------------------------------------------------------------------
%                                                                Train
% --------------------------------------------------------------------

switch opts.networkType
  case 'simplenn', trainfn = @cnn_train ;
  case 'dagnn', trainfn = @cnn_train_dag ;
end

[net, info] = trainfn(net, imdb, getBatch(opts), ...
  'expDir', opts.expDir, ...
  net.meta.trainOpts, ...
  opts.train, ...
  'val', find(imdb.images.set == 3)) ;

% --------------------------------------------------------------------
function fn = getBatch(opts)
% --------------------------------------------------------------------
switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(x,y) ;
  case 'dagnn'
    bopts = struct('numGpus', numel(opts.train.gpus)) ;
    fn = @(x,y) getDagNNBatch(bopts,x,y) ;
end

% --------------------------------------------------------------------
function [images, labels] = getSimpleNNBatch(imdb, batch)
% --------------------------------------------------------------------
images = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;

% --------------------------------------------------------------------
function inputs = getDagNNBatch(opts, imdb, batch)
% --------------------------------------------------------------------
images = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;
if opts.numGpus > 0
  images = gpuArray(images) ;
end
inputs = {'input', images, 'label', labels} ;

% --------------------------------------------------------------------
function imdb = getMnistImdb(opts)
% --------------------------------------------------------------------
% Preapre the imdb structure, returns image data with mean image subtracted
files = {'train-images-idx3-ubyte', ...
         'train-labels-idx1-ubyte', ...
         't10k-images-idx3-ubyte', ...
         't10k-labels-idx1-ubyte'} ;

if ~exist(opts.dataDir, 'dir')
  mkdir(opts.dataDir) ;
end

for i=1:4
  if ~exist(fullfile(opts.dataDir, files{i}), 'file')
    url = sprintf('http://yann.lecun.com/exdb/mnist/%s.gz',files{i}) ;
    fprintf('downloading %s\n', url) ;
    gunzip(url, opts.dataDir) ;
  end
end

f=fopen(fullfile(opts.dataDir, 'train-images-idx3-ubyte'),'r') ;
x1=fread(f,inf,'uint8');
fclose(f) ;
x1=permute(reshape(x1(17:end),28,28,60e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 't10k-images-idx3-ubyte'),'r') ;
x2=fread(f,inf,'uint8');
fclose(f) ;
x2=permute(reshape(x2(17:end),28,28,10e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 'train-labels-idx1-ubyte'),'r') ;
y1=fread(f,inf,'uint8');
fclose(f) ;
y1=double(y1(9:end)')+1 ;

f=fopen(fullfile(opts.dataDir, 't10k-labels-idx1-ubyte'),'r') ;
y2=fread(f,inf,'uint8');
fclose(f) ;
y2=double(y2(9:end)')+1 ;

set = [ones(1,numel(y1)) 3*ones(1,numel(y2))];
data = single(reshape(cat(3, x1, x2),28,28,1,[]));
dataMean = mean(data(:,:,:,set == 1), 4);
data = bsxfun(@minus, data, dataMean) ;

imdb.images.data = data ;
imdb.images.data_mean = dataMean;
imdb.images.labels = cat(2, y1, y2) ;
imdb.images.set = set ;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = arrayfun(@(x)sprintf('%d',x),0:9,'uniformoutput',false) ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function L = ldlrowmodify(L, c2, k)
%LDLROWMODIFY  Function to modify the sparse cholesky factorization
%              L*D*L' = C, when a row and column k of C have changed
%
%  Description
%    L = rowmodify(L, c2, k)
%    Function to modify the Cholesky factorization L*D*L' = C
%    (which is stored in L), when a row and column k of C have
%    changed from c to c2. The change in C is assumed to be such
%    that the sparsity structure of C is remained the same. c is
%    the old kth column of C and c2 is the new column.
%
%    See GPEP_E for usage
%
%  Reference
%    See Davis and Hager 2005 (Row Modification of a sparse
%    Cholesky factorization) section 4 for details of the
%    algorithm.

%  Copyright (c) 2009-2010 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% First evaluate the l_21 vector. That is the kth row vector of L

% Evaluate the old c. NOTE! The C-version (mex-file) of this algorithm
% does not need the old c. However, this m-version is an old
% implementation
    [L2, D] = ldlsplit(L);
    c = L2(k,:)*D(k)*L2(k,:)';

    n = size(L,1);
    if k>1
        deltac12 = c2(1:k-1)-c(1:k-1);
        [L2, D] = ldlsplit(L);

        %deltal12 = D(1:k-1,1:k-1) \ (L2(1:k-1,1:k-1)\deltac12);

        x = deltac12;
        [I,J,ld] = find(L);
        a2 = cumsum(histc(J,1:n));
        Lp=zeros(n,1);
        Lp(1) = 1; Lp(2:end) = a2(1:end-1) + 1; Lp(end+1) = a2(end);
        Li = I;
        Lx = ld;

        for j=1:k-1
            x(j) = x(j);
            for p = Lp(j)+1:Lp(j+1)-1
                if Li(p) >= k
                    break
                end
                x(Li(p)) = x(Li(p))-Lx(p).*x(j);
            end
            x(j) = x(j)./Lx(Lp(j));
        end

        deltal12 = x;

        l12 = L(k,1:k-1)';
        l12b = l12 + deltal12;
        L(k,1:k-1) = l12b;

        d = D(k,k);
        db = d + c2(k)-c(k) - sum( deltal12.*(l12 + l12b).*diag(D(1:k-1,1:k-1)));
        L(k,k) = db;
    else
        d = L(k,k);
        db = d + c2(k)-c(k);
        L(k,k) = db;
    end

    if k < n
        deltac32 = c2(k+1:n)-c(k+1:n);
        l32 = L(k+1:n,k);
        if k>1
            l32b = (deltac32 + l32.*d - L(k+1:n,1:k-1)*D(1:k-1,1:k-1)*deltal12)./db;
        else
            l32b = (deltac32 + l32.*d )./db;
        end

        w = L(:,k).*sqrt(d);
        L(k+1:n,k) = l32b;
        w2 =L(:,k).*sqrt(db);

        % Then update the L_33 part. That is L(k+1:n, k+1:n)
        alpha = 1;
        alpha2 = 1;
        for i = k+1:n    % Li(Lp(k)+1:Lp(k+1)-1)

            if w(i) ~= 0
                beta = alpha + w(i).^2 / L(i,i);
                gamma = w(i) / (beta*L(i,i));
                L(i,i) = (beta/alpha)*L(i,i);
                alpha = beta;

                beta2 = alpha2 - w2(i).^2 /L(i,i);
                gamma2 = w2(i) / (beta2*L(i,i));
                L(i,i) = (beta2/alpha2)*L(i,i);
                alpha2 = beta2;

                if (i == n) return, end

                w(i+1:n) = w(i+1:n) - w(i).*L(i+1:n, i);
                L(i+1:n,i) = L(i+1:n,i) + gamma.*w(i+1:n);

                w2(i+1:n) = w2(i+1:n) - w2(i).*L(i+1:n, i);
                L(i+1:n,i) = L(i+1:n,i) - gamma2.*w2(i+1:n);
            end
        end
    end
end


text:
function spm_eeg_inv_vbecd_disp(action,varargin)
% Display the dipoles as obtained from VB-ECD
%
% FORMAT spm_eeg_inv_vbecd_disp('Init',D)
% Display the latest VB-ECD solution saved in the .inv{} field of the
% data structure D.
%
% FORMAT spm_eeg_inv_vbecd_disp('Init',D, ind)
% Display the ind^th .inv{} cell element, if it is actually a VB-ECD
% solution.
%__________________________________________________________________________
% Copyright (C) 2008-2015 Wellcome Trust Centre for Neuroimaging

% Christophe Phillips
% $Id: spm_eeg_inv_vbecd_disp.m 6501 2015-07-17 14:32:09Z spm $

% Note:
% unfortunately I cannot see how to ensure that when zooming in the image
% the dipole location stays in place...

global st

Fig     = spm_figure('GetWin','Graphics');
colors  = {'y','b','g','r','c','m'};              % 6 possible colors
marker  = {'o','x','+','*','s','d','v','p','h'};  % 9 possible markers
Ncolors = length(colors);
Nmarker = length(marker);

if nargin == 0, action = 'Init'; end

switch lower(action)

%==========================================================================
case 'init'
%==========================================================================
% FORMAT spm_eeg_inv_vbecd_disp('init',D,ind)
% Initialise the variables with GUI
%--------------------------------------------------------------------------

if nargin<2
    D = spm_eeg_load;
else
    D = varargin{1};
end
if nargin<3
    % find the latest inverse produced with vbecd
    Ninv = length(D.inv);
    lind = [];
    for ii=1:Ninv
        if isfield(D.inv{ii},'method') && strcmp(D.inv{ii}.method,'vbecd')
            lind = [lind ii];
        end
    end
    ind = max(lind);
    if ~ind,
        spm('alert*','No VB-ECD solution found with this data file!',...
                'VB-ECD display')
        return
    end
else
    ind = varargin{3};
end

% Stash dipole(s) information in sdip structure
sdip = D.inv{ind}.inverse;

% if the exit flag is not in the structure, assume everything went ok.
if ~isfield(sdip,'exitflag')
    sdip.exitflag = ones(1,sdip.n_seeds);
end

%try
%    Pimg = spm_vol(D.inv{ind}.mesh.sMRI);
%catch
    Pimg = spm_vol(fullfile(spm('dir'), 'canonical', 'single_subj_T1.nii'));
%end

spm_orthviews('Reset');
spm_orthviews('Image', Pimg, [0.0 0.45 1 0.55]);
spm_orthviews('MaxBB');
spm_orthviews('AddContext')
st.callback = 'spm_image(''shopos'');';
% remove clicking in image
for ii=1:3,
    set(st.vols{1}.ax{ii}.ax,'ButtonDownFcn',';');
end

WS = spm('WinScale');

% Build GUI
%==========================================================================
% Location:
%--------------------------------------------------------------------------
uicontrol(Fig,'Style','Frame','Position',[60 25 200 325].*WS, ...
    'DeleteFcn','spm_image(''reset'');');
uicontrol(Fig,'Style','Frame','Position',[70 250 180 90].*WS);
uicontrol(Fig,'Style','Text', 'Position',[75 320 170 016].*WS, ...
    'String','Current Position');
uicontrol(Fig,'Style','Text', 'Position',[75 295 35 020].*WS,'String','mm:');
uicontrol(Fig,'Style','Text', 'Position',[75 275 35 020].*WS,'String','vx:');
uicontrol(Fig,'Style','Text', 'Position',[75 255 75 020].*WS,'String','Img Intens.:');

st.mp = uicontrol(Fig,'Style','Text', 'Position',[110 295 135 020].*WS,'String','');
st.vp = uicontrol(Fig,'Style','Text', 'Position',[110 275 135 020].*WS,'String','');
st.in = uicontrol(Fig,'Style','Text', 'Position',[150 255  85 020].*WS,'String','');

c = 'if get(gco,''Value'')==1, spm_orthviews(''Xhairs'',''off''), else, spm_orthviews(''Xhairs'',''on''); end;';
uicontrol(Fig,'Style','togglebutton','Position',[95 220 125 20].*WS,...
    'String','Hide Crosshairs','Callback',c,'ToolTipString','show/hide crosshairs');

% Dipoles/seeds selection:
%--------------------------------------------------------------------------
uicontrol(Fig,'Style','Frame','Position',[300 25 180 325].*WS);
sdip.hdl.hcl = uicontrol(Fig,'Style','pushbutton','Position',[310 320 100 20].*WS, ...
        'String','Clear all','CallBack','spm_eeg_inv_vbecd_disp(''ClearAll'')');

sdip.hdl.hseed=zeros(sdip.n_seeds,1);
for ii=1:sdip.n_seeds
    if sdip.exitflag(ii)==1
        sdip.hdl.hseed(ii) = uicontrol(Fig,'Style','togglebutton','String',num2str(ii),...
            'Position',[310+rem(ii-1,8)*20 295-fix((ii-1)/8)*20 20 20].*WS,...
            'CallBack','spm_eeg_inv_vbecd_disp(''ChgSeed'')');
    else
        sdip.hdl.hseed(ii) = uicontrol(Fig,'Style','Text','String',num2str(ii), ...
            'Position',[310+rem(ii-1,8)*20 293-fix((ii-1)/8)*20 20 20].*WS) ;
    end
end

uicontrol(Fig,'Style','text','String','Select dipole # :', ...
        'Position',[310 255-fix((sdip.n_seeds-1)/8)*20 110 20].*WS);

txt_box = cell(sdip.n_dip,1);
for ii=1:sdip.n_dip, txt_box{ii} = num2str(ii); end
txt_box{sdip.n_dip+1} = 'all';
sdip.hdl.hdip = uicontrol(Fig,'Style','popup','String',txt_box, ...
        'Position',[420 258-fix((sdip.n_seeds-1)/8)*20 40 20].*WS, ...
        'Callback','spm_eeg_inv_vbecd_disp(''ChgDip'')');

% Dipoles orientation and strength:
%--------------------------------------------------------------------------
uicontrol(Fig,'Style','Frame','Position',[70 120 180 90].*WS);
uicontrol(Fig,'Style','Text', 'Position',[75 190 170 016].*WS, ...
                'String','Dipole orientation & strength');
uicontrol(Fig,'Style','Text', 'Position',[75 165 65 020].*WS, ...
                'String','x-y-z or.:');
uicontrol(Fig,'Style','Text', 'Position',[75 145 75 020].*WS, ...
                'String','theta-phi or.:');
uicontrol(Fig,'Style','Text', 'Position',[75 125 75 020].*WS, ...
                'String','Dip. intens.:');

sdip.hdl.hor1 = uicontrol(Fig,'Style','Text', 'Position', ...
                            [140 165  105 020].*WS,'String','a');
sdip.hdl.hor2 = uicontrol(Fig,'Style','Text', 'Position', ...
                            [150 145  85 020].*WS,'String','b');
sdip.hdl.int  = uicontrol(Fig,'Style','Text', 'Position', ...
                            [150 125  85 020].*WS,'String','c');

st.vols{1}.sdip = sdip;

% First plot = all the seeds that converged !
l_conv = find(sdip.exitflag==1);
if isempty(l_conv)
    error('No seed converged towards a stable solution, nothing to be displayed !')
else
    spm_eeg_inv_vbecd_disp('DrawDip',l_conv,1)
    set(sdip.hdl.hseed(l_conv),'Value',1); % toggle all buttons
end

%==========================================================================
case 'drawdip'
%==========================================================================
% FORMAT spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip,sdip)
% e.g. spm_eeg_inv_vbecd_disp('DrawDip',1,1,sdip)
% e.g. spm_eeg_inv_vbecd_disp('DrawDip',[1:5],1,sdip)
%--------------------------------------------------------------------------

if nargin < 2
    i_seed = 1;
else
    i_seed = varargin{1};
end
if nargin<3
    i_dip = 1;
else
    i_dip = varargin{2};
end

if nargin<4
    if isfield(st.vols{1},'sdip')
        sdip = st.vols{1}.sdip;
    else
        error('I can''t find sdip structure');
    end
else
    sdip = varargin{3};
    st.vols{1}.sdip = sdip;
end

if any(i_seed>sdip.n_seeds) || i_dip>(sdip.n_dip+1)
    error('Wrong i_seed or i_dip index in spm_eeg_inv_vbecd_disp');
end

% Note if i_dip==(sdip.n_dip+1) all dipoles are displayed simultaneously,
% The 3D cut will then be at the mean location of all sources !!!
if i_dip == (sdip.n_dip+1)
    i_dip = 1:sdip.n_dip;
end

% if seed indexes passed is wrong (no convergence) remove the wrong ones
i_seed(sdip.exitflag(i_seed)~=1) = [];
if isempty(i_seed)
    error('You passed the wrong seed indexes...')
end
if size(i_seed,2)==1, i_seed=i_seed'; end

% Display business
%--------------------------------------------------------------------------
loc_mm = sdip.mniloc{i_seed(1)}(:,i_dip);
if length(i_seed)>1
%     unit = ones(1,sdip.n_dip);
    for ii = i_seed(2:end)
        loc_mm = loc_mm + sdip.mniloc{ii}(:,i_dip);
    end
    loc_mm = loc_mm/length(i_seed);
end
if length(i_dip)>1
    loc_mm = mean(loc_mm,2);
end

% Place the underlying image at right cuts

% GRB FIX- removing the line below July 2013
%spm_orthviews('Reposition',loc_mm);

if length(i_dip)>1
    tabl_seed_dip = [kron(ones(length(i_dip),1),i_seed') ...
                        kron(i_dip',ones(length(i_seed),1))];
else
    tabl_seed_dip = [i_seed' ones(length(i_seed),1)*i_dip];
end

% Scaling, according to all dipoles in the selected seed sets.
% The time displayed is the one corresponding to the maximum EEG power !
Mn_j = -1;
l3 = -2:0;
for ii = 1:length(i_seed)
    for jj = 1:sdip.n_dip
        Mn_j = max([Mn_j sqrt(sum(sdip.jmni{ii}(jj*3+l3,sdip.Mtb).^2))]);
    end
end
st.vols{1}.sdip.tabl_seed_dip = tabl_seed_dip;

% Display all dipoles, the 1st one + the ones close enough.
% Run through the 6 colors and 9 markers to differentiate the dipoles.
% NOTA: 2 dipoles coming from the same set will have same colour/marker
ind = 1 ;
dip_h = zeros(9,size(tabl_seed_dip,1),1);
    % each dipole displayed has 9 handles:
    %   3 per view (2*3): for the line, for the circle & for the error
js_m = zeros(3,1);

% Deal with case of multiple i_seed and i_dip displayed.
% make sure dipole from same i_seed have same colour but different marker.

pi_dip = find(diff(tabl_seed_dip(:,2)));
if isempty(pi_dip)
    % i.e. only one dip displayed per seed, use old fashion
    for ii=1:size(tabl_seed_dip,1)
        if ii>1
            if tabl_seed_dip(ii,1)~=tabl_seed_dip(ii-1,1)
                ind = ind+1;
            end
        end
        ic = mod(ind-1,Ncolors)+1;
        im = fix(ind/Ncolors)+1;

        loc_pl = sdip.mniloc{tabl_seed_dip(ii,1)}(:,tabl_seed_dip(ii,2));
        js = sdip.jmni{tabl_seed_dip(ii,1)}(tabl_seed_dip(ii,2)*3+l3,sdip.Mtb);

        vloc = sdip.cov_loc{tabl_seed_dip(ii,1)}(tabl_seed_dip(ii,2)*3+l3,tabl_seed_dip(ii,2)*3+l3);
        dip_h(:,ii) = add1dip(loc_pl,js/Mn_j*20,vloc, ...
                            marker{im},colors{ic},st.vols{1}.ax,Fig,st.bb);
        js_m = js_m+js;
    end
else
    for ii=1:pi_dip(1)
        if ii>1
            if tabl_seed_dip(ii,1)~=tabl_seed_dip(ii-1,1)
                ind = ind+1;
            end
        end
        ic = mod(ind-1,Ncolors)+1;
        for jj=1:sdip.n_dip
            im = mod(jj-1,Nmarker)+1;

            loc_pl = sdip.mniloc{tabl_seed_dip(ii,1)}(:,jj);
            js = sdip.jmni{tabl_seed_dip(ii,1)}(jj*3+l3,sdip.Mtb);
            vloc = sdip.cov_loc{tabl_seed_dip(ii,1)}(jj*3+l3,jj*3+l3);
            js_m = js_m+js;
            dip_h(:,ii+(jj-1)*pi_dip(1)) = ...
                add1dip(loc_pl,js/Mn_j*20,vloc, ...
                        marker{im},colors{ic},st.vols{1}.ax,Fig,st.bb);
        end
    end
end
st.vols{1}.sdip.ax = dip_h;

% Display dipoles orientation and strength
js_m = js_m/size(tabl_seed_dip,1);
[th,phi,Ijs_m] = cart2sph(js_m(1),js_m(2),js_m(3));
Njs_m = round(js_m'/Ijs_m*100)/100;
Angle = round([th phi]*1800/pi)/10;

set(sdip.hdl.hor1,'String',[num2str(Njs_m(1)),' ',num2str(Njs_m(2)), ...
        ' ',num2str(Njs_m(3))]);
set(sdip.hdl.hor2,'String',[num2str(Angle(1)),'  ',num2str(Angle(2)),' ']);
set(sdip.hdl.int,'String',Ijs_m);

% Change the colour of toggle button of dipoles actually displayed
for ii=tabl_seed_dip(:,1)
    set(sdip.hdl.hseed(ii),'BackgroundColor',[.7 1 .7]);
end

%==========================================================================
case 'clearall'
%==========================================================================
% Clears all dipoles, and reset the toggle buttons
%--------------------------------------------------------------------------

if isfield(st.vols{1},'sdip')
    sdip = st.vols{1}.sdip;
else
    error('I can''t find sdip structure');
end

disp('Clears all dipoles')
spm_eeg_inv_vbecd_disp('ClearDip');
for ii=1:st.vols{1}.sdip.n_seeds
    if sdip.exitflag(ii)==1
        set(st.vols{1}.sdip.hdl.hseed(ii),'Value',0);
    end
end
set(st.vols{1}.sdip.hdl.hdip,'Value',1);

%==========================================================================
case 'chgseed'
%==========================================================================
% Changes the seeds displayed
%--------------------------------------------------------------------------

% disp('Change seed')

sdip = st.vols{1}.sdip;
if isfield(sdip,'tabl_seed_dip')
    prev_seeds = p_seed(sdip.tabl_seed_dip);
else
    prev_seeds = [];
end

l_seed = zeros(sdip.n_seeds,1);
for ii=1:sdip.n_seeds
    if sdip.exitflag(ii)==1
        l_seed(ii) = get(sdip.hdl.hseed(ii),'Value');
    end
end
l_seed = find(l_seed);

% Modify the list of seeds displayed
if isempty(l_seed)
    % Nothing left displayed
    i_seed=[];
elseif isempty(prev_seeds)
    % Just one dipole added, nothing before
    i_seed=l_seed;
elseif length(prev_seeds)>length(l_seed)
    % One seed removed
    i_seed = prev_seeds;
    for ii=1:length(l_seed)
        p = find(prev_seeds==l_seed(ii));
        if ~isempty(p)
            prev_seeds(p) = [];
        end % prev_seeds is left with the index of the one removed
    end
    i_seed(i_seed==prev_seeds) = [];
    % Remove the dipole & change the button colour
    spm_eeg_inv_vbecd_disp('ClearDip',prev_seeds);
    set(sdip.hdl.hseed(prev_seeds),'BackgroundColor',[.7 .7 .7]);
else
    % One dipole added
    i_seed = prev_seeds;
    for ii=1:length(prev_seeds)
        p = find(prev_seeds(ii)==l_seed);
        if ~isempty(p)
            l_seed(p) = [];
        end % l_seed is left with the index of the one added
    end
    i_seed = [i_seed ; l_seed];
end

i_dip = get(sdip.hdl.hdip,'Value');
spm_eeg_inv_vbecd_disp('ClearDip');
if ~isempty(i_seed)
    spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip);
end

%==========================================================================
case 'chgdip'
%==========================================================================
% Changes the dipole index for the first seed displayed
%--------------------------------------------------------------------------

disp('Change dipole')
sdip = st.vols{1}.sdip;

i_dip = get(sdip.hdl.hdip,'Value');
if isfield(sdip,'tabl_seed_dip')
    i_seed = p_seed(sdip.tabl_seed_dip);
else
    i_seed = [];
end

if ~isempty(i_seed)
    spm_eeg_inv_vbecd_disp('ClearDip')
    spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip);
end

%==========================================================================
case 'cleardip'
%==========================================================================
% FORMAT spm_eeg_inv_vbecd_disp('ClearDip',seed_i)
% e.g. spm_eeg_inv_vbecd_disp('ClearDip')
%       clears all displayed dipoles
% e.g. spm_eeg_inv_vbecd_disp('ClearDip',1)
%       clears the first dipole displayed
%--------------------------------------------------------------------------

if nargin>2
    seed_i = varargin{1};
else
    seed_i = 0;
end
if isfield(st.vols{1},'sdip')
        sdip = st.vols{1}.sdip;
    else
        return; % I don't do anything, as I can't find sdip strucure
end

if isfield(sdip,'ax')
    Nax = size(sdip.ax,2);
    else
        return; % I don't do anything, as I can't find axes info
end

if seed_i==0 % removes everything
    for ii=1:Nax
        for jj=1:9
            delete(sdip.ax(jj,ii));
        end
    end
    for ii=sdip.tabl_seed_dip(:,1)
        set(sdip.hdl.hseed(ii),'BackgroundColor',[.7 .7 .7]);
    end
    sdip = rmfield(sdip,'tabl_seed_dip');
    sdip = rmfield(sdip,'ax');
elseif seed_i<=Nax % remove one seed only
    l_seed = find(sdip.tabl_seed_dip(:,1)==seed_i);
    for ii=l_seed
        for jj=1:9
            delete(sdip.ax(jj,ii));
        end
    end
    sdip.ax(:,l_seed) = [];
    sdip.tabl_seed_dip(l_seed,:) = [];
else
    error('Trying to clear unspecified dipole');
end
st.vols{1}.sdip = sdip;

%==========================================================================
case 'redrawdip'
%==========================================================================
% spm_eeg_inv_vbecd_disp('RedrawDip')
% redraw everything, useful when zooming into image
%--------------------------------------------------------------------------

% spm_eeg_inv_vbecd_disp('ClearDip')
% spm_eeg_inv_vbecd_disp('ChgDip')

% disp('Change dipole')
sdip = st.vols{1}.sdip;

i_dip = get(sdip.hdl.hdip,'Value');
if isfield(sdip,'tabl_seed_dip')
    i_seed = p_seed(sdip.tabl_seed_dip);
else
    i_seed = [];
end

if ~isempty(i_seed)
    spm_eeg_inv_vbecd_disp('ClearDip')
    spm_eeg_inv_vbecd_disp('DrawDip',i_seed,i_dip);
end

%==========================================================================
otherwise
%==========================================================================
    warning('Unknown action string');
end
% warning(sw);
return

%==========================================================================
% dh = add1dip(loc,js,vloc,mark,col,ax,Fig,bb)
%==========================================================================
function dh = add1dip(loc,js,vloc,mark,col,ax,Fig,bb)
% Plots the dipoles on the 3 views, with an error ellipse for location
% Then returns the handle to the plots

global st

is = inv(st.Space);
loc = is(1:3,1:3)*loc(:) + is(1:3,4);
% taking into account the zooming/scaling only for the location
% NOT for the dipole's amplitude.
% Amplitude plotting is quite arbitrary anyway and up to some scaling
% defined for better viewing...

loc(1,:) = loc(1,:) - bb(1,1)+1;
loc(2,:) = loc(2,:) - bb(1,2)+1;
loc(3,:) = loc(3,:) - bb(1,3)+1;
% +1 added to be like John's orthview code

% prepare error ellipse
vloc = is(1:3,1:3)*vloc*is(1:3,1:3);
[V,E] = eig(vloc);
VE = V*diag(sqrt(diag(E))); % use std
% VE = V*E;   % or use variance ???

dh = zeros(9,1);
figure(Fig)
% Transverse slice, # 1
%----------------------
set(Fig,'CurrentAxes',ax{1}.ax)
set(ax{1}.ax,'NextPlot','add')
dh(1) = plot(loc(1),loc(2),[mark,col],'LineWidth',1);
dh(2) = plot(loc(1)+[0 js(1)],loc(2)+[0 js(2)],col,'LineWidth',2);
% add error ellipse
[uu,ss,vv] = svd(VE([1 2],:));
[phi] = cart2pol(uu(1,1),uu(2,1));
e = diag(ss);
t = (-1:.02:1)*pi;
x = e(1)*cos(t)*cos(phi)-e(2)*sin(t)*sin(phi)+loc(1);
y = e(2)*sin(t)*cos(phi)+e(1)*cos(t)*sin(phi)+loc(2);
dh(3) = plot(x,y,[':',col],'LineWidth',.5);
set(ax{1}.ax,'NextPlot','replace')

% Coronal slice, # 2
%----------------------
set(Fig,'CurrentAxes',ax{2}.ax)
set(ax{2}.ax,'NextPlot','add')
dh(4) = plot(loc(1),loc(3),[mark,col],'LineWidth',1);
dh(5) = plot(loc(1)+[0 js(1)],loc(3)+[0 js(3)],col,'LineWidth',2);
% add error ellipse
[uu,ss,vv] = svd(VE([1 3],:));
[phi] = cart2pol(uu(1,1),uu(2,1));
e = diag(ss);
t = (-1:.02:1)*pi;
x = e(1)*cos(t)*cos(phi)-e(2)*sin(t)*sin(phi)+loc(1);
y = e(2)*sin(t)*cos(phi)+e(1)*cos(t)*sin(phi)+loc(3);
dh(6) = plot(x,y,[':',col],'LineWidth',.5);
set(ax{2}.ax,'NextPlot','replace')

% Sagital slice, # 3
%----------------------
set(Fig,'CurrentAxes',ax{3}.ax)
set(ax{3}.ax,'NextPlot','add')
% dh(5) = plot(dim(2)-loc(2),loc(3),[mark,col],'LineWidth',2);
% dh(6) = plot(dim(2)-loc(2)+[0 -js(2)],loc(3)+[0 js(3)],col,'LineWidth',2);
dh(7) = plot(bb(2,2)-bb(1,2)-loc(2),loc(3),[mark,col],'LineWidth',1);
dh(8) = plot(bb(2,2)-bb(1,2)-loc(2)+[0 -js(2)],loc(3)+[0 js(3)],col,'LineWidth',2);
% add error ellipse
[uu,ss,vv] = svd(VE([2 3],:));
[phi] = cart2pol(uu(1,1),uu(2,1));
e = diag(ss);
t = (-1:.02:1)*pi;
x = -(e(1)*cos(t)*cos(phi)-e(2)*sin(t)*sin(phi))+bb(2,2)-bb(1,2)-loc(2);
y = e(2)*sin(t)*cos(phi)+e(1)*cos(t)*sin(phi)+loc(3);
dh(9) = plot(x,y,[':',col],'LineWidth',.5);
set(ax{3}.ax,'NextPlot','replace')

return

%==========================================================================
% pr_seed = p_seed(tabl_seed_dip)
%==========================================================================
function pr_seed = p_seed(tabl_seed_dip)
% Gets the list of seeds used in the previous display

ls = sort(tabl_seed_dip(:,1));
if length(ls)==1
    pr_seed = ls;
else
    pr_seed = ls([find(diff(ls)) ; length(ls)]);
end

%

% OLD STUFF

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use it with arguments or not:
% - spm_eeg_inv_vbecd_disp('Init')
%        The routine asks for the dipoles file and image to display
% - spm_eeg_inv_vbecd_disp('Init',sdip)
%        The routine will use the avg152T1 canonical image
% - spm_eeg_inv_vbecd_disp('Init',sdip,P)
%        The routines dispays the dipoles on image P.
%
% If multiple seeds have been used, you can select the seeds to display
% by pressing their index.
% Given that the sources could have different locations, the slices
% displayed will be the 3D view at the *average* or *mean* locations of
% selected sources.
% If more than 1 dipole was fitted at a time, then selection of source 1
% to N is possible through the pull-down selector.
%
% The location of the source/cut is displayed in mm and voxel, as well as
% the underlying image intensity at that location.
% The cross hair position can be hidden by clicking on its button.
%
% Nota_1: If the cross hair is manually moved by clicking in the image or
%       changing its coordinates, the dipole displayed will NOT be at
%       the right displayed location. That's something that needs to be improved...
%
% Nota_2: Some seeds may have not converged within the limits fixed,
%       these dipoles are not displayed...
%
% Fields needed in sdip structure to plot on an image:
%       + n_seeds: nr of seeds set used, i.e. nr of solutions calculated
%       + n_dip: nr of fitted dipoles on the EEG time series
%       + loc: location of fitted dipoles, cell{1,n_seeds}(3 x n_dip)
%               remember that loc is fixed over the time window.
%       + j: sources amplitude over the time window,
%            cell{1,n_seeds}(3*n_dip x Ntimebins)
%       + Mtb: index of maximum power in EEG time series used

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % First point to consider
% loc_mm = sdip.loc{i_seed(1)}(:,i_dip);
%
% % PLace the underlying image at right cuts
% spm_orthviews('Reposition',loc_mm);
% % spm_orthviews('Reposition',loc_vx);
% % spm_orthviews('Xhairs','off')
%
% % if i_seed = set, Are there other dipoles close enough ?
% tabl_seed_dip=[i_seed(1) i_dip]; % table summarising which set & dip to use.
% if length(i_seed)>1
%   unit = ones(1,sdip.n_dip);
%   for ii = i_seed(2:end)'
%         d2 = sqrt(sum((sdip.loc{ii}-loc_mm*unit).^2));
%         l_cl = find(d2<=lim_cl);
%         if ~isempty(l_cl)
%             for jj=l_cl
%                 tabl_seed_dip = [tabl_seed_dip ; [ii jj]];
%             end
%         end
%   end
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get(sdip.hdl.hseed(1),'Value')
% for ii=1:sdip.n_seeds, delete(hseed(ii)); end
% h1 = uicontrol(Fig,'Style','togglebutton','Position',[600 25 10 10].*WS)
% h2 = uicontrol(Fig,'Style','togglebutton','Position',[620 100 20 20].*WS,'String','1')
% h2 = uicontrol(Fig,'Style','checkbox','Position',[600 100 10 10].*WS)
% h3 = uicontrol(Fig,'Style','radiobutton','Position',[600 150 20 20].*WS)
% h4 = uicontrol(Fig,'Style','radiobutton','Position',[700 150 20 20].*WS)
% delete(h2),delete(h3),delete(h4),
% delete(hdip)

author: spm


text:
function F = filmStrip( I, overlap, delta, border )
% Used to display R stacks of T images as a "filmstrip".
%
% See examples below to see what is meant by "filmstrip".
%
% USAGE
%  F = filmStrip( I, overlap, delta, border )
%
% INPUTS
%  I          - MxNxTxR or MxNx1xTxR or MxNx3xTxR array
%               (of bw or color images). R can equal 1.
%  overlap    - amount of overlap between successive frames
%  delta      - amount to shift each successive frame upward
%  border     - width of black border around each frame
%
% OUTPUTS
%  F       - filmstrip
%
% EXAMPLE - one filmstrip
%  load images;
%  F1 = filmStrip( video(:,:,1:15), 10, 2, 5 );   figure(1); im(F1); % one
%  F2 = filmStrip( videos(:,:,:,1:10), 5, 2, 3 ); figure(2); im(F2); % many
%
% See also MONTAGE2
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% convert I to be of type double, and have dimensions MxNxCxTxR
I = double(I); I = I/max(I(:)); sizI = size(I);
if(~any(sizI(3)==[1 3])); I=reshape(I,[sizI(1:2),1,sizI(3:end)]); end
I = padarray( I, [border border 0 0 0], 0, 'both' );
[mRows, nCols, nColor, nFrame, nStrip] = size(I);

% size of final filmstip object
sizF1 = [mRows+delta*(nFrame-1), nFrame*nCols-overlap*(nFrame-1), nColor];
sizF=sizF1;  sizF(1)=sizF1(1)*nStrip - ((nFrame-1-2)*delta)*(nStrip-1);
mRowsF1 = sizF1(1);

for i=1:nStrip

  % Create i-th film strip
  Fi = -ones( sizF1 );  row = 1;  col = sizF1(2);
  for f=nFrame:-1:1
    Fi( row:(row+mRows-1), (col-nCols+1):col, : ) = I(:,:,:,f,i);
    row = row + delta;  col = col - nCols + overlap;
  end

  % stop if creating single filmstrip
  if( nStrip==1 ); F=Fi; break; end

  % merge with the previous film strips
  if( i==1 ); F = -ones( sizF );  row2=1;  end
  Fc = F( row2:(row2+mRowsF1-1), : );
  locs=(Fc<0);  Fc(locs) = Fi(locs);
  F( row2:(row2+mRowsF1-1), :  ) = Fc;
  row2 = row2 + mRowsF1 - ((nFrame-1-2)*delta);
end
F(F<0)=1;

author: zouchuhang

repo: LayoutNet


text:
function [ b ] = polyshift( a, shift )
%POLYSHIFT Compute coefficients of a shifted 1-D polynomial
%
% The input A is a polynomial of length n and represents the polynomial of
% order (n-1):
% a(1) + a(2)*x + a(3)*(x^2) + ... + a(n)*(x^(n-1))
%
% The polynomial shifted by the scalar input SHIFT is:
% a(1) + a(2)*(x+shift) + a(3)*((x+shift)^2) + ... + a(n)*((x+shift)^(n-1))
% =
% b(1) + b(2)*x + b(3)*(x^2) + ... + b(n)*(x^(n-1))
%
% This function is a direct computation of the coefficients B of the
% shifted polynomial.  "Synthetic division" would be more efficient, but it
% doesn't matter much when only computing a few shifts for small
% polynomials.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

b=zeros(size(a));
for j = 1:numel(a)
    for k=j:numel(a)
        b(j) = b(j) + (a(k)*nchoosek(k-1,j-1)*shift^(k-j));
    end
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function value = i4_huge ( )

%*****************************************************************************80
%
%% I4_HUGE returns a "huge" I4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, a huge integer.
%
  value = 2147483647;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function segment = p03_boundary_segment ( segment_index, m, segment_length )

%*****************************************************************************80
%
%% P03_BOUNDARY_SEGMENT returns a boundary segment in problem 03.
%
%  Discussion:
%
%    For boundary segment #1, the value of SEGMENT_LENGTH should be
%    at least 5.  Values of 4*N+1 will result in an "even" mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer SEGMENT_INDEX, the index of the boundary segment.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer SEGMENT_LENGTH, the number of points in the segment.
%
%    Output, SEGMENT(M,SEGMENT_LENGTH), points on the boundary segment.
%
  center = [ 0.0, 0.0 ];
  r2 = 0.4;

  if ( segment_index == 1 )

    n1 = round ( ( segment_length - 1 ) / 4 );
    n2 = round ( ( 2 * ( segment_length - 1 ) ) / 4 ) - n1;
    n3 = round ( ( 3 * ( segment_length - 1 ) ) / 4 ) - n1 - n2;
    n4 = segment_length - 1 - n1 - n2 - n3;

    s(1:2,1) = [ -1.0, -1.0 ]';
    s(1:2,2) = [  1.0, -1.0 ]';
    s(1:2,3) = [  1.0,  1.0 ]';
    s(1:2,4) = [ -1.0,  1.0 ]';

    j = 0;

    for i = 1 : n1
      j = j + 1;
      segment(1:2,j) = ( ( n1 - i + 1 ) * s(1:2,1)   ...
                       + (      i - 1 ) * s(1:2,2) ) ...
                       / ( n1         );
    end

    for i = 1 : n2
      j = j + 1;
      segment(1:2,j) = ( ( n2 - i + 1 ) * s(1:2,2)   ...
                       + (      i - 1 ) * s(1:2,3) ) ...
                       / ( n2         );
    end

    for i = 1 : n3
      j = j + 1;
      segment(1:2,j) = ( ( n3 - i + 1 ) * s(1:2,3)   ...
                       + (      i - 1 ) * s(1:2,4) ) ...
                       / ( n3         );
    end

    for i = 1 : n4
      j = j + 1;
      segment(1:2,j) = ( ( n4 - i + 1 ) * s(1:2,4)   ...
                       + (      i - 1 ) * s(1:2,1) ) ...
                       / ( n4         );
    end

    j = j + 1;
    segment(1:2,j) = s(1:2,1);

  elseif ( segment_index == 2 )

    for i = 1 : segment_length
      angle = 2.0 * pi * ( segment_length - i ) / ( segment_length - 1 );
      segment(1,i) = center(1) + r2 * cos ( angle );
      segment(2,i) = center(2) + r2 * sin ( angle );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P03_BOUNDARY_SEGMENT - Fatal error!\n' );
    fprintf ( 1, '  Illegal SEGMENT_INDEX = %d\n', segment_index );
    error ( 'P03_BOUNDARY_SEGMENT - Fatal error!' );

  end

  return
end