text:
% pdfbtrainthmt.m
% written by: Duncan Po
% Date: August 24, 2002
% Using the EM algorithm, train a model for the provided data
% Usage: model = pdfbtrainthmt(tree, levndir, mD, ns,zeromean)
%        [model, stateprob] = pdfbtrainthmt(tree, levndir, mD, ns,zeromean)
%        model = pdfbtrainthmt(tree, levndir, mD, initmodel)
%        [model, stateprob] = pdfbtrainthmt(tree, levndir, mD, initmodel)
% Inputs:   tree        - data in tree structure
%           levndir     - the number of subbands in each level (e.g. [2 2 3 3])
%           mD          - convergence value
%           ns          - number of states in the desired model
%           zeromean    - 'yes' for zero mean model and 'no' for non-zeromean model
%           initmodel   - an initial model can be provided to speed up the training
% Output:   model       - the model generated
%           stateprob   - state probabilities

function [model, stateprob] = pdfbtrainthmt(tree, levndir, mD, ns,zeromean)

if nargout == 2
    needstateprob = 1;
else
    needstateprob = 0;
end;

nlevel = length(tree);

if nargin == 4
    initmodel.nstates = ns.nstates;
    realns = initmodel.nstates;
    initmodel.nlevels = ns.nlevels;
    initmodel.zeromean = ns.zeromean;
    initmodel.rootprob = zeros(1,realns);
    for l = 1:realns
        initmodel.rootprob(l) = ns.rootprob(l);
    end;
    for l = 1:nlevel-1
        for k = 1:2.^(levndir(l+1)-levndir(1))
            initmodel.transprob{l}{k} = zeros(realns);
            for m = 1:realns
                for n = 1:realns
                    initmodel.transprob{l}{k}(m,n) = ns.transprob{l}{k}(m,n);
                end;
            end;
        end;
    end;
    for l = 1:nlevel
        for k = 1:2.^(levndir(l)-levndir(1))
            if strcmp(ns.zeromean, 'yes') == 0
                initmodel.mean{l}{k} = zeros(1,realns);
                for m = 1:realns
                    initmodel.mean{l}{k}(m) = ns.mean{l}{k}(m);
                end;
            end;
            initmodel.stdv{l}{k} = zeros(1, realns);
            for m = 1:realns
                initmodel.stdv{l}{k}(m) = ns.stdv{l}{k}(m);
            end;
        end;
    end;
else
    model.nstates = -2;
    model.nlevels = -1;
    model.zeromean = 0;
    model.rootprob = zeros(1,ns);
    for l = 1:nlevel-1
        for k = 1:2.^(levndir(l+1)-levndir(1))
            model.transprob{l}{k} = zeros(ns);
        end;
    end;
    for l = 1:nlevel
        for k = 1:2.^(levndir(l)-levndir(1))
            if strcmp(zeromean, 'yes') == 0
                model.mean{l}{k} = zeros(1,ns);
            end;
            model.stdv{l}{k} = zeros(1, ns);
        end;
    end;
end;

if ~exist('initmodel', 'var')
    if needstateprob == 0
        pdfbtrain_thmt(ns, nlevel, levndir, zeromean, tree, mD, model);
    else
        for l = 1:nlevel
            numofel = length(tree{l});
            stateprob{l} = zeros(numofel, ns);
        end;
        pdfbtrain_thmt(ns, nlevel, levndir, zeromean, tree, mD, model, stateprob);
    end;
else
    model = initmodel;
    if needstateprob == 0
        pdfbprotrain_thmt(ns, nlevel, levndir, tree ,mD, model);
    else
        ns = initmodel.nstates;
        for l = 1:nlevel
            numofel = length(tree{l});
            stateprob{l} = ones(numofel, ns);
        end;
        pdfbprotrain_thmt(ns, nlevel, levndir, tree, mD, model, stateprob);
    end;
end;


text:
function [ vX, mX ] = SolveLsL1Prox( mA, vB, paramLambda, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX, mX ] = SolveLsL1Prox( mA, vB, lambdaFctr, numIterations )
% Solve L1 Regularized Least Squares Using Proximal Gradient (PGM) Method.
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia PGM - https://en.wikipedia.org/wiki/Proximal_gradient_method.
% Remarks:
%   1.  Using vanilla PGM.
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     23/08/2017
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

mAA = mA.' * mA;
vAb = mA.' * vB;
vX  = pinv(mA) * vB; %<! Dealing with "Fat Matrix"

stepSize = 1 / (2 * (norm(mA, 2) ^ 2));
% stepSize = 1 / sum(mA(:) .^ 2); %<! Faster to calculate, conservative (Hence slower)

mX = zeros([size(vX, 1), numIterations]);
mX(:, 1) = vX;

for ii = 2:numIterations

    vG = (mAA * vX) - vAb;
    vX = ProxL1(vX - (stepSize * vG), stepSize * paramLambda);

    mX(:, ii) = vX;

end

end

function [ vX ] = ProxL1( vX, lambdaFactor )

% Soft Thresholding
vX = max(vX - lambdaFactor, 0) + min(vX + lambdaFactor, 0);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
classdef (Abstract) DisplayItem < handle
    % Abstract class for items in frost.Animator.Display
    %
    % @author omar @date 2017-06-01
    %
    % Copyright (c) 2017, UMICH Biped Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/BSD-3-Clause

    properties (SetAccess = private, GetAccess = protected)
        ax
        model
    end

    properties (SetAccess = private, GetAccess = public)
        name
    end

    methods
        function obj = DisplayItem(ax, model, name)
            obj.ax = ax;
            obj.model = model;
            obj.name = name;
        end
    end

    methods (Abstract)
        update(obj, x);
    end
end

author: ayonga

repo: frost-dev

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


text:
function bearpath = bearroot()
bearpath = fileparts(mfilename('fullpath'));
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [Fs spec] = PTSpec2d(Y, F, psd)
%% [Fs spec] = PTSpec2d(Y, F, psd)
%   computes standard fft on input data Y. F is sample frequency in Hz.

% ----------------------------------------------------------------------------------
% "THE BEER-WARE LICENSE" (Revision 42):
% <brian.white@queensu.ca> wrote this file. As long as you retain this notice you
% can do whatever you want with this stuff. If we meet some day, and you think
% this stuff is worth it, you can buy me a beer in return. -Brian White
% ----------------------------------------------------------------------------------

if psd
    % N = length(Y);
    % [psdx,Fs] = periodogram(Y,[], N-1, F*1000,'psd'); % power psd
    % [spec,Fs] = pspectrum(Y, F*1000, 'FrequencyResolution', 10);

    N = length(Y);
    Fs = ((F*1000)*(0:(N/2))/N);
    Y = Y.*hann(N)';
    Y = fft(Y);
    psdx = abs(Y).^2 / (F*1000*N); % (1/(F*N)) * abs(Y).^2  is exactly same as  abs(Y).^2 / (F*N), to
    psdx(2:end-1) = 2*psdx(2:end-1);
    psdx = psdx(1:N/2+1);
%     % scale to dB
     spec = 10 * log10(psdx)';
else
    N = length(Y);
    Fs = ((F*1000)*(0:(N/2))/N);
    Y = Y.*hann(N)';
    Y = fft(Y);
    spec = abs(Y) / N; % (1/(F*N)) * abs(Y).^2  is exactly same as  abs(Y).^2 / (F*N), to
    spec = spec(1:N/2+1)';
end

end

repo: PIDtoolbox


text:
function aedat = RescaleTime(aedat, factor, packetTimeStamps, dataTimeStamps)

%{
Stretches or squeezes time by multiplying all timestamps by a factor.
This affects both data.<dataType>.timeStamp and info.packetTimeStamps,
by default, overridden by the optional flags.
%}

dbstop if error

if ~exist('packetTimeStamps', 'var') || packetTimeStamps

    if ~isfield(aedat, 'info') || ~isfield(aedat.info, 'packetTimeStamps')
        disp('No packet timestamps found')
    else
        aedat.info.packetTimeStamps = uint64(double(aedat.info.packetTimeStamps * factor));
    end
end

if ~exist('dataTimeStamps', 'var') || dataTimeStamps

    if ~isfield(aedat, 'data')
        disp('No data found')
    else

        % Special
        if isfield(aedat.data, 'special')
            aedat.data.special.timeStamp = uint64(double(aedat.data.special.timeStamp * factor));
        end

        % Polarity
        if isfield(aedat.data, 'polarity')
            aedat.data.polarity.timeStamp = uint64(double(aedat.data.polarity.timeStamp * factor));
        end

        % Frames
        % This assumes that timestamps have been simplified to aedat2 standard, if
        % they came from aedat3 file
        if isfield(aedat.data, 'frame')
            aedat.data.frame.timeStampStart = uint64(double(aedat.data.frame.timeStampStart * factor));
            aedat.data.frame.timeStampEnd = uint64(double(aedat.data.frame.timeStampEnd * factor));
        end

        % Imu6
        if isfield(aedat.data, 'imu6')
            aedat.data.imu6.timeStamp = uint64(double(aedat.data.imu6.timeStamp * factor));
        end

        if isfield(aedat.data, 'sample')
            aedat.data.sample.timeStamp = uint64(double(aedat.data.sample.timeStamp * factor));
        end

        if isfield(aedat.data, 'ear')
            aedat.data.ear.timeStamp = uint64(double(aedat.data.ear.timeStamp * factor));
        end

        if isfield(aedat.data, 'point1D')
            aedat.data.point1D.timeStamp = uint64(double(aedat.data.point1D.timeStamp * factor));
        end

        if isfield(aedat.data, 'point2D')
            aedat.data.point2D.timeStamp = uint64(double(aedat.data.point2D.timeStamp * factor));
        end

        if isfield(aedat.data, 'point3D')
            aedat.data.point3D.timeStamp = uint64(double(aedat.data.point3D.timeStamp * factor));
        end

        aedat = FindFirstAndLastTimeStamps(aedat);

        if isfield(aedat.info, 'packetTimeStamps')
            aedat.info.packetTimeStamps = uint64(double(aedat.info.packetTimeStamps * factor));
        end
    end
end

author: panpanfei


text:
function test04 ( dim_num, level_max )

%*****************************************************************************80
%
%% TEST04 sums the weights and compares them to 2^DIM_NUM.
%
%  Discussion:
%
%    This routine gets the sparse grid indices and determines the
%    corresponding sparse grid abscissas.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the level.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04:\n' );
  fprintf ( 1, '  Compute the weights of a Clenshaw Curtis sparse grid .\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  As a simple test, sum these weights.\n' );
  fprintf ( 1, '  They should sum to exactly 2^DIM_NUM.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
%
%  Determine the number of points.
%
  point_num = sparse_grid_cfn_size ( dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of unique points in the grid = %d\n', point_num );
%
%  Compute the weights and points.
%
  [ grid_weight, grid_point ] = sparse_grid_cc ( dim_num, level_max, point_num );
%
%  Sum the weights.
%
  weight_sum = sum ( grid_weight(1:point_num) );

  weight_sum_exact = 2.0^dim_num;

  weight_sum_error = abs ( weight_sum - weight_sum_exact );

  fprintf ( 1, '\n' ) ;
  fprintf ( 1, '    Weight sum     Exact sum    Difference\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %12e  %12e  %12e\n', ...
    weight_sum, weight_sum_exact, weight_sum_error );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function nrm = norm(X)
%NORM Norm of a ttensor.
%
%   See also TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if prod(size(X)) > prod(size(X.core))
    V = cell(ndims(X),1);
    for n = 1:ndims(X)
        V{n} = X.u{n}'*X.u{n};
    end
    Y = ttm(X.core,V);
    tmp = innerprod(Y, X.core);
    nrm = sqrt(tmp);
else
    nrm = norm(full(X));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function varargout = daqstopbutton(fig,obj,varargin);
%DAQSTOPBUTTON        Add a stop/start button to a data acquisition application
%
% DAQSTOPBUTTON(FIG,OBJ) adds a start/stop button to figure FIG.  This
% button can be used to start and stop data acquisition object OBJ.
% DAQSTOPBBUTTON will also delete OBJ when FIG is closed (i.e., it sets
% FIG's CloseRequestFcn to delete the object)
%
% DAQSTOPBUTTON(FIG,OBJ,'P1','V1','P2','V2', ...) specifies Property-Value
% pairs for configuring properties of the start/stop button.  Any valid
% property of a togglebutton can be specified.
%
% HBUTTON = DAQSTOPBUTTON(...) returns a handle to the start/stop button
%
% Note: This button does not "listen" to determine if your object changes
%   state.  For instance, if you have a finite number of samples
%   acquired after starting, the button will not reset automatically
%   when your object stops running.  Don't despair, since it is easy
%   enough to do on your own!
%
% Example:
%      fh = figure;                         % Create a figure
%      ai = analoginput('winsound');        % Create an input object
%      addchannel(ai,1);                    % Add a channel
%      set(ai,'TriggerRepeat',inf);         % Configure to run infinitely
%      set(ai,'TimerFcn','plot(peekdata(ai,500))'); % Each timer event will plot recent data
%      hButton = daqstopbutton(fh,ai);      % Add the stopbutton

%   Scott Hirsch
%   shirsch@mathworks.com
%   Copyright 2003 The MathWorks, Inc

%Error checking
msg = nargchk(2,inf,nargin);
error(msg)

if ~all(isvalid(obj))
    error('Second input argument must be valid daq object')
end;

% Create the button
hButton = uicontrol(fig,'style','togglebutton',varargin{:});

%Check current state of the object
val = strcmp(obj.Running,'On');

if val==1       %Already running
    set(hButton,'String','Stop');
    set(hButton,'Value',1)
else
    set(hButton,'String','Start');
    set(hButton,'Value',0)
end;

set(hButton,'Callback',{@localStartStopObject,obj})

% Configure the CloseRequestFcn of the figure
setappdata(fig,'DaqStopButtonObject',obj)
cr = get(fig,'CloseRequestFcn');
cr = ['obj=getappdata(gcf,''DaqStopButtonObject'');stop(obj);delete(obj);' cr];
set(fig,'CloseRequestFcn',cr);

% Return a handle to the button
if nargout
    varargout{1} = hButton;
end;

function localStartStopObject(hButton,action,obj)
% Callback for the start/stop button
val = get(hButton,'Value');
if val==1       %Pushed in
    set(hButton,'String','Stop');
    if all(isvalid(obj))
        start(obj)
    end;
else
    set(hButton,'String','Start');
    if all(isvalid(obj))
        stop(obj)
    end;
end;

author: Sable

repo: mcbench-benchmarks


text:
% Kernel Affine Projection (KAP) algorithm with coherence criterion
%
% C. Richard, J.C.M. Bermudez, and P. Honeine, "Online Prediction of Time
% Series Data With Kernels," IEEE Transactions on Signal Processing,
% vol. 57, no. 3, pp. 1058-1067, March 2009,
% http://dx.doi.org/10.1109/TSP.2008.2009895
%
% Remark: memories are initialized empty in this implementation
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef kap < kernel_adaptive_filter

    properties (GetAccess = 'public', SetAccess = 'private')
        mu0 = .95; % coherence criterion threshold
        eta = .5; % step size
        eps = 1E-2; % regularization
        p = 20; % memory length
        kerneltype = 'gauss'; % kernel type
        kernelpar = 1; % kernel parameter
    end

    properties (GetAccess = 'public', SetAccess = 'private')
        memx = []; % input memory
        memy = []; % output memory
        dict = []; % dictionary
        modict = []; % modulus of the dictionary elements
        alpha = []; % expansion coefficients
    end

    methods
        function kaf = kap(parameters) % constructor
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(parameters)'
                    if ismember(fn,fieldnames(kaf))
                        kaf.(fn{1}) = parameters.(fn{1});
                    end
                end
            end
        end

        function y_est = evaluate(kaf,x) % evaluate the algorithm
            if size(kaf.dict,1)>0
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                y_est = k'*kaf.alpha;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(kaf,x,y) % train the algorithm
            if (length(kaf.memy) < kaf.p)
                kaf.memx = [kaf.memx; x]; % grow the memory
                kaf.memy = [kaf.memy; y]; % grow the memory
            else
                kaf.memx = [kaf.memx(2:end,:); x]; % sliding memory
                kaf.memy = [kaf.memy(2:end); y]; % sliding memory
            end

            if size(kaf.dict,2)==0 % initialize
                k = kernel(x,x,kaf.kerneltype,kaf.kernelpar);
                kaf.dict = x;
                kaf.modict = sqrt(k);
                kaf.alpha = 0;
            else
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                kx = kernel(x,x,kaf.kerneltype,kaf.kernelpar);
                C = k./(sqrt(kx)*kaf.modict); % coherence
                if (max(C) <= kaf.mu0) % coherence criterion
                    kaf.dict = [kaf.dict; x]; % order increase
                    kaf.alpha = [kaf.alpha; 0]; % order increase
                end
            end

            H = kernel(kaf.memx,kaf.dict,kaf.kerneltype,kaf.kernelpar);
            kaf.alpha = kaf.alpha + ...
                kaf.eta*H'/...
                (kaf.eps*eye(size(H,1)) + H*H')*...
                (kaf.memy - H*kaf.alpha);
        end
    end
end

repo: kafbox


text:
function trans = homothecy(point, ratio)
%HOMOTHECY create a homothecy as an affine transform
%
%   TRANS = homothecy(POINT, K);
%   POINT is the center of the homothecy, K is its factor.
%
%   See also:
%   transforms2d, transformPoint, createTranslation
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 20/01/2005.
%

%   HISTORY
%   22/04/2009: copy to createHomothecy and deprecate

% deprecation warning
warning('geom2d:deprecated', ...
    '''homothecy'' is deprecated, use ''createHomothecy'' instead');

% call current implementation
trans = createHomothecy(point, ratio);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% axcopy() - Copy a Matlab figure axis and its graphic objects to a new pop-up window
%            using the left mouse button.
%
% Usage:  >> axcopy
%         >> axcopy(fig)
%         >> axcopy('noticks')
%         >> axcopy(fig, command)
%
% Notes:
%   1) Clicking the left mouse button on a Matlab figure axis copies the graphic objects
%      in the axis to a new (pop-up) figure window.
%   2) Option 'noticks' does not make x and y tickloabelmodes 'auto' in the pop-up.
%   2) The command option is an optional string that is evaluated in the new window
%      when it pops up. This allows the user to customize the pop-up display.
%   3) Deleting the pop-up window containing the copied axis leaves the selected axis
%      as the current graphic axis (gca).

% Authors: Tzyy-Ping Jung & Scott Makeig, SCCN/INC/UCSD, La Jolla, 2000

% Copyright (C) 2000 Tzyy-Ping Jung & Scott Makeig, SCCN/INC/UCSD, 3-16-00
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% requires copyaxes.m
% 4-2-00 added 'noticks' -sm
% 01-25-02 reformated help & license -ad
% 02-16-02 debuged exist('fig') -> exist('fig') == 1 -ad
% 02-16-02 added the ignore parent size option -ad
% 03-11-02 remove size option and added command option -ad

function axcopy_modified(fig, command)
if (exist('fig') == 1) && strcmp(fig,'noticks')
   noticks = 1;
   if nargin> 1
     shift
   else
     fig = [];
   end
end
if ~(exist('fig') ==1) || isempty(fig) || fig == 0
   fig = gcf;
end
if ~strcmpi(get(fig, 'type'), 'axes')
    hndl= findobj('parent',fig,'type','axes');
else
    hndl = fig;
end;
offidx = [];
if exist('command') ~= 1
    comstr = 'copyaxis_modified';
else
   command_dbl = double(command);
   comstr = double(['copyaxis_modified(''' char(command_dbl) ''')']);
end;
for a=1:length(hndl)             % make all axes visible
    allobjs = findobj('parent',hndl(a));
    for index = 1:length(allobjs)
        if isempty(get(allobjs(index), 'ButtonDownFcn'))
            set(allobjs(index), 'ButtonDownFcn', char(comstr));
        end;
    end;
end
if ~strcmpi(get(fig, 'type'), 'axes')
    figure(fig);
else
    figure(get(fig, 'parent'));
end;
if exist('command') ~= 1
    set(hndl(a),'ButtonDownFcn','copyaxis_modified');
else
    % set(hndl,'ButtonDownFcn',['copyaxis(''' command ''')']);
    comstr = double(['copyaxis_modified(''' char(command_dbl) ''')']);
    set(hndl,'ButtonDownFcn',char(comstr));
end;
%set(hndl,'ButtonDownFcn','copyaxis');
%if ~exist('noticks')
  %axis on
  %set(gca,'xticklabelmode','auto')
  %set(gca,'yticklabelmode','auto')
%end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%% State Space MPC Tutorial
% This document explains how to use the setup function and online
% controller returned by ssmpcsetup.

%% State-space MPC set-up
% The online controller has to be set-up before use. To set-up the state
% space MPC controller, the user has to supply the state-space model,
% represented by A, B, C and D matrices (Note: it should be in
% discrete-time), the predication horizon, P, moving horizon, M and
% performance weights, Q and R. The default initial state and input are set
% to zero.
%
% SSMPC=MPCSETUP(A,B,C,D,P,M,Q,R,X0,U0);

%% Online controller
% The returned function handle from the MPC setup program is an online
% controller, SSMPC. The controller is called by proving two inputs:
% current measurement, Y and future reference, Ref. On return, it produces
% the optimal input, U for next step:
%
% U = SSMPC(Y,Ref);

%% A two-CSTR example
% A two-CSTR (Continuous Stired Reaction Tank) process is shown as follows.
%
% <html>
% <img src="2cstrPlant.png" width="600" height="400">
% </html>
%
% A linear state space of the model is developed for the plant. The model
% has six states, two inputs (two cooling water flow rates) and two
% disturbances (cooling water temperatures). The measured outputs are two
% tank temperatures. Further details of the process can be find in
% 1. Cao, Y and Yang, ZJ, "Multiobjective process controllability
% analysis", _Computers and Chemical Engineering_ , 28(2004), 83--90.
% 2. Al Seyab, RK and Cao, Y, "Differential recurrent neural network based
% predictive control", _Computers and Chemical Engineering_ , to appear.
%
% The discrete model with sampling rate 0.1 s is as follows:
A=[ 0.1555  -13.7665   -0.0604         0         0         0
    0.0010    1.0008    0.0068         0         0         0
         0    0.0374    0.9232         0         0         0
    0.0015   -0.1024   -0.0003    0.1587  -13.6705   -0.0506
         0    0.0061         0    0.0006    0.9929    0.0057
         0    0.0001         0         0    0.0366    0.9398];
Bu=[0.0001       0
         0       0
   -0.0036       0
         0  0.0001
         0       0
         0 -0.0028];
Bd=[      0         0
          0         0
     0.0013         0
          0         0
          0         0
          0    0.0008];
C=[0 362.995 0 0 0 0
   0 0 0 0 362.995 0];
D=zeros(2,2);

%% MPC parameters
% The MPC controller is configured with following parameters.
% Prediction horizon and moving horizon
p=10;
m=3;
% Performance wights
Q=1.5*eye(2*p);
R=eye(2*m);

%% MPC set-up
% The MPC controller is set-up by calling SSMPCSETUP:
ssmpc=mpcsetup(A,Bu,C,D,p,m,Q,R);

%% Simulation
% 150 seconds (1500 sampling intervals) simulation is conducted with
% several setpoint changes and random cooling water temperature changes
% within positive and negative 1 degree.
% Simulation length and variables for results
N=1500;
x0=zeros(6,1);
Y=zeros(N,2);
U=zeros(N,2);
% Predefined reference
T=zeros(N,2);
T(10:N,1)=1;
T(351:N,1)=3;
T(600:N,1)=5;
T(1100:N,1)=3;
T(100:N,2)=2;
T(451:N,2)=1;
T(700:N,2)=3;
T(1200:N,2)=4;
% Simulation
%%
for k=1:N
    % Process disturbances
    w=Bd*(rand(2,1)-0.5)*2;
    % Measurements noise
    v=0.01*randn(2,1);
    % actual measurement
    y=C*x0+v;
    % online controller
    u=ssmpc(y,T(k:end,:)');
    % plant update
    x0=A*x0+Bu*u+w;
    % save results
    Y(k,:)=y';
    U(k,:)=u';
end

%% Results
% The simulation results are summarized in two sub-plots.
t=(0:N-1)*0.1;
subplot(211)
plot(t,Y,t,T,'r--','linewidth',2)
title('output and setpoint')
ylabel('temp, C^\circ')
legend('T_1','T_2','Ref','location','southeast')
subplot(212)
stairs(t,U,'linewidth',2)
legend('u_1','u_2','location','southeast')
title('input')
ylabel('flow rate, m^3/s')
xlabel('time, s')


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% show LK tracking results
function LK_show(im1, im2, xFI, BB1, xFJ, BB2)

% Create a new image showing the two images side by side.
im3 = appendimages(im1,im2);

% Show a figure with lines joining the accepted matches.
h = figure(2);
set(h, 'Position', [100 100 2*size(im3,2) 2*size(im3,1)]);
axis equal;
colormap('gray');
imagesc(im3);
hold on;
cols1 = size(im1,2);
loc1 = xFI(1:2,:)';
loc2 = xFJ(1:2,:)';
plot(loc1(:,1), loc1(:,2), 'ro', 'LineWidth', 2);
plot(loc2(:,1)+cols1, loc2(:,2), 'yo', 'LineWidth', 2);
for i = 1: size(loc1,1)
    if isnan(loc2(i,1)) == 0
        line([loc1(i,1) loc2(i,1)+cols1], [loc1(i,2) loc2(i,2)], 'Color', 'c', 'LineWidth', 2);
    end
end

rectangle('Position', [BB1(1) BB1(2) BB1(3)-BB1(1) BB1(4)-BB1(2)], 'EdgeColor', 'g', 'LineWidth', 4);
rectangle('Position', [BB2(1)+cols1 BB2(2) BB2(3)-BB2(1) BB2(4)-BB2(2)], 'EdgeColor', 'y', 'LineWidth', 4);
axis off;
hold off;

% Return a new image that appends the two images side-by-side.

function im = appendimages(image1, image2)

% Select the image with the fewest rows and fill in enough empty rows
%   to make it the same height as the other image.
rows1 = size(image1,1);
rows2 = size(image2,1);

if (rows1 < rows2)
     image1(rows2,1) = 0;
else
     image2(rows1,1) = 0;
end

% Now append both images side-by-side.
im = [image1 image2];

author: yuxng


text:
function varargout = patchGraph(nodes, edges, faces) %#ok<INUSL>
%PATCHGRAPH Transform 3D graph (mesh) into a patch handle
%
%   [PX, PY, PZ] = PATCHGRAPH(NODES, EDGES, FACES)
%   Transform the graph defined as a set of nodes, edges and faces in a
%   patch which can be drawn usind matlab function 'patch'.
%   The result is a set of 3 array of size [NV*NF], with NV being the
%   number of vertices per face, and NF being the total number of faces.
%   each array contains one coordinate of vertices of each patch.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 28/06/2004.
%

if iscell(faces)
    p = zeros(length(faces), 1);
    for i = 1:length(faces)
        p(i) = patch( ...
            'Faces', faces{i}, ...
            'Vertices', nodes, ...
            'FaceColor', 'r', ...
            'EdgeColor', 'none') ;
    end
else
    p = patch( ...
        'Faces', faces, ...
        'Vertices', nodes, ...
        'FaceColor', 'r', ...
        'EdgeColor', 'none') ;
end

if nargout>0
    varargout{1}=p;
end


text:
function saveash5(data, path_to_file, options)

% save the variable stored in data in a h5 file named path_to_file (includes the
% path to the directory). The file is saved with the same class as data.
% options.groupname: the name of the group (default: /mov)
% options.append: append to existing file if it exists (default: True)

if ~exist('path_to_file','var'); path_to_file = 'data.h5'; end
[pathstr, name, ~] = fileparts(path_to_file);
path_to_file = fullfile(pathstr,[name,'.h5']); % make sure file has an h5 extension
data_type = class(data);

sizY = size(data);
nd = ndims(data)-1;

defoptions.groupname = '/mov';
defoptions.append = true;

if ~exist('options','var'); options = defoptions; end
if ~isfield(options,'groupname'); options.groupname = defoptions.groupname; end
if ~isfield(options,'groupname'); options.append = defoptions.append; end

if ~exist(path_to_file,'file') || ~options.append
    if ~options.append && exist('path_to_file','file')
        warning('Rewriting existing file')
    end
    h5create(path_to_file,options.groupname,[sizY(1:nd),Inf],'Chunksize',[sizY(1:nd),100],'Datatype',data_type);
    start_point = ones(1,nd+1);
else
    fprintf('Appending to existing file \n')
    info = h5info(path_to_file);
    options.groupname = ['/',info.Datasets.Name];
    start_point = [ones(1,nd),info.Datasets.Dataspace.Size(end)+1];
end

h5write(path_to_file,options.groupname,data,start_point,size(data));

author: flatironinstitute

repo: NoRMCorre


text:
function lll = dnetLowerBound(model)

% DNETLOWERBOUND Computes lower bound on log likelihood for an DNET model.
% FORMAT
% DESC computes the variational lower bound on the log likelihood
% of a mixtures of probabilistic PCA model.
% ARG model : the model for which log likelihood is to be computed.
% RETURN lll : the lower bound on the log likelihood computed for the model.
%
% SEEALSO : dnetCreate, modelLogLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

lll = 0.5*model.N*model.d*log(model.beta/(2*pi)) ...
      - model.N*log(model.M);

lll = lll - 0.5*model.alpha*sum(sum(model.A.*model.A)) ...
      - 0.5*model.alpha*sum(model.b.*model.b);

lll = lll - sum(sum(xlogy(model.w)));

% Get projections of latent samples.
Ypred = dnetOut(model);
if model.N > model.M
  for i = 1:model.M
    diffY = model.y - repmat(Ypred(i, :), model.N, 1);
    diffY = diffY.*diffY.*repmat(model.w(:,i), 1, model.d);
    lll = lll - 0.5*model.beta*sum(sum(diffY));
  end
else
  for i = 1:model.N
    diffY = repmat(model.y(i, :), model.M, 1) - Ypred;
    diffY =diffY.*diffY.*repmat(model.w(i,:)', 1, model.d);
    lll = lll - 0.5*model.beta*sum(sum(diffY));
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ all_bboxes_gt ] = get_grount_truth_bboxes_from_voc( voc_path, image_set, voc_year, with_hard_samples, cache_dir )

if ~exist('cache_dir','var')
    cache_dir = '.';
end

VOCopts            = initVOCOpts( voc_path, voc_year );
fullimage_set_name = ['voc_', voc_year, '_' image_set];
addpath([voc_path, filesep, 'VOCcode']);
cache_file = [cache_dir, filesep, 'gt_bbox_', fullimage_set_name];
if with_hard_samples, cache_file = [cache_file, '_with_hard_samples']; end
cache_file = [cache_file, '.mat'];

if exist(cache_file,'file')
    all_bboxes_gt = loadGroundTruthBBoxes(cache_file);
else
    class_to_id     = containers.Map(VOCopts.classes, 1:length(VOCopts.classes));
    image_ids       = textread(sprintf(VOCopts.imgsetpath, image_set), '%s');
    num_imgs        = length(image_ids);
    all_bboxes_gt   = cell(num_imgs, 1);
    for img_idx = 1:num_imgs
        if mod(img_idx,500) == 0, fprintf('Load gt bboxes::%s %d/%d\n', fullimage_set_name, img_idx, num_imgs); end
        all_bboxes_gt{img_idx} = getGroundTruthBBoxes(sprintf(VOCopts.annopath, image_ids{img_idx}), class_to_id, with_hard_samples);
    end
    saveGroundTruthBBoxes(cache_file, all_bboxes_gt);
end

end

function gt_bboxes = getGroundTruthBBoxes(filename, class_to_id, with_hard_samples)
try
    voc_rec       = PASreadrecord(filename);
    valid_objects = 1:length(voc_rec.objects(:));
    is_difficult  = cat(1, voc_rec.objects(:).difficult);
   if ~with_hard_samples
       valid_objects = valid_objects(~is_difficult);
       is_difficult  = is_difficult(~is_difficult);
   end

   gt_bboxes = single(cat(1, voc_rec.objects(valid_objects).bbox));
   gt_class_idx = class_to_id.values({voc_rec.objects(valid_objects).class});
   gt_class_idx = single(cat(1, gt_class_idx{:}));
   gt_bboxes    = single([gt_bboxes, gt_class_idx, is_difficult(:)]);
catch
    gt_bboxes = zeros(0, 6, 'single');
end

end

function all_bboxes_gt = loadGroundTruthBBoxes(filename)
load(filename, 'all_bboxes_gt');
end

function saveGroundTruthBBoxes(filename, all_bboxes_gt)
save(filename, 'all_bboxes_gt', '-v7.3');
end

author: gidariss

repo: mrcnn-object-detection


text:
  % 1 x 1  % 1 x 1% plot_info
% plot Free energy and other infomation
a_min = 1e-5;
NX = 2;
NY = 2;
alpha_plot = 0;

FE  = Info.FE;

tend = length(FE);

% Start and End of iteration number for plot
tmode = 2;
t0 = 10;
t3 = 600;
switch	tmode
case	1
	t1 = t0;
	t2 = tend;
case	2
	t1 = t3;
	t2 = tend;
case	3
	t1 = t0;
	t2 = t3;
end

Hmode = 1;
Fmode = 2;
% Fmode = 1 :Remove sudden drop of free energy due to weight pruning

if Fmode == 1
	% Remove sudden drop of free energy due to weight pruning
	Fdif = diff(FE);

	ixz  = find( Fdif < 0);
	ixz  = [ixz, ixz + 1, ixz + 2];
	ixz  = ixz(:);

	ix = setdiff2(1:tend,ixz);
else
	ix = 1:tend;
end

% Free energy
nfig = 1;
subplot(NX,NY,nfig)
%hold on
plot(ix,FE(ix))

title('Free energy')
xlim([t1 t2])

%return

% Model entropy

if	Hmode < 2 &isfield(Info,'H')
	nfig = nfig + 1;
	subplot(NX,NY,nfig)

	if isfield(parm,'T')
		switch	Hmode
		case 0
			H   = -Info.H * parm.T/log(parm.T);
		case 1
			H   = -Info.H /log(parm.T);
		end
	else
		H   = -Info.H ;
	end

	plot(ix,H(ix))
	hold on
	title('Effective parameter number')
	xlim([t1 t2])
end

if	Hmode == 2 & isfield(Info,'M')
	nfig = nfig + 1;
	subplot(NX,NY,nfig)

	plot(ix,Info.M(ix), '-b')
	title('Number of input dimension')
	xlim([t1 t2])
end

if	isfield(Info,'LP')
	nfig = nfig + 1;
	subplot(NX,NY,nfig)
	LP  = Info.LP;
	plot(ix,LP(ix))
	title('Expected log-likelihood')
	xlim([t1 t2])
end

if	isfield(Info,'Err')
	nfig = nfig + 1;
	subplot(NX,NY,nfig)
	Err  = Info.Err;
	plot(ix,Err(ix))
	title('Error')
	xlim([t1 t2])
end

if alpha_plot == 0, return; end;

% Hyper variance parameter
if isfield(Model,'A')
	if isfield(Model,'ix_act')
		N = size(Model.A ,1);
		A = zeros(N,Model.M_all);
		A(:,Model.ix_act) = Model.A;
	else
		A =  Model.A;
	end

    nfig = nfig + 1;
	subplot(NX,NY,nfig)
	plot(A')
	title('Hyper variance parameter')

    nfig = nfig + 1;
	subplot(NX,NY,nfig)
	hist(Model.A, 100)
	title('Hyper variance parameter')

	M=length(Model.A);
	K=sum(Model.A < a_min);
	fprintf('# of active dim = %d, # of near zero = %d\n', M, K)
end

if alpha_plot == 1, return; end;

if isfield(Info, 'A')
    AA  = Info.A ;
    nfig = nfig + 1;
	subplot(NX,NY,nfig)
	plot(AA')
	title('Hyper variance histry')
end

author: KamitaniLab

repo: GenericObjectDecoding


text:
function [tt]=round(tt,varargin)
%Approximate QTT-Tucker with another one with specified accuracy
%   [QTT]=ROUND(QTT,EPS) Approximate QTT-Tucker tensor with relative
%   accuracy EPS
%
%   [QTT]=ROUND(QTT,EPS,RMAX) Approximate QTT-Tucker tensor with relative
%   accuracy
%   EPS and maximal rank RMAX. RMAX can be array of ranks or a number
%
% Please see @qtt_tucker/round2 for a more accurate version of the algorithm
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
d=tt.dphys;
core=tt.core;
tuck=tt.tuck;
eps=varargin{1};
tolcorr = 0;
for i=1:d
    tolcorr = tolcorr+tuck{i}.d;
end;
rmax = [];
if (nargin==3)
    rmax = varargin{2};
end;
ismatrix = 0;
if (isa(tuck{1}, 'tt_matrix'))
    ismatrix = 1;
    curn = cell(d,1);
    curm = cell(d,1);
end;
for i=1:d
    if (ismatrix)
        curn{i} = tuck{i}.n;
        curm{i} = tuck{i}.m;
        tuck{i} = tt_tensor(tuck{i});
    end;
   [tuck{i},rm]=qr(tuck{i},'lr');
   core{i}=ten_conv(core{i},2,rm.');
end
if (isempty(rmax))
    core=round(core,eps*sqrt(d)/sqrt(tolcorr));
else
    core=round(core,eps*sqrt(d)/sqrt(tolcorr),rmax);
end;
%Round the core --- we know the result comes
%with rl orthogonality? -< No, we don't
[core, nrm] = qr(core, 'lr');
core{d} = core{d}*nrm;
rtt=rank(core);
n=size(core);
for i=d:-1:1
   cr=reshape(core{i},[rtt(i),n(i),rtt(i+1)]);
   cr=permute(cr,[2,1,3]); cr=reshape(cr,n(i),rtt(i)*rtt(i+1));
   [u,s,v]=svd(cr,'econ');
   s=diag(s);
   r=my_chop2(s,norm(s)*eps/sqrt(tolcorr));
   if (~isempty(rmax))
       r = min(r,rmax);
   end;
   u=u(:,1:r); s=s(1:r); v=v(:,1:r);
   tuck{i}=tuck{i}*(u*diag(s));
   if (isempty(rmax))
       tuck{i}=round(tuck{i},eps*sqrt(tuck{i}.d)/sqrt(tolcorr));
   else
       tuck{i}=round(tuck{i},eps*sqrt(tuck{i}.d)/sqrt(tolcorr), rmax);
   end;
   [tuck{i},rm]=qr(tuck{i},'lr');
   cr=rm*v';
   cr=reshape(cr,[r,rtt(i),rtt(i+1)]);
   % Shift QR to the next core block
   if (i>1)
       cr=permute(cr,[1,3,2]);
       cr = reshape(cr, r*rtt(i+1), rtt(i));
       [cr, rv] = qr(cr, 0);
       cr2 = core{i-1};
       rtuck2 = size(cr2, 2);
       cr2 = reshape(cr2, rtt(i-1)*rtuck2, rtt(i));
       cr2 = cr2*(rv.');
       rtt(i) = size(cr, 2);
       core{i-1} = reshape(cr2, rtt(i-1), rtuck2, rtt(i));
       core{i} = reshape(cr.', rtt(i), r, rtt(i+1));
   else
       core{i}=permute(cr,[2,1,3]);
   end;
end
if (ismatrix)
    for i=1:d
        tuck{i} = tt_matrix(tuck{i}, curn{i}, curm{i});
    end;
end;
tt.core=core;
tt.tuck=tuck;
return
end

author: oseledets

repo: TT-Toolbox


text:
function M = dunzip(Z)
% DUNZIP - decompress gzipped stream of bytes
% FORMAT M = dzip(Z)
% Z  -  compressed variable to decompress (uint8 vector)
% M  -  decompressed output
%
% See also DZIP

% Carefully tested, but no warranty; use at your own risk.
% Michael Kleder, Nov 2005
% Modified by Guillaume Flandin, May 2008

import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
a   = java.io.ByteArrayInputStream(Z);
b   = java.util.zip.InflaterInputStream(a);
isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
c   = java.io.ByteArrayOutputStream;
isc.copyStream(b,c);
M   = c.toByteArray;

author: brainstorm-tools


text:
function point = polygonPoint(poly, pos)
%POLYGONPOINT Extract a point from a polygon
%
%   POINT = polygonPoint(POLYGON, POS)
%
%
%   Example
%   polygonPoint
%
%   See also
%   polygons2d, polylinePoint
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2009-04-30,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

% eventually copy first point at the end to ensure closed polygon
if sum(poly(end, :) == poly(1,:))~=2
    poly = [poly; poly(1,:)];
end

% number of points to compute
nPoints = length(pos(:));

% number of vertices in polygon
Nv = size(poly, 1)-1;

% allocate memory results
point = zeros(nPoints, 2);

% iterate on points
for i = 1:nPoints
    % compute index of edge (between 0 and Nv)
    ind = floor(pos(i));

    % special case of last point of polyline
    if ind==Nv
        point(i,:) = poly(end,:);
        continue;
    end

    % format index to ensure being on polygon
    ind = min(max(ind, 0), Nv-1);

    % position on current edge
    t = min(max(pos(i)-ind, 0), 1);

    % parameters of current edge
    x0 = poly(ind+1, 1);
    y0 = poly(ind+1, 2);
    dx = poly(ind+2,1)-x0;
    dy = poly(ind+2,2)-y0;

    % compute position of current point
    point(i, :) = [x0+t*dx, y0+t*dy];
end


text:
function z = svec( x, nrm )

if nargin < 2 || isempty( nrm ) || isequal( nrm, 'fro' ),
    nrm = 2;
elseif ~isnumeric( nrm ) || length( nrm ) ~= 1 || nrm < 1,
    error( 'Second argument must be a number between 1 and Inf, or ''fro''.' );
end

if ~isreal( x ) && nrm ~= 2,
    z = vec( x );
    return
else
    [ xR, y ] = bcompress( x );
    if isempty( y ),
        z = cvx( 0 );
    else
        z = y .* norms( xR, nrm, 2 );
    end
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [w,run] = train_mis(x,w,lambda)
% Modified iterative scaling
% x is premultiplied by y

% Written by Thomas P Minka

if nargin < 3
  lambda = 0;
end
if lambda > 0
  error('must have lambda = 0')
end
[d,n] = size(x);
flops(0);
step = 1/max(sum(abs(x),1));

i1 = (x > 0);
x1 = abs(x).*i1;
x2 = abs(x).*(1-i1);
flops(flops + 3*d*n);
if nargout > 1
  run.w = [];
  run.flops = [];
  run.e = [];
end
for iter = 1:10000
  old_w = w;
  % s1 = 1-sigma
  s1 = 1./(1+exp(w'*x));
  delta = (x1*s1')./(x2*s1');
  w = w + step*0.5*log(delta);
  if iter == 1
    % same for every iteration
    % use spmul because x1,x2 have structural zeros
    fl = flops_mul(w',x)+n*(flops_exp+2) + ...
	flops_spmul(x1,s1')+flops_spmul(x2,s1')+d + ...
	d*(flops_exp+2);
  end
  flops(flops + fl);

  if nargout > 1 & rem(iter,100) == 1
    run.w(:,end+1) = w;
    run.flops(end+1) = flops;
    run.e(end+1) = logProb(x,w) -0.5*lambda*w'*w;
  end
  if rem(iter,1000) == 0
    fprintf('MIS iter %d\n', iter)
  end

  if max(abs(w - old_w)) < 1e-6
    break
  end
end
if iter == 10000
  warning('not enough iters')
end
if nargout > 1
  figure(2)
  plot(run.e)
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function imStack=analyze2mrLoadRetTSeries(inFileRoot,outFileRoot,nVols,firstVolIndex,doRotate,scaleFact,flipFlag)
% function imStack=analyze2mrLoadRetTSeries(inFileRoot,outFileRoot,nVols,firstVolIndex,doRotate,scaleFact,flipFlag)
% Converts from analyze functional image data to mrLoadRet TSeries format
% Analyze functional data are stored as N individual volumes, each with S slices
% mrLoadRet has S individual files with N acquisitions in each.
% the doRotate param allows to you roate the functional data by doRotate*90 degrees
% Try to read the the first volume to see if it's there and get the dimensions of all the rest

% TODO: 1:  Rotate these so that they are in the same orientation as the anatomicals
% (Maybe make this a flag)
% 2: Why can't we see the anatomicals? Claims not to be able to read from /Raw/Anatomy/Inplane
% 3: Why can't we load in the corAnal? We do this by hand perfectly well but it won't do it from mrLoadRet
% 4: Do we want to resample the functional data? Or is this the moment to go to mrLoadRet3?
%

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

if (~exist('scaleFact','var'))
    scaleFact=[1 1]; % No interpolation. (Scaling=1)
end

if (length(scaleFact)~=2) % If we just get a scalar for the scale factor, assume that it applies in both dimensions
    scaleFact=repmat(scaleFact(1),2);
end

if (~exist('doRotate','var'))
    doRotate=1; % This is on by default. Rotates 1*90 degrees
end

if (~exist('flipFlag','var'))
    flipFlag=1; % This is on by default. Flips up/down after rotation
end
suffix=sprintf('%03d',firstVolIndex);

fileName=[inFileRoot,suffix,'.hdr'];

V=spm_vol(fileName);
im=spm_read_vols(V);

[y,x,nSlices]=size(im);
fprintf('Read in a volume of size %d, %d, %d',y,x,nSlices);

% Pre-allocate memory: This is a monster and will fail on many machines.
fprintf('\nTrying to allocate an array with %d elements...\n',y*x*nSlices*nVols);

if (mod(doRotate,2)) % When we rotate by 180 degrees the x and y dimensions remain unchanged
    funcVol=zeros(y,x,nSlices,nVols);
else
    funcVol=zeros(x,y,nSlices,nVols);
end
 fprintf('Rotating by %d x 90, flipFlag=%d',doRotate,flipFlag);

for t=0:(nVols-1)
    thisImIndex=t+firstVolIndex;
    suffix=sprintf('%03d',thisImIndex);
    fileName=[inFileRoot,suffix];
    V=spm_vol(fileName);
    im=spm_read_vols(V);

    % Do the rotation and scaling

    if (mod(doRotate,2)) % When we rotate by 180 degrees the x and y dimensions remain unchanged
        im2=zeros(y,x,nSlices);
    else
        im2=zeros(x,y,nSlices);
    end
    fprintf('\nVol=%d',thisImIndex);

    for thisSlice=1:nSlices
        imSlice=squeeze(im(:,:,thisSlice));
    %s    imSlice=imresize(imSlice,[scaleFact(1)*y,scaleFact(2)*x],'nearest');

        im2(:,:,thisSlice)=rot90(imSlice,doRotate);
        if (flipFlag)

            im2(:,:,thisSlice)=flipud(im2(:,:,thisSlice));
        end

    end % next imSlice

    funcVol(:,:,:,t+1)=im2;
    %fprintf('.');

end
size(funcVol);
[y, x, nSlices, nVols]=size(funcVol);

% Now write them out in a different format
fprintf('\nDone reading data: Writing now...\n');
for t=1:nSlices
    suffix=int2str(t);
    fileName=[outFileRoot,suffix,'.dat'];
    tSeries=squeeze(funcVol(:,:,t,:));
    [a,b,c]=size(tSeries);
    fprintf('\nSize before: %d by %d by %d vols ',a,b,c);

    tSeries=squeeze(shiftdim(tSeries,2));
    [a,b,c]=size(tSeries);
    fprintf('\nSize after:%d slices by %d by %d',a,b*scaleFact(1),c*scaleFact(2));

    savetSeriesDat3(fileName,tSeries,y,x,scaleFact);

    fprintf('_');
end
fprintf('\nDone\n');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function center = findCenter( PSF )
%
%       center = findCenter( PSF );
%
%  Given (possibly) several PSF images, this function sets the center
%  of the PSFs (location of point source) to be the location of max
%  entry of the PSFs.
%
%  Input:
%         PSF  -  cell array containing the PSF image(s)
%
%  Output:
%         center - array [row_index, col_index] containing the
%                  location of the center of the PSF.
%

%  J.Nagy 1/8/02

center = cell(size(PSF));

%
%  Note that this should work for 2D as well as 3D images.  We simply
%  find the maximum entry in the center of the image.
%
for k = 1:size(PSF, 3)
  for i = 1:size(PSF, 1)
    for j = 1:size(PSF, 2)
       P = PSF{i,j,k};
       [m, n, l] = size(P);
       idx_row = floor((m+1)/2):ceil((m+1)/2);
       idx_col = floor((n+1)/2):ceil((n+1)/2);
       idx_depth = floor((l+1)/2):ceil((l+1)/2);
       P2 = zeros(size(P));
       P2(idx_row, idx_col, idx_depth) = P(idx_row, idx_col, idx_depth);
       [ci, cj, ck] = ind2sub(size(P2), find( P2 == max( P2(:) ) ));
       c = [min(ci), min(cj), min(ck)];
       center{i,j,k} = c(1:length(size(P2)));
    end
  end
end

repo: IRtools


text:
function xform = mrAlignMI(sessionDir, vAnatFile, inplaneFile, initXFM, nuCorrect, forceSave)
% function xform = mrAlignMI(sessionDir, vAnatFile, inplaneFile, initXFM, nuCorrect, forceSave)
%
% * sessionDir is where the relevent mrSESSION lives. Defaults to pwd.
%
% * vAnatFile: default is a file browser if it's not found in the usual places.
%
% * inplaneFile: we currently need the first and last slice from the original
% inplane data to compute the coarse alignment. In the future, this info
% should be stored in the Inplane/anat.mat file.
%
% * initXFM : define [x y z] translations (first 3 params of the
% terminal output) from the automatically  computed coarse
% alignment. This is usuful when the algorithm  converges to an
% obviously wrong location (local minima) by allowing a different
% starting position.
%
% forceSave: flag to force saving computed xofrm to the mrSESSION file.
% Default is to show the computed alignment and then ask the user if they
% want to save.
%
% Tries to compute an inplane-to-vAnatomy alignment fully automatically.
% The core algorithm is quite simple (once you get past the file management
% stuff). It generates a coarse alignment by:
%  1. assuming that the vAnatomy is in the standard vAnatomy orientation
%     (sagittal, left-is-left, yxz, etc), and
%  2. using the original inplane header info to bring the inplane into a
%     standard orientation.
% It then uses a rigid body mutual information registration algorithm (from
% spm2) to refine the alignmnet.
%
% It seems to work perfectly when you have non-surface coil, whole-brain or
% most-brain data. YMMV with other data varieties.
%
% HISTORY:
% 2005.05.27 RFD (bob@sirl.stanford.edu) wrote it.
% 2005.06.27 RFD: slight changes to allow gzipped I-files.
% 2005.10.60 SOD: modified so it can refine existing mrSESSION.alignment's,
% put in a correction for an intensity gradient which seems
% to improve the alignment for surface-coil data, and put in
% another check that i foud useful. They come unfortenately with
% two more dialogs.
% 2006.04.03 SOD: added initXFM option

if(isempty(which('spm_coreg')))
    error('Requires spm2 tools! Try adding .../matlab/toolbox/mri/spm2 to your path.');
end
if(~exist('sessionDir','var') || isempty(sessionDir))
    sessionDir = pwd;
end
if(~exist('forceSave','var') || isempty(forceSave))
    forceSave = false;
end
if (~exist('nuCorrect','var') || isempty(nuCorrect))
    nuCorrect = [];
end
loadSession(sessionDir);
global HOMEDIR;
HOMEDIR = sessionDir;
global mrSESSION;
if(~exist('vAnatFile','var') || isempty(vAnatFile))
    % Try to find it
    vAnatFile = getVAnatomyPath;
end
if(~exist('inplaneFile','var') || isempty(inplaneFile))
    % Try to find it
    inplaneDir = fullfile(sessionDir,'Raw','Anatomy','Inplane');
    d = dir(fullfile(inplaneDir, 'I*001*'));
    if ~isempty(d), % ie any Ifiles are found
        inplaneFile = fullfile(inplaneDir, d(1).name);
        if(~exist(inplaneFile,'file'))
            [f,p] = uigetfile({'*.dcm','Dicom';'*.001','I-file';'*.gz','gzipped';'*.*','all'}, ...
                'Select a slice from raw Inplane files...', inplaneFile);
            if(isequal(f,0) || isequal(p,0))
                disp('User canceled.'); return;
            else
                inplaneFile = fullfile(p,f);
            end
        end;
    end
end
inplaneAnatFile = fullfile(sessionDir,'Inplane','anat.mat');
if(~exist(inplaneAnatFile,'file'))
    [f,p] = uigetfile({'*.mat','mat';'*.*','all'}, ...
        'Select the Inplane anat.mat file...', inplaneAnatFile);
    if(isequal(f,0) || isequal(p,0))
        disp('User canceled.'); return;
    else
        inplaneAnatFile = fullfile(p,f);
    end
end
disp(inplaneAnatFile);
disp(vAnatFile);

% check if alignment exists and give option to use this as a
% starting point allowing refinements of an existing alignment
% only do this if no initXFM is defined
useStartingEstimate = 0;
if notDefined('initXFM'),
    if isfield(mrSESSION,'alignment'),
        bn = questdlg(['mrSESSION.alignment exists. Do you want to to use ' ...
            'it as a starting estimate?'],...
            'starting estimate','Yes','No','No');
        if strmatch(bn,'Yes'),
            useStartingEstimate = 1;
        end;
    end;
end;

[v.img,v.mm] = readVolAnat(vAnatFile);
sz = size(v.img);
% Build a very rough xform to axial, ac-pc space
hsz = sz./2;
% ***FIXME: I'm not sure the v.mm is correct here. I've only tested it on
% isotropic data so far, so I can't tell.
vAnatAcpcXform = [0 0 v.mm(3) -hsz(3); 0 -v.mm(2) 0 hsz(1); -v.mm(1) 0 0 hsz(2); 0 0 0 1];
VF.uint8 = uint8(v.img);
VF.mat = vAnatAcpcXform;

ip = load(inplaneAnatFile);
% Correct for intensity gradient.
% yet another anoying dialog, sorry...
% For surface-coil data the results look better with correction - at
% least for the subject I tested.
if isempty(nuCorrect),
    bn = questdlg('Do you want to correct for an intensity gradient?',...
        'intensity gradient','Yes','No','No');
    if strmatch(bn,'Yes'),
        nuCorrect = true;
    else
        nuCorrect = false;
    end
end
if nuCorrect,
    [GradInt GradNoise]  = regEstFilIntGrad(ip.anat);
    ip.anat              = regCorrIntGradWiener(ip.anat,GradInt,GradNoise);
end;

% **** FIXE ME: we should compute and store the cannonical xform when we
% create the anat.mat file. Also, we can do a bit better for the initial
% alignment by taking the crop into account in the translations.
%sz = size(ip.anat);
% We assume that the brain is centered in the crop region.
%cropOffset = ip.inplanes.crop(1,:)-(ip.inplanes.fullSize-diff(ip.inplanes.crop))/2;
%cropOffset = ip.inplanes.crop(1,[2,1]); % OK
% apparently ip.inplanes.crop will not exist for older versions
if ~isfield(ip,'inplanes'),
    sz               = size(ip.anat);
    ip.inplanes.crop = [1 1;sz(1) sz(2)];
end
% not used?!
%cropOffset = [ip.inplanes.crop(1,[2,1])./2 0];%4 (best?)
%cropOffset = [eye(3) cropOffset'; 0 0 0 1];

% If Ifiles exist use them for an initial xform otherwise center on 0.
if ~ieNotDefined('inplaneFile'),
    xformToScanner = inv(computeXformFromIfile(inplaneFile,ip.inplanes.crop));
    xformToScanner(1:3,4) = xformToScanner(1:3,4)+[10 -20 -20]';
    % GE convention is right-to-left, but we want l-to-r (Talairach convention)
    %lrFlip = [-1 0 0 ip.inplanes.fullSize(1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
    %lrFlip = [-1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
    VG.mat = xformToScanner;
else
    % default xform (center on 0,0,0)
    disp(sprintf('[%s]:No Ifiles found centering volume on 0.',mfilename));
    voxSize = mrSESSION.inplanes.voxelSize;
    VG.mat = [diag(voxSize),-(size(ip.anat)'.*voxSize'./2); 0 0 0 1];
    % GE convention is right-to-left, but we want l-to-r (Talairach
    % convention)
    % Careful! If you rerun you don't want to keep flipping your image
    % all the time
    VG.mat(1,:) =  -VG.mat(1,:);
    %VG.mat(2,:) =  -VG.mat(2,:);
    %VG.mat(3,:) =  -VG.mat(3,:);
    VG.mat
end;
%ip.anat = mrAnatHistogramClip(ip.anat, 0.4, 0.99);
ip.anat = mrAnatHistogramClip(ip.anat, 0.2, 0.99);
VG.uint8 = uint8(ip.anat*255+0.5);

if useStartingEstimate,
    revAlignment = spm_imatrix(VF.mat*mrSESSION.alignment/VG.mat);
    flags.params = revAlignment(1:6);
    disp('flags.params :');disp(flags.params);
    flags.sep    = [8 4 2 1];
    rotTrans     = spm_coreg(VG,VF,flags);
    oldAlignment = mrSESSION.alignment;
    disp('Old alignment :');disp(oldAlignment);
else
    % allow manually reset of the starting position
    if ~exist('initXFM','var')
        initXFM =[];
    end
    % max 6 params if less fill with zeros
    if numel(initXFM) < 6,
        flags.params = [initXFM(:)' zeros(1,6-numel(initXFM))];
    else
        flags.params = initXFM;
    end;
    % for dire need
    flags.sep    = [16 8 4 2];
    rotTrans     = spm_coreg(VG,VF,flags);
end;
%xform = spm_matrix(rotTrans)*VG.mat\VF.mat;
% volCoords = xform*[inplaneCoords]
alignment = VF.mat\spm_matrix(rotTrans)*VG.mat;
disp('Alignment matrix :');disp(alignment);

% output
xform = alignment;

sz = size(v.img);
bb = [1,1,1; sz];
ipVol = mrAnatResliceSpm(ip.anat, inv(alignment), bb, [1 1 1], [1 1 1 0 0 0], 0);
ipVol(isnan(ipVol)) = 0;
if(forceSave)
    % Save an image so user can check the alignment
    volIm = makeMontage(v.img); volIm = volIm./255;
    ipIm = makeMontage(ipVol);
    ipIm(ipIm<0)=0;ipIm(ipIm>1)=1;
    ipIm = cat(3,ipIm,volIm,volIm);
    imwrite(ipIm, fullfile(sessionDir,'align.png'));
    % just save it to existing mrSESSION
    mrSESSION.alignment = alignment;
    saveSession;
    % Delete old files that were built using previous alignment.
    cleanAllFlats; cleanGray; cleanVolume;
else
    disp('checking alignment...');
    % another check that i (SOD) found usefull
    % coronal (rotated...)
    % overlayVolumes(shiftdim(ipVol,2),shiftdim(v.img,2));
    % sagittal
    overlayVolumes(ipVol,v.img);
    % horizontal
    overlayVolumes(shiftdim(ipVol,1),shiftdim(v.img,1));

    volIm = makeMontage(v.img,40:9:sz(3)-40); volIm = volIm./255;
    if ~useStartingEstimate,
        ipVolInitial = mrAnatResliceSpm(ip.anat, inv(VF.mat\eye(4)*VG.mat), ...
            bb, [1 1 1], [1 1 1 0 0 0], 0);
        titlestring  = 'Initial coarse alignment';
    else
        ipVolInitial = mrAnatResliceSpm(ip.anat, inv(oldAlignment), ...
            bb, [1 1 1], [1 1 1 0 0 0], 0);
        titlestring  = 'Initial alignment from mrSESSION.alignment';
    end;
    ipVolInitial(isnan(ipVolInitial)) = 0;
    %figure;imagesc(makeMontage(ipVolInitial));axis equal;
    ipIm = makeMontage(ipVolInitial,40:9:sz(3)-40);
    ipIm(ipIm<0)=0;ipIm(ipIm>1)=1;
    figure;
    subplot(2,1,1);
    image(cat(3,ipIm,volIm,volIm)); axis equal tight off;
    title(titlestring);
    subplot(2,1,2);
    ipIm = makeMontage(ipVol,40:9:sz(3)-40);
    ipIm(ipIm<0)=0;ipIm(ipIm>1)=1;
    image(cat(3,ipIm,volIm,volIm)); axis equal tight off;
    title('Final alignment');

    [f,p] = uigetfile('*.mat','Save alignment to mrSESSION file?',fullfile(sessionDir,'mrSESSION.mat'));
    if(isequal(f,0)||isequal(p,0))
        error('Cancelled- alignment NOT SAVED.');
    else
        mrSESSION.alignment = alignment;
        saveSession;
        % Delete old files that were built using previous alignment.
        cleanAllFlats; cleanGray; cleanVolume;
    end
end
return;


text:
% The layout of the configuration tree and the types of configuration items
% have been kept compatible to a configuration system and job manager
% implementation in SPM5 (Statistical Parametric Mapping, Copyright (C)
% 2005 Wellcome Department of Imaging Neuroscience). This code has been
% completely rewritten based on an object oriented model of the
% configuration tree.

author: spm


text:
function [n,V,p] = affine_fit(X,W)
    %Computes the plane that fits best (lest square of the normal distance
    %to the plane) a set of sample points.
    %INPUTS:
    %
    %X: a N by 3 matrix where each line is a sample point
    %
    %OUTPUTS:
    %
    %n : a unit (column) vector normal to the plane
    %V : a 3 by 2 matrix. The columns of V form an orthonormal basis of the
    %plane
    %p : a point belonging to the plane
    %
    %NB: this code actually works in any dimension (2,3,4,...)
    %Author: Adrien Leygue
    %Date: August 30 2013

    if nargin<2
      W = 1/size(X,1)*ones(size(X,1),1);
    end
    W = W(:)./sum(W);

    %the mean of the samples belongs to the plane
    p = sum(bsxfun(@times,W,X));

    %The samples are reduced:
    R = bsxfun(@times,bsxfun(@minus,X,p),W);
    %Computation of the principal directions if the samples cloud
    [V,D] = eig(R'*R);
    %Extract the output from the eigenvectors
    n = V(:,1);
    V = V(:,2:end);
end

% Copyright (c) 2013, Adrien Leygue
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

author: alecjacobson

repo: gptoolbox


text:
function [mi] = calculate_mutual_information_array(data)
% FUNCTION [MI_ARRAY] = CALCULATE_MUTUAL_INFORMATION_ARRAY(DATA)
% calculates the mutual information between all pairs of variables
% Data must be discrete, and take values 1,2,...,size
% data(i,m) is the node i in the case m.

[num_nodes num_examples] = size(data);

node_sizes = max(data');
for i = 1:num_nodes
  for ic = 1:node_sizes(i) % I CLASS ic
    px(i,ic) = sum(data(i,:)==ic);
    for j = 1:num_nodes    % J CLASS jc
      for jc = 1:node_sizes(j)
        pxy(i,ic,j,jc) = sum( (data(i,:)==ic) & (data(j,:)==jc) );
      end
      mi(i,j) = 0;
    end
  end
end

for i = 1:num_nodes
  for ic = 1:node_sizes(i)
    for j = 1:num_nodes
      for jc = 1:node_sizes(j)
        if( pxy(i,ic,j,jc)~=0 & px(i,ic)~=0 & px(j,jc)~= 0)
          mi(i,j) = mi(i,j) + pxy(i,ic,j,jc)*log2( num_examples*pxy(i,ic,j,jc)/(px(i,ic)*px(j,jc)) )/num_examples;
        end
      end
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [cfg] = ft_interactiverealign(cfg)

% FT_INTERACTIVEREALIGN allows the user to interactively translate, rotate and scale an
% individual geometrical object to a template geometrical object. It can for example be used
% to align EEG electrodes to a model of the scalp surface.
%
% Use as
%   [cfg] = ft_interactiverealign(cfg)
%
% The configuration structure should contain the individuals geometrical object that
% has to be realigned
%   cfg.individual.elec           = structure, see FT_READ_SENS
%   cfg.individual.grad           = structure, see FT_READ_SENS
%   cfg.individual.opto           = structure, see FT_READ_SENS
%   cfg.individual.headmodel      = structure, see FT_PREPARE_HEADMODEL
%   cfg.individual.headshape      = structure, see FT_READ_HEADSHAPE
%   cfg.individual.mri            = structure, see FT_READ_MRI
%   cfg.individual.mesh           = structure, see FT_PREPARE_MESH
% You can specify the style with which the objects are displayed using
%   cfg.individual.headmodelstyle = 'vertex', 'edge', 'surface' or 'both' (default = 'edge')
%   cfg.individual.headshapestyle = 'vertex', 'edge', 'surface' or 'both' (default = 'vertex')
%
% The configuration structure should also contain the geometrical object of a
% template that serves as target
%   cfg.template.axes             = string, 'yes' or 'no (default = 'no')
%   cfg.template.elec             = structure, see FT_READ_SENS
%   cfg.template.grad             = structure, see FT_READ_SENS
%   cfg.template.opto             = structure, see FT_READ_SENS
%   cfg.template.headmodel        = structure, see FT_PREPARE_HEADMODEL
%   cfg.template.headshape        = structure, see FT_READ_HEADSHAPE
%   cfg.template.mri              = structure, see FT_READ_MRI
%   cfg.template.mesh             = structure, see FT_PREPARE_MESH
% You can specify the style with which the objects are displayed using
%   cfg.template.headmodelstyle   = 'vertex', 'edge', 'surface' or 'both' (default = 'edge')
%   cfg.template.headshapestyle   = 'vertex', 'edge', 'surface' or 'both' (default = 'vertex')
%
% You can specify one or multiple individual objects which will all be realigned and
% one or multiple template objects.
%
% See also FT_VOLUMEREALIGN, FT_ELECTRODEREALIGN, FT_DETERMINE_COORDSYS,
% FT_READ_SENS, FT_READ_HEADMODEL, FT_READ_HEADSHAPE

% Undocumented options, primarily to support FT_DEFACEVOLUME
%   cfg.showalpha
%   cfg.showlight
%   cfg.showapply
%   cfg.rotate
%   cfg.scale
%   cfg.translate
%   cfg.transformorder

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble provenance

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'required', {'individual', 'template'});
cfg.individual = ft_checkconfig(cfg.individual, 'renamed', {'vol', 'headmodel'});
cfg.individual = ft_checkconfig(cfg.individual, 'renamed', {'volstyle', 'headmodelstyle'});
cfg.template   = ft_checkconfig(cfg.template, 'renamed', {'vol', 'headmodel'});
cfg.template   = ft_checkconfig(cfg.template, 'renamed', {'volstyle', 'headmodelstyle'});

% allow to pass an initial transformation, this is used by FT_DEFACEVOLUME
cfg.rotate         = ft_getopt(cfg, 'rotate', [0 0 0]);
cfg.scale          = ft_getopt(cfg, 'scale', [1 1 1]);
cfg.translate      = ft_getopt(cfg, 'translate', [0 0 0]);
cfg.transformorder = ft_getopt(cfg, 'transformorder', {'rotate', 'translate', 'scale'});

% get the other options
cfg.showalpha                 = ft_getopt(cfg, 'showalpha'); % default is set below
cfg.showlight                 = ft_getopt(cfg, 'showlight', 'yes');
cfg.showapply                 = ft_getopt(cfg, 'showapply', 'yes');
cfg.unit                      = ft_getopt(cfg, 'unit', 'mm');
cfg.individual.elec           = ft_getopt(cfg.individual, 'elec', []);
cfg.individual.elecstyle      = ft_getopt(cfg.individual, 'elecstyle', {}); % key-value pairs
cfg.individual.grad           = ft_getopt(cfg.individual, 'grad', []);
cfg.individual.gradstyle      = ft_getopt(cfg.individual, 'gradstyle', {}); % key-value pairs
cfg.individual.opto           = ft_getopt(cfg.individual, 'opto', []);
cfg.individual.optostyle      = ft_getopt(cfg.individual, 'optostyle', {}); % key-value pairs
cfg.individual.headshape      = ft_getopt(cfg.individual, 'headshape', []);
if ~isempty(cfg.individual.headshape) && isfield(cfg.individual.headshape, 'tri')
  cfg.individual.headshapestyle = ft_getopt(cfg.individual, 'headshapestyle', 'surface');
else
  cfg.individual.headshapestyle = ft_getopt(cfg.individual, 'headshapestyle', 'vertex');
end
cfg.individual.headmodel      = ft_getopt(cfg.individual, 'headmodel', []);
cfg.individual.headmodelstyle = ft_getopt(cfg.individual, 'headmodelstyle', 'edge');
cfg.individual.mri            = ft_getopt(cfg.individual, 'mri', []);
cfg.individual.mristyle       = ft_getopt(cfg.individual, 'mristyle', {});
cfg.individual.mesh           = ft_getopt(cfg.individual, 'mesh', []);
cfg.individual.meshstyle      = ft_getopt(cfg.individual, 'meshstyle', {});

cfg.template.axes             = ft_getopt(cfg.template, 'axes', 'no');
cfg.template.elec             = ft_getopt(cfg.template, 'elec', []);
cfg.template.elecstyle        = ft_getopt(cfg.template, 'elecstyle', {}); % key-value pairs
cfg.template.grad             = ft_getopt(cfg.template, 'grad', []);
cfg.template.gradstyle        = ft_getopt(cfg.template, 'gradstyle', {}); % key-value pairs
cfg.template.opto             = ft_getopt(cfg.template, 'opto', []);
cfg.template.optostyle        = ft_getopt(cfg.template, 'optostyle', {}); % key-value pairs
cfg.template.headshape        = ft_getopt(cfg.template, 'headshape', []);
if ~isempty(cfg.template.headshape) && isfield(cfg.template.headshape, 'tri')
  cfg.template.headshapestyle   = ft_getopt(cfg.template, 'headshapestyle', 'surface');
else
  cfg.template.headshapestyle   = ft_getopt(cfg.template, 'headshapestyle', 'vertex');
end
cfg.template.headmodel        = ft_getopt(cfg.template, 'headmodel', []);
cfg.template.headmodelstyle   = ft_getopt(cfg.template, 'headmodelstyle', 'edge');
cfg.template.mri              = ft_getopt(cfg.template, 'mri', []);
cfg.template.mristyle         = ft_getopt(cfg.template, 'mristyle', {});
cfg.template.mesh             = ft_getopt(cfg.template, 'mesh', []);
cfg.template.meshstyle        = ft_getopt(cfg.template, 'meshstyle', {});

% convert the string that describes the style to a cell-array
cfg.template.headshapestyle   = updatestyle(cfg.template.headshapestyle);
cfg.individual.headshapestyle = updatestyle(cfg.individual.headshapestyle);
cfg.template.headmodelstyle   = updatestyle(cfg.template.headmodelstyle);
cfg.individual.headmodelstyle = updatestyle(cfg.individual.headmodelstyle);

% ensure that these are keyval cell-arrays, not cfg structures
for fn1={'individual', 'template'}
  fn1 = fn1{1};
  for fn2={'elecstyle', 'gradstyle', 'headshapestyle', 'headmodelstyle', 'mristyle', 'meshstyle'}
    fn2 = fn2{1};
    style = cfg.(fn1).(fn2);
    if isstruct(style)
      cfg.(fn1).(fn2) = ft_cfg2keyval(style);
    end
  end
end

if ~isempty(cfg.individual.headshape) && isfield(cfg.individual.headshape, 'color')
  if isfield(cfg.individual.headshape, 'tri') && size(cfg.individual.headshape.tri,1)==size(cfg.individual.headshape.color,1)
    cfg.individual.headshapestyle = ft_setopt(cfg.individual.headshapestyle, 'facecolor', cfg.individual.headshape.color);
  elseif size(cfg.individual.headshape.pos,1)==size(cfg.individual.headshape.color,1)
    cfg.individual.headshapestyle = ft_setopt(cfg.individual.headshapestyle, 'vertexcolor', cfg.individual.headshape.color);
  end
end

if ~isempty(cfg.template.headshape) && isfield(cfg.template.headshape, 'color')
  if isfield(cfg.template.headshape, 'tri') && size(cfg.template.headshape.tri,1)==size(cfg.template.headshape.color,1)
    cfg.template.headshapestyle = ft_setopt(cfg.template.headshapestyle, 'facecolor', cfg.template.headshape.color);
  elseif size(cfg.template.headshape.pos,1)==size(cfg.template.headshape.color,1)
    cfg.template.headshapestyle = ft_setopt(cfg.template.headshapestyle, 'vertexcolor', cfg.template.headshape.color);
  end
end

if isempty(cfg.showalpha)
  % only show the global alpha level if not specified for any of the individual objects
  cfg.showalpha = true;
  for fn1={'individual', 'template'}
    fn1 = fn1{1};
    for fn2={'elecstyle', 'gradstyle', 'headshapestyle', 'headmodelstyle', 'mristyle', 'meshstyle'}
      fn2 = fn2{1};
      style = cfg.(fn1).(fn2);
      cfg.showalpha = cfg.showalpha & ~any(endsWith(style(1:2:end), 'alpha'));
    end
  end
end

template   = struct(cfg.template);
individual = struct(cfg.individual);

% ensure that they are consistent with the latest FieldTrip version
for fn1={'individual', 'template'}
  fn1 = fn1{1};
  for fn2={'elec', 'grad', 'opto'}
    fn2 = fn2{1};
    if ~isempty(cfg.(fn1).(fn2))
      cfg.(fn1).(fn2) = ft_datatype_sens(cfg.(fn1).(fn2));
    end
  end
end

% ensure that they are consistent with the latest FieldTrip version
if ~isempty(template.headshape)
  template.headshape = fixpos(template.headshape);
end
if ~isempty(individual.headshape)
  individual.headshape = fixpos(individual.headshape);
end

% convert the coordinates of all geometrical objects into mm
fn = {'elec', 'grad', 'headshape', 'headmodel', 'mri', 'mesh'};
hasindividual = false(size(fn));
originalunit = cell(size(fn));
for i=1:length(fn)
  if ~isempty(individual.(fn{i}))
    hasindividual(i) = true;
    originalunit{i} = individual.(fn{i}).unit;
    individual.(fn{i}) = ft_convert_units(individual.(fn{i}), cfg.unit); % ensure that the units are known and all the same
  end
end
hastemplate = false(size(fn));
for i=1:length(fn)
  if ~isempty(template.(fn{i}))
    hastemplate(i) = true;
    template.(fn{i}) = ft_convert_units(template.(fn{i}), cfg.unit); % ensure that the units are known and all the same
  end
end

% determine the coordinate system of the template objects
coordsys = [];
for i=1:length(fn)
  if ~isempty(template.(fn{i}))
    if isfield(template.(fn{i}), 'coordsys')
      if isempty(coordsys)
        % remember the first coordinate system
        coordsys = template.(fn{i}).coordsys;
      end
      % ensure that all template objects have the same coordinate system as the first one
      assert(isequal(coordsys, template.(fn{i}).coordsys));
    end
  end
end

% ensure that the headshape surface is triangulated
if ~isempty(template.headshape)
  if ~isfield(template.headshape, 'tri') || isempty(template.headshape.tri)
    template.headshape.tri = projecttri(template.headshape.pos);
  end
end
if ~isempty(individual.headshape)
  if ~isfield(individual.headshape, 'tri') || isempty(individual.headshape.tri)
    individual.headshape.tri = projecttri(individual.headshape.pos);
  end
end

ft_info('Use the mouse to rotate the geometry, and click "redisplay" to update the light.');
ft_info('Close the figure when you are done.');

% open a figure
fig = figure;
set(fig, 'CloseRequestFcn',    @cb_quit);
set(fig, 'windowkeypressfcn',  @cb_keyboard);
set(gca, 'position', [0.05 0.15 0.75 0.75]);

% add the data and the settings to the figure
setappdata(fig, 'individual',     individual);
setappdata(fig, 'template',       template);
setappdata(fig, 'transform',      eye(4));
setappdata(fig, 'init',           true);
setappdata(fig, 'cleanup',        false);
setappdata(fig, 'rotate',         cfg.rotate);
setappdata(fig, 'scale',          cfg.scale);
setappdata(fig, 'translate',      cfg.translate);
setappdata(fig, 'coordsys',       coordsys); % can be unknown
setappdata(fig, 'unit',           cfg.unit);
setappdata(fig, 'labels',         true);
setappdata(fig, 'axes',           true);
setappdata(fig, 'grid',           true);
setappdata(fig, 'camlight',       true);
setappdata(fig, 'cfg',            cfg);

if isempty(coordsys) || strcmp(coordsys, 'unknown')
  ft_notice('the template coordinate system is unknown, selecting the viewpoint is not possible');
  ft_uilayout(gcf, 'tag', 'viewpointbtn', 'Visible', 'off');
else
  [labelx, labely, labelz] = coordsys2label(coordsys, 2, 0);
  ft_notice('the template coordinate system is "%s"', coordsys);
  ft_info('the positive X-axis is pointing to %s', labelx);
  ft_info('the positive Y-axis is pointing to %s', labely);
  ft_info('the positive Z-axis is pointing to %s', labelz);
end

% add the GUI elements
axmax = 150 * ft_scalingfactor('mm', cfg.unit);
axis([-axmax axmax -axmax axmax -axmax axmax]);
cb_creategui(gcf);
cb_redraw(gcf);
cb_help(cfg);
rotate3d on

cleanup = false;
while ~cleanup
  uiwait(fig);
  cfg.m   = getappdata(fig, 'transform');
  cleanup = getappdata(fig, 'cleanup');
end

% remember the details of the transformation
cfg.m         = getappdata(fig, 'transform');
cfg.rotate    = getappdata(fig, 'rotate');
cfg.scale     = getappdata(fig, 'scale');
cfg.translate = getappdata(fig, 'translate');

if istrue(cfg.showapply)
  % do not return the individual steps if the apply button was available and potentially pressed
  % in that case they only represent the last iteration, not the whole transformation
  cfg = removefields(cfg, {'rotate', 'scale', 'translate'});
end

delete(fig);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble provenance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_creategui(h, eventdata, handles)

% define the position of each GUI element
fig = getparent(h);

% get the initial values
rotate    = getappdata(fig, 'rotate');
scale     = getappdata(fig, 'scale');
translate = getappdata(fig, 'translate');
cfg       = getappdata(fig, 'cfg');

% constants
CONTROL_WIDTH   = 0.04;
CONTROL_HEIGHT  = 0.05;
CONTROL_HOFFSET = 0.75;
CONTROL_VOFFSET = 0.60;

% rotateui
uicontrol('tag', 'rotateui', 'parent',  fig, 'units', 'normalized', 'style', 'text', 'string', 'rotate',  'callback', [])
uicontrol('tag', 'rx',       'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', rotate(1), 'callback', @cb_redraw)
uicontrol('tag', 'ry',       'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', rotate(2), 'callback', @cb_redraw)
uicontrol('tag', 'rz',       'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', rotate(3), 'callback', @cb_redraw)
ft_uilayout(fig, 'tag', 'rotateui', 'width',  2*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET,                 'vpos',  CONTROL_VOFFSET+3*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'rx',       'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+3*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'ry',       'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+3*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'rz',       'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+3*CONTROL_HEIGHT);

% scaleui
uicontrol('tag', 'scaleui', 'parent',  fig, 'units', 'normalized', 'style', 'text', 'string', 'scale',  'callback', [])
uicontrol('tag', 'sx',      'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', scale(1), 'callback', @cb_redraw)
uicontrol('tag', 'sy',      'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', scale(2), 'callback', @cb_redraw)
uicontrol('tag', 'sz',      'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', scale(3), 'callback', @cb_redraw)
ft_uilayout(fig, 'tag', 'scaleui', 'width',  2*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET,                 'vpos',  CONTROL_VOFFSET+2*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'sx',      'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+2*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'sy',      'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+2*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'sz',      'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+2*CONTROL_HEIGHT);

% translateui
uicontrol('tag', 'translateui', 'parent',  fig, 'units', 'normalized', 'style', 'text', 'string', 'translate',  'callback', [])
uicontrol('tag', 'tx',          'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', translate(1), 'callback', @cb_redraw)
uicontrol('tag', 'ty',          'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', translate(2), 'callback', @cb_redraw)
uicontrol('tag', 'tz',          'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', translate(3), 'callback', @cb_redraw)
ft_uilayout(fig, 'tag', 'translateui', 'width',  2*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET,                 'vpos',  CONTROL_VOFFSET+1*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'tx',          'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+3*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+1*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'ty',          'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+4*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+1*CONTROL_HEIGHT);
ft_uilayout(fig, 'tag', 'tz',          'width',  CONTROL_WIDTH,   'height',  CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+5*CONTROL_WIDTH, 'vpos',  CONTROL_VOFFSET+1*CONTROL_HEIGHT);

% alpha ui
uicontrol('tag', 'alphaui', 'parent',  fig, 'units', 'normalized', 'style', 'text', 'string', 'alpha', 'value', [], 'callback', [], 'Visible', istrue(cfg.showalpha));
uicontrol('tag', 'alpha',   'parent',  fig, 'units', 'normalized', 'style', 'edit', 'string', '0.6',   'value', [], 'callback', @cb_redraw, 'Visible', istrue(cfg.showalpha));
ft_uilayout(fig, 'tag', 'alphaui',  'BackgroundColor', [0.8 0.8 0.8], 'width',  3*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-2*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'alpha',    'BackgroundColor', [0.8 0.8 0.8], 'width',  3*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-2*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET+3*CONTROL_WIDTH);

% control buttons
uicontrol('tag', 'viewpointbtn',  'parent',  fig, 'units', 'normalized', 'style', 'popup',      'string', 'top|bottom|left|right|front|back', 'value', 1, 'callback', @cb_viewpoint);
uicontrol('tag', 'camlightbtn',   'parent',  fig, 'units', 'normalized', 'style', 'popup',      'string', 'none|camlight|uniform',            'value', 1, 'callback', @cb_camlight, 'Visible', istrue(cfg.showlight));
uicontrol('tag', 'axesbtn',       'parent',  fig, 'units', 'normalized', 'style', 'checkbox',   'string', 'axes',         'value', getappdata(fig, 'axes'),     'callback', @cb_axes);
uicontrol('tag', 'labelsbtn',     'parent',  fig, 'units', 'normalized', 'style', 'checkbox',   'string', 'label',        'value', getappdata(fig, 'labels'),   'callback', @cb_labels);
uicontrol('tag', 'gridbtn',       'parent',  fig, 'units', 'normalized', 'style', 'checkbox',   'string', 'grid',         'value', getappdata(fig, 'grid'),     'callback', @cb_grid);
uicontrol('tag', 'redisplaybtn',  'parent',  fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'redisplay',    'value', [],                          'callback', @cb_redisplay);
uicontrol('tag', 'applybtn',      'parent',  fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'apply',        'value', [],                          'callback', @cb_apply,    'Visible', istrue(cfg.showapply));
uicontrol('tag', 'quitbtn',       'parent',  fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'quit',         'value', 1,                           'callback', @cb_quit);
ft_uilayout(fig, 'tag', 'viewpointbtn',                                     'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-3*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'camlightbtn',                                      'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-4*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'axesbtn',                                          'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-5*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'labelsbtn',                                        'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-6*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'gridbtn',                                          'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-7*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'redisplaybtn',   'BackgroundColor', [0.8 0.8 0.8], 'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-8*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'applybtn',       'BackgroundColor', [0.8 0.8 0.8], 'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-9*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);
ft_uilayout(fig, 'tag', 'quitbtn',        'BackgroundColor', [0.8 0.8 0.8], 'width',  6*CONTROL_WIDTH, 'height',  CONTROL_HEIGHT, 'vpos',  CONTROL_VOFFSET-10*CONTROL_HEIGHT, 'hpos',  CONTROL_HOFFSET);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_redraw(h, eventdata, handles)

fig            = getparent(h);
individual     = getappdata(fig, 'individual');
template       = getappdata(fig, 'template');
transform      = getappdata(fig, 'transform');
coordsys       = getappdata(fig, 'coordsys');
unit           = getappdata(fig, 'unit');
cfg            = getappdata(fig, 'cfg');

% disable all GUI elements
set(findall(fig, 'type', 'UIControl'), 'Enable', 'off')

% get the transformation details
rx = str2double(get(findobj(fig, 'tag', 'rx'), 'string'));
ry = str2double(get(findobj(fig, 'tag', 'ry'), 'string'));
rz = str2double(get(findobj(fig, 'tag', 'rz'), 'string'));
sx = str2double(get(findobj(fig, 'tag', 'sx'), 'string'));
sy = str2double(get(findobj(fig, 'tag', 'sy'), 'string'));
sz = str2double(get(findobj(fig, 'tag', 'sz'), 'string'));
tx = str2double(get(findobj(fig, 'tag', 'tx'), 'string'));
ty = str2double(get(findobj(fig, 'tag', 'ty'), 'string'));
tz = str2double(get(findobj(fig, 'tag', 'tz'), 'string'));

% create the transformation matrix
R = rotate   ([rx ry rz]);
S = scale    ([sx sy sz]);
T = translate([tx ty tz]);
% combine the transformation from the GUI with the one that has been previously applied
H = combine_transform(R, S, T, cfg.transformorder);
transform = H * transform;

axis vis3d
cla
hold on

if getappdata(fig, 'init')
  % start from a fresh view
  [az, el] = view(3);
  setappdata(fig, 'init', false);
else
  % remember the current view
  [az, el] = view();
  % FIXME it would be possible to move the plotting from all template objects here
  % FIXME so that they do not have to be plotted again upon each update
  % FIXME however, that requires that the "cla" above is implemented more specifically
end

% the "individual" struct is a local copy, so it is safe to change it here
if ~isempty(individual.headmodel)
  individual.headmodel = ft_transform_geometry(transform, individual.headmodel);
end
if ~isempty(individual.elec)
  individual.elec = ft_transform_geometry(transform, individual.elec);
end
if ~isempty(individual.grad)
  individual.grad = ft_transform_geometry(transform, individual.grad);
end
if ~isempty(individual.headshape)
  individual.headshape = ft_transform_geometry(transform, individual.headshape);
end
if ~isempty(individual.mri)
  individual.mri = ft_transform_geometry(transform, individual.mri);
end
if ~isempty(individual.mesh)
  individual.mesh = ft_transform_geometry(transform, individual.mesh);
end

% plot all the template and individual objects
if istrue(template.axes)
  ft_plot_axes([], 'unit', unit, 'coordsys', coordsys);
end

if ~isempty(template.mri)
  ft_plot_ortho(template.mri.anatomy, 'transform', template.mri.transform, 'unit', template.mri.unit, 'style', 'intersect', 'intersectmesh', individual.headshape, template.mristyle{:});
end

if ~isempty(individual.mri)
  ft_plot_ortho(individual.mri.anatomy, 'transform', individual.mri.transform, 'unit', individual.mri.unit, 'style', 'intersect', 'intersectmesh', template.headshape, individual.mristyle{:});
end

if ~isempty(template.elec)
  if isfield(template.elec, 'line')
    tmpbnd = [];
    tmpbnd.pos = template.elec.chanpos;
    tmpbnd.tri = template.elec.line;
    ft_plot_mesh(tmpbnd, 'vertexcolor', 'b', 'facecolor', 'none', 'edgecolor', 'b', 'vertexsize', 10)
  else
    ft_plot_sens(template.elec, template.elecstyle{:});
  end
end

if ~isempty(individual.elec)
  if isfield(individual.elec, 'line')
    tmpbnd = [];
    tmpbnd.pos = individual.elec.chanpos;
    tmpbnd.tri = individual.elec.line;
    ft_plot_mesh(tmpbnd, 'vertexcolor', 'r', 'facecolor', 'none', 'edgecolor', 'r', 'vertexsize', 10)
  else
    ft_plot_sens(individual.elec, individual.elecstyle{:});
  end
end

if ~isempty(template.grad)
  ft_plot_sens(template.grad, template.gradstyle{:});
end

if ~isempty(individual.grad)
  ft_plot_sens(individual.grad, individual.gradstyle{:});
end

if isstruct(template.headmodel)
  ft_plot_headmodel(template.headmodel, template.headmodelstyle{:});
end

if isstruct(individual.headmodel)
  ft_plot_headmodel(individual.headmodel, individual.headmodelstyle{:});
end

if isstruct(template.headshape) && isfield(template.headshape, 'pos') && ~isempty(template.headshape.pos)
  ft_plot_headshape(template.headshape, template.headshapestyle{:});
end

if isstruct(individual.headshape) && isfield(individual.headshape, 'pos') && ~isempty(individual.headshape.pos)
  ft_plot_headshape(individual.headshape, individual.headshapestyle{:})
end

if isstruct(template.mesh) && isfield(template.mesh, 'pos') && ~isempty(template.mesh.pos)
  % there can be multiple meshes as a struct-array
  ft_plot_mesh(template.mesh, template.meshstyle{:});
end

if isstruct(individual.mesh) && isfield(individual.mesh, 'pos') && ~isempty(individual.mesh.pos)
  % there can be multiple meshes as a struct-array
  ft_plot_mesh(individual.mesh, individual.meshstyle{:});
end

if istrue(cfg.showalpha)
  % only apply the global alpha level when it is not explicitly set for the individual objects
  alpha(str2double(get(findobj(fig, 'tag', 'alpha'), 'string')));
end

cb_camlight(h, []);
cb_axes(h, []);
cb_labels(h, []);
cb_grid(h, []);

% restore the current view
view(az, el);

% re-enable all GUI elements
set(findall(fig, 'type', 'UIControl'), 'Enable', 'on')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_camlight(h, eventdata)
fig = getparent(h);
delete(findall(fig, 'Type', 'light'));  % remove all existing lights
switch get(findall(fig, 'tag', 'camlightbtn'), 'value')
  case 1 % none
    lighting none
  case 2 % camlight
    lighting gouraud
    camlight
  case 3 % uniform
    lighting gouraud
    uniformlight
end
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_axes(h, eventdata)
fig = getparent(h);
if get(findall(fig, 'tag', 'axesbtn'), 'value')
  axis on
else
  axis off
end
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_labels(h, eventdata)
fig  = getparent(h);
unit = getappdata(fig, 'unit');
if get(findall(fig, 'tag', 'labelsbtn'), 'value')
  xlabel(sprintf('x (%s)', unit))
  ylabel(sprintf('y (%s)', unit))
  zlabel(sprintf('z (%s)', unit))
else
  xlabel('')
  ylabel('')
  zlabel('')
end
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_grid(h, eventdata)
fig = getparent(h);
if get(findall(fig, 'tag', 'gridbtn'), 'value')
  grid on
else
  grid off
end
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_apply(h, eventdata, handles)

fig       = getparent(h);
transform = getappdata(fig, 'transform');
cfg       = getappdata(fig, 'cfg');

% get the transformation details
rx = str2double(get(findobj(fig, 'tag', 'rx'), 'string'));
ry = str2double(get(findobj(fig, 'tag', 'ry'), 'string'));
rz = str2double(get(findobj(fig, 'tag', 'rz'), 'string'));
sx = str2double(get(findobj(fig, 'tag', 'sx'), 'string'));
sy = str2double(get(findobj(fig, 'tag', 'sy'), 'string'));
sz = str2double(get(findobj(fig, 'tag', 'sz'), 'string'));
tx = str2double(get(findobj(fig, 'tag', 'tx'), 'string'));
ty = str2double(get(findobj(fig, 'tag', 'ty'), 'string'));
tz = str2double(get(findobj(fig, 'tag', 'tz'), 'string'));

% create the transformation matrix
R = rotate   ([rx ry rz]);
S = scale    ([sx sy sz]);
T = translate([tx ty tz]);
% combine the transformation from the GUI with the one that has been previously applied
H = combine_transform(R, S, T, cfg.transformorder);
transform = H * transform;

set(findobj(fig, 'tag', 'rx'), 'string',  0);
set(findobj(fig, 'tag', 'ry'), 'string',  0);
set(findobj(fig, 'tag', 'rz'), 'string',  0);
set(findobj(fig, 'tag', 'sx'), 'string',  1);
set(findobj(fig, 'tag', 'sy'), 'string',  1);
set(findobj(fig, 'tag', 'sz'), 'string',  1);
set(findobj(fig, 'tag', 'tx'), 'string',  0);
set(findobj(fig, 'tag', 'ty'), 'string',  0);
set(findobj(fig, 'tag', 'tz'), 'string',  0);

setappdata(fig, 'transform',  transform);
setappdata(fig, 'rotate',     [rx ry rz]);
setappdata(fig, 'scale',      [sx sy sz]);
setappdata(fig, 'translate',  [tx ty tz]);

if ~getappdata(fig, 'cleanup')
  cb_redraw(h);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_help(h, eventdata)

fprintf('\n');
fprintf('==================================================================================\n');
fprintf('Press "h" to show this help.\n');
fprintf('Press "q" or close the window when you are done.\n');
fprintf('Press "v" to update the light position.\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_keyboard(h, eventdata)

fig = getparent(h);

if isempty(eventdata)
  % determine the key that corresponds to the uicontrol element that was activated
  key = get(fig, 'userdata');
else
  % determine the key that was pressed on the keyboard
  key = parsekeyboardevent(eventdata);
end

% get focus back to figure
if ~strcmp(get(h, 'type'), 'figure')
  set(h, 'enable', 'off');
  drawnow;
  set(h, 'enable', 'on');
end

if isempty(key)
  % this happens if you press the apple key
  key = '';
end

% the following code is largely shared by FT_SOURCEPLOT, FT_VOLUMEREALIGN, FT_INTERACTIVEREALIGN, FT_MESHREALIGN, FT_ELECTRODEPLACEMENT
switch key
  case {'' 'shift+shift' 'alt-alt' 'control+control' 'command-0'}
    % do nothing

  case 'q'
    cb_quit(h);

  case 'h'
    cb_help(h);

  case 'v' % camlight angle reset
    cb_camlight(h);

  otherwise
    % do nothing

end % switch key

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_redisplay(h, eventdata)
fig = getparent(h);
setappdata(fig, 'init', true);
cb_redraw(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_viewpoint(h, eventdata)

fig       = getparent(h);
coordsys  = getappdata(fig, 'coordsys');

% get the index of the option that was selected
val = get(h, 'value');
viewpoint = {'top', 'bottom', 'left', 'right', 'front', 'back'};
setviewpoint(gca, coordsys, viewpoint{val});

uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_quit(h, eventdata)
fig = getparent(h);
setappdata(fig, 'cleanup', true); % ensure that it will apply the current transformation
cb_apply(h);
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = getparent(h)
p = h;
while p~=0
  h = p;
  p = get(h, 'parent');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function style = updatestyle(style)
if ischar(style)
  switch style
    case 'vertex'
      style = {'vertexcolor', 'k', 'edgecolor', 'none', 'facecolor', 'none'};
    case 'edge'
      style = {'vertexcolor', 'none', 'edgecolor', 'k', 'facecolor', 'none'};
    case 'surface'
      style = {'vertexcolor', 'none', 'edgecolor', 'none', 'facecolor', 'skin'};
    case 'both'
      style = {'vertexcolor', 'none', 'edgecolor', 'k', 'facecolor', 'skin'};
    otherwise
      ft_error('unsupported style "%s"', style);
  end % switch
end % if

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip