text:
function [Vseq,Dseq] = eigenshuffle(Asequence)
% eigenshuffle: Consistent sorting for an eigenvalue/vector sequence
% [Vseq,Dseq] = eigenshuffle(Asequence)
%
% Includes munkres.m (by gracious permission from Yi Cao)
% to choose the appropriate permutation. This greatly
% enhances the speed of eigenshuffle over my previous
% release.
%
% http://www.mathworks.com/matlabcentral/fileexchange/20652
%
% Arguments: (input)
%  Asequence - an array of eigenvalue problems. If
%      Asequence is a 3-d numeric array, then each
%      plane of Asequence must contain a square
%      matrix that will be used to call eig.
%
%      Eig will be called on each of these matrices
%      to produce a series of eigenvalues/vectors,
%      one such set for each eigenvalue problem.
%
% Arguments: (Output)
%  Vseq - a 3-d array (pxpxn) of eigenvectors. Each
%      plane of the array will be sorted into a
%      consistent order with the other eigenvalue
%      problems. The ordering chosen will be one
%      that maximizes the energy of the consecutive
%      eigensystems relative to each other.
%
%  Dseq - pxn array of eigen values, sorted in order
%      to be consistent with each other and with the
%      eigenvectors in Vseq.
%
% Example:
%  Efun = @(t) [1 2*t+1 t^2 t^3;2*t+1 2-t t^2 1-t^3; ...
%               t^2 t^2 3-2*t t^2;t^3 1-t^3 t^2 4-3*t];
%
%  Aseq = zeros(4,4,21);
%  for i = 1:21
%    Aseq(:,:,i) = Efun((i-11)/10);
%  end
%  [Vseq,Dseq] = eigenshuffle(Aseq);
%
% To see that eigenshuffle has done its work correctly,
% look at the eigenvalues in sequence, after the shuffle.
%
% t = (-1:.1:1)';
% [t,Dseq']
% ans =
%        -1     8.4535           5      2.3447     0.20181
%      -0.9     7.8121      4.7687      2.3728     0.44644
%      -0.8     7.2481        4.56      2.3413     0.65054
%      -0.7     6.7524      4.3648      2.2709      0.8118
%      -0.6     6.3156      4.1751      2.1857     0.92364
%      -0.5     5.9283      3.9855      2.1118     0.97445
%      -0.4     5.5816      3.7931      2.0727     0.95254
%      -0.3     5.2676      3.5976      2.0768       0.858
%      -0.2     4.9791      3.3995      2.1156     0.70581
%      -0.1     4.7109         3.2      2.1742     0.51494
%         0     4.4605           3      2.2391     0.30037
%       0.1     4.2302         2.8      2.2971    0.072689
%       0.2     4.0303      2.5997      2.3303    -0.16034
%       0.3     3.8817      2.4047      2.3064    -0.39272
%       0.4     3.8108      2.1464      2.2628    -0.62001
%       0.5     3.8302      1.8986      2.1111    -0.83992
%       0.6     3.9301      1.5937      1.9298     -1.0537
%       0.7     4.0927      1.2308       1.745     -1.2685
%       0.8     4.3042     0.82515      1.5729     -1.5023
%       0.9     4.5572     0.40389      1.4272     -1.7883
%         1     4.8482  -8.0012e-16     1.3273     -2.1755
%
% Here, the columns are the shuffled eigenvalues.
% See that the second eigenvalue goes to zero, but
% the third eigenvalue remains positive. We can plot
% eigenvalues and see that they have crossed, near
% t = 0.35 in Efun.
%
% plot(-1:.1:1,Dseq')
%
% For a better appreciation of what eigenshuffle did,
% compare the result of eig directly on Efun(.3) and
% Efun(.4). Thus:
%
% [V3,D3] = eig(Efun(.3))
% V3 =
%     -0.74139      0.53464     -0.23551       0.3302
%      0.64781       0.4706     -0.16256      0.57659
%    0.0086542     -0.44236     -0.89119      0.10006
%     -0.17496     -0.54498      0.35197      0.74061
%
% D3 =
%     -0.39272            0            0            0
%            0       2.3064            0            0
%            0            0       2.4047            0
%            0            0            0       3.8817
%
% [V4,D4] = eig(Efun(.4))
% V4 =
%     -0.73026      0.19752      0.49743      0.42459
%      0.66202      0.21373      0.35297      0.62567
%     0.013412     -0.95225      0.25513      0.16717
%     -0.16815    -0.092308     -0.75026      0.63271
%
% D4 =
%     -0.62001            0            0            0
%            0       2.1464            0            0
%            0            0       2.2628            0
%            0            0            0       3.8108
%
% With no sort or shuffle applied, look at V3(:,3). See
% that it is really closest to V4(:,2), but with a sign
% flip. Since the signs on the eigenvectors are arbitrary,
% the sign is changed, and the most consistent sequence
% will be chosen. By way of comparison, see how the
% eigenvectors in Vseq have been shuffled, the signs
% swapped appropriately.
%
% Vseq(:,:,14)
% ans =
%       0.3302      0.23551     -0.53464      0.74139
%      0.57659      0.16256      -0.4706     -0.64781
%      0.10006      0.89119      0.44236   -0.0086542
%      0.74061     -0.35197      0.54498      0.17496
%
% Vseq(:,:,15)
% ans =
%      0.42459     -0.19752     -0.49743      0.73026
%      0.62567     -0.21373     -0.35297     -0.66202
%      0.16717      0.95225     -0.25513    -0.013412
%      0.63271     0.092308      0.75026      0.16815
%
% See also: eig
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 3.0
% Release date: 2/18/09

% Is Asequence a 3-d array?
Asize = size(Asequence);
if (Asize(1)~=Asize(2))
  error('Asequence must be a (pxpxn) array of eigen-problems, each of size pxp')
end
p = Asize(1);
if length(Asize)<3
  n = 1;
else
  n = Asize(3);
end

% the initial eigenvalues/vectors in nominal order
Vseq = zeros(p,p,n);
Dseq = zeros(p,n);
for i = 1:n
  [V,D] = eig(Asequence(:,:,i));
  D = diag(D);
  % initial ordering is purely in decreasing order.
  % If any are complex, the sort is in terms of the
  % real part.
  [junk,tags] = sort(real(D),1,'descend');

  Dseq(:,i) = D(tags);
  Vseq(:,:,i) = V(:,tags);
end

% was there only one eigenvalue problem?
if n < 2
  % we can quit now, having sorted the eigenvalues
  % as best as we could.
  return
end

% now, treat each eigenproblem in sequence (after
% the first one.)
for i = 2:n
  % compute distance between systems
  V1 = Vseq(:,:,i-1);
  V2 = Vseq(:,:,i);
  D1 = Dseq(:,i-1);
  D2 = Dseq(:,i);
  dist = (1-abs(V1'*V2)).*sqrt( ...
    distancematrix(real(D1),real(D2)).^2+ ...
    distancematrix(imag(D1),imag(D2)).^2);

  % Is there a best permutation? use munkres.
  % much faster than my own mintrace, munkres
  % is used by gracious permission from Yi Cao.
  reorder = munkres(dist);

  Vseq(:,:,i) = Vseq(:,reorder,i);
  Dseq(:,i) = Dseq(reorder,i);

  % also ensure the signs of each eigenvector pair
  % were consistent if possible
  S = squeeze(real(sum(Vseq(:,:,i-1).*Vseq(:,:,i),1))) < 0;
  Vseq(:,S,i) = -Vseq(:,S,i);
end

% =================
% end mainline
% =================
% begin subfunctions
% =================

function d = distancematrix(vec1,vec2)
% simple interpoint distance matrix
[vec1,vec2] = ndgrid(vec1,vec2);
d = abs(vec1 - vec2);

function [assignment,cost] = munkres(costMat)
% MUNKRES   Munkres (Hungarian) Algorithm for Linear Assignment Problem.
%
% [ASSIGN,COST] = munkres(COSTMAT) returns the optimal column indices,
% ASSIGN assigned to each row and the minimum COST based on the assignment
% problem represented by the COSTMAT, where the (i,j)th element represents the cost to assign the jth
% job to the ith worker.
%

% This is vectorized implementation of the algorithm. It is the fastest
% among all Matlab implementations of the algorithm.

% Examples
% Example 1: a 5 x 5 example
%{
[assignment,cost] = munkres(magic(5));
disp(assignment); % 3 2 1 5 4
disp(cost); %15
%}
% Example 2: 400 x 400 random data
%{
n=400;
A=rand(n);
tic
[a,b]=munkres(A);
toc                 % about 2 seconds
%}
% Example 3: rectangular assignment with inf costs
%{
A=rand(10,7);
A(A>0.7)=Inf;
[a,b]=munkres(A);
%}
% Reference:
% "Munkres' Assignment Algorithm, Modified for Rectangular Matrices",
% http://csclab.murraystate.edu/bob.pilgrim/445/munkres.html

% version 2.0 by Yi Cao at Cranfield University on 10th July 2008

assignment = zeros(1,size(costMat,1));
cost = 0;

costMat(costMat~=costMat)=Inf;
validMat = costMat<Inf;
validCol = any(validMat,1);
validRow = any(validMat,2);

nRows = sum(validRow);
nCols = sum(validCol);
n = max(nRows,nCols);
if ~n
    return
end

maxv=10*max(costMat(validMat));

dMat = zeros(n) + maxv;
dMat(1:nRows,1:nCols) = costMat(validRow,validCol);

%*************************************************
% Munkres' Assignment Algorithm starts here
%*************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   STEP 1: Subtract the row minimum from each row.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
minR = min(dMat,[],2);
minC = min(bsxfun(@minus, dMat, minR));

%**************************************************************************
%   STEP 2: Find a zero of dMat. If there are no starred zeros in its
%           column or row start the zero. Repeat for each zero
%**************************************************************************
zP = dMat == bsxfun(@plus, minC, minR);

starZ = zeros(n,1);
while any(zP(:))
    [r,c]=find(zP,1);
    starZ(r)=c;
    zP(r,:)=false;
    zP(:,c)=false;
end

while 1
%**************************************************************************
%   STEP 3: Cover each column with a starred zero. If all the columns are
%           covered then the matching is maximum
%**************************************************************************
    if all(starZ>0)
        break
    end
    coverColumn = false(1,n);
    coverColumn(starZ(starZ>0))=true;
    coverRow = false(n,1);
    primeZ = zeros(n,1);
    [rIdx, cIdx] = find(dMat(~coverRow,~coverColumn)==bsxfun(@plus,minR(~coverRow),minC(~coverColumn)));
    while 1
        %**************************************************************************
        %   STEP 4: Find a noncovered zero and prime it.  If there is no starred
        %           zero in the row containing this primed zero, Go to Step 5.
        %           Otherwise, cover this row and uncover the column containing
        %           the starred zero. Continue in this manner until there are no
        %           uncovered zeros left. Save the smallest uncovered value and
        %           Go to Step 6.
        %**************************************************************************
        cR = find(~coverRow);
        cC = find(~coverColumn);
        rIdx = cR(rIdx);
        cIdx = cC(cIdx);
        Step = 6;
        while ~isempty(cIdx)
            uZr = rIdx(1);
            uZc = cIdx(1);
            primeZ(uZr) = uZc;
            stz = starZ(uZr);
            if ~stz
                Step = 5;
                break;
            end
            coverRow(uZr) = true;
            coverColumn(stz) = false;
            z = rIdx==uZr;
            rIdx(z) = [];
            cIdx(z) = [];
            cR = find(~coverRow);
            z = dMat(~coverRow,stz) == minR(~coverRow) + minC(stz);
            rIdx = [rIdx(:);cR(z)];
            cIdx = [cIdx(:);stz(ones(sum(z),1))];
        end
        if Step == 6
            % *************************************************************************
            % STEP 6: Add the minimum uncovered value to every element of each covered
            %         row, and subtract it from every element of each uncovered column.
            %         Return to Step 4 without altering any stars, primes, or covered lines.
            %**************************************************************************
            [minval,rIdx,cIdx]=outerplus(dMat(~coverRow,~coverColumn),minR(~coverRow),minC(~coverColumn));
            minC(~coverColumn) = minC(~coverColumn) + minval;
            minR(coverRow) = minR(coverRow) - minval;
        else
            break
        end
    end
    %**************************************************************************
    % STEP 5:
    %  Construct a series of alternating primed and starred zeros as
    %  follows:
    %  Let Z0 represent the uncovered primed zero found in Step 4.
    %  Let Z1 denote the starred zero in the column of Z0 (if any).
    %  Let Z2 denote the primed zero in the row of Z1 (there will always
    %  be one).  Continue until the series terminates at a primed zero
    %  that has no starred zero in its column.  Unstar each starred
    %  zero of the series, star each primed zero of the series, erase
    %  all primes and uncover every line in the matrix.  Return to Step 3.
    %**************************************************************************
    rowZ1 = find(starZ==uZc);
    starZ(uZr)=uZc;
    while rowZ1>0
        starZ(rowZ1)=0;
        uZc = primeZ(rowZ1);
        uZr = rowZ1;
        rowZ1 = find(starZ==uZc);
        starZ(uZr)=uZc;
    end
end

% Cost of assignment
rowIdx = find(validRow);
colIdx = find(validCol);
starZ = starZ(1:nRows);
vIdx = starZ <= nCols;
assignment(rowIdx(vIdx)) = colIdx(starZ(vIdx));
cost = trace(costMat(assignment>0,assignment(assignment>0)));

function [minval,rIdx,cIdx]=outerplus(M,x,y)
[nx,ny]=size(M);
minval=inf;
for r=1:nx
    x1=x(r);
    for c=1:ny
        M(r,c)=M(r,c)-(x1+y(c));
        if minval>M(r,c)
            minval=M(r,c);
        end
    end
end
[rIdx,cIdx]=find(M==minval);


text:
function varargout = tan(varargin)
%TAN   Tangent of a DISKFUN.
%
% See also DISKFUN/SIN, DISKFUN/COS, DISKFUN/TANH

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

[varargout{1:nargout}] = tan@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
close all;
clearvars;
clc;
rng('default');

% Signal space
N = 1024;
% Number of measurements
M = 512;
% Sparsity level
K = 16;
% Number of signals
S = 1000;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
% Generate bi-uniform signals
X = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
PhiM = double(Phi);
% Measurement vectors
Y = Phi * X;

% Number of atoms per iteration
L = 4;
% Number of vectors in MMV set
T = 20;

fprintf('\nRunning Matlab-GOMP MMV...');
tic;
solver = spx.pursuit.joint.GOMP(PhiM, K);
solver.L = L;
begins = 1:T:S;
ends = begins+T-1;
for t=1:numel(begins)
    YY = Y(:, begins(t) : ends(t));
    result = solver.solve(YY);
end
t2=toc;

fprintf('\nRunning C GOMP MMV...');
tic;
res_norm_bnd = 0;
options.verbose = 0;
Z = spx.fast.gomp_mmv(Phi, Y, K, L, T, res_norm_bnd, options);
t1=toc;

% display summary  %
fprintf('\n\nSpeed summary for %d signals, dictionary size %d x %d:\n', S, M, N);
fprintf('Algorithm               Total time\n');
fprintf('--------------------------------------------------------\n');
fprintf('Matlab-GOMP MMV   %5.2f seconds\n', t1);
fprintf('C GOMP MMV        %5.2f seconds\n', t2);

author: indigits

repo: sparse-plex


text:
function map = yellowred(n)
%CHROMAJS returns a divergent lightyellow-red color map
%
%   Usage: m = yellowred([n])
%
%   Input parameters:
%       n - optional length of colormap (default uses the figure default length)
%
%   Output parameters:
%       m - colormap [n 3]
%
%   YELLOWRED(N) returns an N-by-3 matrix containing a divergent colormap.
%   Without a given N the same length as the current figure's colormap is used.
%   For details on the colormap have a look at: https://bit.ly/2vC3Ogr
%
%   To change the colormap current figure run: colormap(yellowred)
%
% See also: generate_colormap, moreland, plot_sound_field

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking input parameters =======================================
nargmin = 0;
nargmax = 1;
narginchk(nargmin,nargmax);
if nargin < nargmax
    n = size(get(gcf,'colormap'),1);
end

%% ===== Computation =====================================================
table = [ 255, 255, 224
255, 223, 184
255, 188, 148
255, 151, 119
255, 105,  98
238,  66,  86
210,  31,  71
176,   6,  44
139,   0,   0];
map = generate_colormap(table,n);

author: sfstoolbox

repo: sfs-matlab


text:
function dim_num = p00_dim_num ( prob )

%*****************************************************************************80
%
%% P00_DIM_NUM returns the spatial dimension for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem index.
%
%    Output, integer DIM_NUM, the spatial dimension of the
%    dependent variables.
%
  if ( prob == 1 )
    dim_num = p01_dim_num ( );
  elseif ( prob == 2 )
    dim_num = p02_dim_num ( );
  elseif ( prob == 3 )
    dim_num = p03_dim_num ( );
  elseif ( prob == 4 )
    dim_num = p04_dim_num ( );
  elseif ( prob == 5 )
    dim_num = p05_dim_num ( );
  elseif ( prob == 6 )
    dim_num = p06_dim_num ( );
  elseif ( prob == 7 )
    dim_num = p07_dim_num ( );
  elseif ( prob == 8 )
    dim_num = p08_dim_num ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_DIM_NUM - Fatal error!\n' );
    fprintf ( 1, '  Unexpected input value of PROB.\n' );
    error ( 'P00_DIM_NUM - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ons_secs = tapas_physio_preprocess_phys_timeseries(ons_secs, sqpar, doNormalize)
% Preprocess raw cardiac/respiratory time series (e.g. amplitude
% normalization, padding for acquisition duration)
%
%   ons_secs = tapas_physio_preprocess_phys_timeseries(ons_secs)
%
% IN
%   ons_secs    raw time series (c(ardiac), r(espiratory), t(ime), cpulse
%               (from file))
%
% OUT
%   ons_secs    padded and normalized time series for preprocessing
%
% EXAMPLE
%   tapas_physio_preprocess_phys_timeseries
%
%   See also

% Author: Lars Kasper
% Created: 2016-02-28
% Copyright (C) 2016 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

if nargin < 3
    doNormalize = 1;
end

hasCardiacData = ~isempty(ons_secs.c);
hasRespData = ~isempty(ons_secs.r);
hasDetectedCardiacPulses = ~isempty(ons_secs.cpulse);
hasAcquisitionCodes = ~isempty(ons_secs.acq_codes);

%% Normalize cardiac/respiratory time series to max 1

if doNormalize
    if hasCardiacData
        maxAbsC = max(abs(ons_secs.c));
        ons_secs.c_scaling = maxAbsC ;
        ons_secs.c = ons_secs.c/maxAbsC;
    end

    if hasRespData
        maxAbsR = max(abs(ons_secs.r));
        ons_secs.r_scaling = maxAbsR ;
        ons_secs.r = ons_secs.r/maxAbsR;
    end

end

% since resampling might have occured during read-in, dt is recalculated
ons_secs.dt = ons_secs.t(2) - ons_secs.t(1);

%% Padding of cardiac/respiratory time series, if start too late or duration too short

durationScan = sqpar.TR*(sqpar.Nscans+sqpar.Ndummies);

t = ons_secs.t;

tStartScan = 0;
if t(1) > tStartScan
    nSamples = ceil(-t(1)/ons_secs.dt);
    paddingStart = zeros(nSamples, 1);

    ons_secs.t = [(nSamples:-1:1)'*ons_secs.dt + t(1);t];

    if hasCardiacData
        ons_secs.c = [paddingStart; ons_secs.c];
    end

    if hasRespData
        ons_secs.r = [paddingStart; ons_secs.r];
    end

    if hasAcquisitionCodes
        ons_secs.acq_codes = [paddingStart; ons_secs.acq_codes];
    end
end

t = ons_secs.t;
durationPhysLog = t(end) - t(1);

if durationPhysLog < durationScan

    nSamples = ceil((durationScan - durationPhysLog)/ons_secs.dt);
    paddingEnd = zeros(nSamples, 1);

    ons_secs.t = [t; (1:nSamples)'*ons_secs.dt + t(end)];
    if hasCardiacData
        ons_secs.c = [ons_secs.c; paddingEnd];
    end

    if hasRespData
        ons_secs.r = [ons_secs.r; paddingEnd];
    end

    if hasAcquisitionCodes
        ons_secs.acq_codes = [ons_secs.acq_codes; paddingEnd];
    end
end

%% Remove onset time of log file to simplify plotting
tStartLog = t(1);
ons_secs.t_start = tStartLog;

ons_secs.t = ons_secs.t - tStartLog;

if hasDetectedCardiacPulses
    ons_secs.cpulse = ons_secs.cpulse - tStartLog;
end

author: translationalneuromodeling

repo: tapas


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (frontier)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% maximize - if 1, does max-product (not yet supported), else sum-product [0]
% filter -   if 1, do filtering, else smoothing [0]
%
% e.g., engine = enter_evidence(engine, ev, 'maximize', 1)

maximize = 0;
filter = 0;

% parse optional params
args = varargin;
nargs = length(args);
if nargs > 0
  for i=1:2:nargs
    switch args{i},
     case 'maximize', maximize = args{i+1};
     case 'filter', filter = args{i+1};
     otherwise,
      error(['invalid argument name ' args{i}]);
    end
  end
end

assert(~maximize);

[ss T] = size(evidence);
bnet = bnet_from_engine(engine);
onodes = find(~isemptycell(evidence));
cnodes = unroll_set(bnet.cnodes(:), ss, T);
pot_type = determine_pot_type(bnet, onodes);

CPDpot = convert_dbn_CPDs_to_pots(bnet, evidence, pot_type);

[engine.fwdback, loglik, engine.fwd_frontier, engine.back_frontier] = ...
    enter_soft_evidence(engine, CPDpot, onodes, pot_type, filter);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ xy, elem ] = triangulate_rectangle ( xl, xr, xn, yb, yt, yn, base )

%*****************************************************************************80
%
%% TRIANGULATE_RECTANGLE makes a triangular grid of a rectangle.
%
%  Discussion:
%
%    The rectangle is presumed to lie between (XL,YB) and (XR,YT).
%
%    We subdivide the rectangle into XN regions in the X direction,
%    and YN regions in the Y direction, and then split each quadrilateral,
%    creating 2 * XN * YN triangular elements.
%
%    The locations of the NN=(XN+1)*(YN+1) nodes are stored in the 2 by NN
%    real array XY.
%
%    The triples of node indices that make up each triangle are stored int
%    the 3 by NE integer array ELEM.  Here, NE = 2 * XN * YN, and the
%    nodes are stored in counterclockwise order.  The node indices will
%    begin at 0 if the input quantity BASE is set to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XL, XR, the left and right X limits.
%
%    Input, integer XN, the number of elements along the X direction.
%
%    Input, real YB, YT, the bottom and top Y limits.
%
%    Input, integer YN, the number of elements along the Y direction.
%
%    Input, integer BASE, the indexing base:
%    0, the first node has index 0.
%    1, the first node has index 1.
%
  if ( nargin < 7 )
    base = 1;
  end
%
%  Generate 1D data.
%
  x1d = linspace ( xl, xr, xn + 1 );
  y1d = linspace ( yb, yt, yn + 1 );
%
%  Create (X,Y) table.
%  This can be done faster, using meshgrid, but then you lose control
%  of ordering and arranging.
%
  xyn = ( xn + 1 ) * ( yn + 1 );
  xy = zeros ( 2, xyn );

  k = 0;
  for j = 1 : yn + 1
    for i = 1 : xn + 1
      k = k + 1;
      xy(1,k) = x1d(i);
      xy(2,k) = y1d(j);
    end
  end
%
%  Create connectivity.
%
  en = 2 * xn * yn;
  elem = zeros ( 3, en );
  e = 0;

  for j = 1 : yn
    for i = 1 : xn
      sw = ( j - 1 ) * ( xn + 1 ) + i;
      se = ( j - 1 ) * ( xn + 1 ) + i + 1;
      ne = ( j     ) * ( xn + 1 ) + i + 1;
      nw = ( j     ) * ( xn + 1 ) + i;
      e = e + 1;
      elem(1,e) = sw;
      elem(2,e) = se;
      elem(3,e) = nw;
      e = e + 1;
      elem(1,e) = ne;
      elem(2,e) = nw;
      elem(3,e) = se;
    end
  end
%
%  Shift base if requested.
%
  if ( base == 0 )
    elem = elem - 1;
  end

  return
end


text:
function [  ] = DngToJpg( dngFile,qual,gamma )
%DNGTOJPG Summary of this function goes here
%   Detailed explanation goes here

    [raw, info] = Load_Data_and_Metadata_from_DNG(dngFile);
    if exist('gamma','var')
        srgb = run_pipeline(raw,info,'raw','srgb');
        srgb=srgb.^gamma;
        imwrite(srgb, [dngFile(1:end-4),'_gamma.JPG'], 'quality', qual);
    else
        srgb = run_pipeline(raw,info,'raw','tone');
        imwrite(srgb, [dngFile(1:end-4), '.JPG'], 'quality', qual);
    end
end

author: AbdoKamel

repo: sidd-ground-truth-image-estimation


text:
function [GB] = bit2GB(bit)
% Convert computery things from bits to gigabytes.
% Chad A. Greene 2012
GB = bit*2^-33;

author: Sable

repo: mcbench-benchmarks


text:
function corr_coef = RegreesionLIVE(objectiveValues,mos)
%this script is used to calculate the pearson linear correlation
%coefficient and root mean sqaured error after regression

%get the objective scores computed by the IQA metric and the subjective
%scores provided by the dataset
% matData = load('VSIOnLIVE.mat');
% VSIOnLIVE = matData.VSIOnLIVE;
% objectiveValues = VSIOnLIVE(:,1);
% mos = VSIOnLIVE(:,2);

%plot objective-subjective score pairs
% p = plot(objectiveValues,mos,'+');
% set(p,'Color','blue','LineWidth',1);

%initialize the parameters used by the nonlinear fitting function
beta(1) = max(mos);
beta(2) = min(mos);
beta(3) = mean(objectiveValues);
beta(4) = 0.1;
beta(5) = 40;

%fitting a curve using the data
[bayta ehat,J] = nlinfit(objectiveValues,mos,@logistic,beta);
%given a ssim value, predict the correspoing mos (ypre) using the fitted curve
[ypre junk] = nlpredci(@logistic,objectiveValues,bayta,ehat,J);

RMSE = sqrt(sum((ypre - mos).^2) / length(mos));%root meas squared error
corr_coef = corr(mos, ypre, 'type','Pearson'); %pearson linear coefficient

%draw the fitted curve
% t = min(objectiveValues):0.01:max(objectiveValues);
% [ypre junk] = nlpredci(@logistic,t,bayta,ehat,J);
% hold on;
% p = plot(t,ypre);
% set(p,'Color','black','LineWidth',2);
% legend('Images in LIVE','Curve fitted with logistic function', 'Location','NorthEast');
% xlabel('Objective score by VSI');
% ylabel('MOS');

author: HuiZeng


text:
function [iT] =  inverse_field(T)
    % inverse displacement field
    % size(T) = [N, M, 2] or
    % size(T) = [N, M, K, 3]
    %
    sz = size(T);
    if ndims(T) == 3
        T(:,:, 1) = imgaussfilt(T(:,:,1), 1.3);
        T(:,:, 2) = imgaussfilt(T(:,:,2), 1.3);

        [n1, n2] = ndgrid(1:sz(1), 1:sz(2));
        fB1 = griddata(T(:,:,1) + n1, T(:,:,2) + n2, n1, n1, n2, 'cubic') - n1;
        fB2 = griddata(T(:,:,1) + n1, T(:,:,2) + n2, n2, n1, n2, 'cubic') - n2;
        iT = cat(3, fB1, fB2);

%         [n1, n2] = ndgrid(1:sz(1), 1:sz(2));
%         F = scatteredInterpolant(fl(T(:,:,1) + n1), fl(T(:,:,2) + n2), fl(n1), 'natural', 'nearest');
%         fB1 = F(n1, n2) - n1;
%         F = scatteredInterpolant(fl(T(:,:,1) + n1), fl(T(:,:,2) + n2), fl(n2), 'natural', 'nearest');
%         fB2 = F(n1, n2) - n2;
%         iT = cat(3, fB1, fB2);
    elseif false%ndims(T) == 4 && size(T,4) == 3 && isa(T, 'double')
        sgm = 1.3;
%         sgm = 0.5;
        T = cat(4, imgaussfilt3(T(:,:,:, 1), sgm), imgaussfilt3(T(:,:,:, 2), sgm), imgaussfilt3(T(:,:,:, 3), sgm));
        iT = mex_inverse_3d_displacements_double(T, 1);
    elseif ndims(T) == 4
        fB = cell(3, 1);
        n = cell(3, 1);
        [n{1}, n{2}, n{3}] = ndgrid(1:sz(1), 1:sz(2), 1:sz(3));
        parfor i = 1:3
%             tic
%             F = scatteredInterpolant(fl(T(:,:,:,1) + n{1}), fl(T(:,:,:,2) + n{2}), fl(T(:,:,:,3) + n{3}), fl(n{i}));
%             tmp = F({1:sz(1), 1:sz(2), 1:sz(3)}) - n{i};
%             toc
%             imagesc(tmp(:,:, 4)); pause;
%             tic
            fB{i} = griddata(T(:,:,:,1) + n{1}, T(:,:,:,2) + n{2}, T(:,:,:,3) + n{3}, ...
                n{i}, n{1}, n{2}, n{3}, 'linear') - n{i};
%             toc
%             imagesc(fB{i}(:,:,4)); pause;
%             imagesc(abs(fB{i}(:,:,4) - tmp(:,:,4))); colorbar; pause;
        end

%         fB1 = griddata(T(:,:,:,1) + n1, T(:,:,:,2) + n2, T(:,:,:,3) + n3, ...
%             n1, n1, n2, n3, 'linear') - n1;
%         fB2 = griddata(T(:,:,:,1) + n1, T(:,:,:,2) + n2, T(:,:,:,3) + n3, ...
%             n2, n1, n2, n3, 'linear') - n2;
%         fB3 = griddata(T(:,:,:,1) + n1, T(:,:,:,2) + n2, T(:,:,:,3) + n3, ...
%             n3, n1, n2, n3, 'linear') - n3;
        iT = cat(4, fB{1}, fB{2}, fB{3});
    end
    iT(isnan(iT))=0;
end

repo: pTVreg


text:
function display(obj)
% Display a NIFTI-1 object
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

disp(' ');
disp([inputname(1),' = '])
disp(' ');
disp(obj)
disp(' ')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef SuperEllipsePrinter < handle

    properties (Access = private)
        mesh
        outPutName
        iteration
    end

    methods (Access = public)

        function obj = SuperEllipsePrinter(cParams)
            obj.init(cParams)
        end

        function print(obj)
            dI.mesh            = obj.mesh;
            dI.outName         = obj.outPutName;
            dI.pdim            = '2D';
            dI.ptype           = 'MICRO';
            ps = PostProcessDataBaseCreator(dI);
            dB = ps.getValue();
            dB.printers = 'Density';
            postCase = 'Density';
            postProcess = Postprocess(postCase,dB);
            d.x = ones(size(obj.mesh.coord(:,1),1),1);
            postProcess.print(obj.iteration,d);
        end

        function captureImage(obj)
            f = obj.outPutName;
            it = obj.iteration;
            imagePath = '/home/alex/git-repos/MicroStructurePaper/';
            outPutNameImage = fullfile(imagePath,[f,num2str(it)]);
            inputFileName = fullfile('Output',f,[f,num2str(it),'.flavia.res']);
            sI.fileName = f;
            sI.outPutImageName = outPutNameImage;
            sI.inputFileName = inputFileName;
            imageCapturer = GiDImageCapturer(sI);
            imageCapturer.capture();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.outPutName = cParams.outPutName;
            obj.mesh       = cParams.mesh;
            obj.iteration  = 0;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% patchMarchDistMapIterative
% Below is a demonstration of the features of the |patchMarchDistMapIterative| function

%%
clear; close all; clc;

%% Syntax
% |[varargout]=patchMarchDistMapIterative(varargin);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function out = spm_shoot_norm(job)
% Spatially normalise and smooth fMRI/PET data to MNI space, using Shoot deformation fields
% FORMAT out = spm_shoot_norm(job)
% job - a structure generated by the configuration file
%   job.template - Shoot template for aligning to MNI space. Aligns to population
%                  average if no template is provided.
%   job.subj(n)  - Subject n
%       subj(n).def    - Shoot deformation field
%       subj(n).images - Images for this subject
%   job.vox      - Voxel sizes for spatially normalised images
%   job.bb       - Bounding box for spatially normalised images
%   job.preserve - How to transform
%                  0 = preserve concentrations
%                  1 = preserve integral (cf "modulation")
%
% Normally, Shoot generates deformations that align with the average-
% shaped template.  This routine includes the option to compose the
% shoot deformations with an affine transform derived from an affine
% registration of the template (the final one generated by Shoot),
% with the TPM data released with SPM.
%
% Note that trilinear interpolation is used, and no masking is done.  It
% is therefore essential that the images are realigned and resliced
% before they are spatially normalised.  Alternatively, contrast images
% generated from unsmoothed native-space fMRI/PET data can be spatially
% normalised for a 2nd level analysis.
%
% Two "preserve" options are provided.  One of them should do the
% equavalent of generating smoothed "modulated" spatially normalised
% images.  The other does the equivalent of smoothing the modulated
% normalised fMRI/PET, and dividing by the smoothed Jacobian determinants.
%
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_shoot_norm.m 7496 2018-11-23 11:14:43Z john $

% Hard coded stuff, that should maybe be customisable
tpm  = fullfile(spm('Dir'),'tpm','TPM.nii');
Mmni = spm_get_space(tpm);

% Shoot template
if ~isempty(job.template{1})
    Nt     = nifti(job.template{1});
    do_aff = true;
else
    Nt     = nifti(tpm);
    do_aff = false;
end

% Deal with desired bounding box and voxel sizes.
%--------------------------------------------------------------------------
bb   = job.bb;
vox  = job.vox;
Mt   = Nt.mat;
dimt = size(Nt.dat);

if any(isfinite(bb(:))) || any(isfinite(vox))
    [bb0,vox0] = spm_get_bbox(Nt, 'old');

    msk = ~isfinite(vox); vox(msk) = vox0(msk);
    msk = ~isfinite(bb);   bb(msk) =  bb0(msk);

    bb  = sort(bb);
    vox = abs(vox);

    % Adjust bounding box slightly - so it rounds to closest voxel.
    bb(:,1) = round(bb(:,1)/vox(1))*vox(1);
    bb(:,2) = round(bb(:,2)/vox(2))*vox(2);
    bb(:,3) = round(bb(:,3)/vox(3))*vox(3);
    dim = round(diff(bb)./vox+1);
    of  = -vox.*(round(-bb(1,:)./vox)+1);
    mat = [vox(1) 0 0 of(1) ; 0 vox(2) 0 of(2) ; 0 0 vox(3) of(3) ; 0 0 0 1];
    if det(Mt(1:3,1:3)) < 0
        mat = mat*[-1 0 0 dim(1)+1; 0 1 0 0; 0 0 1 0; 0 0 0 1];
    end
else
    dim = dimt(1:3);
    mat = Mt;
end

if isfield(job.data,'subj') || isfield(job.data,'subjs')
    if do_aff
        [pth,nam,ext] = fileparts(Nt.dat.fname); %#ok
        if exist(fullfile(pth,[nam '_2mni.mat']),'file')
            load(fullfile(pth,[nam '_2mni.mat']),'mni');
        else
            % Affine registration of Shoot Template with MNI space.
            %--------------------------------------------------------------
            fprintf('** Affine registering "%s" with MNI space **\n', nam);
            clear mni
            mni.affine = Mmni/spm_klaff(Nt,tpm,1);
            mni.code   = 'MNI152';
            save(fullfile(pth,[nam '_2mni.mat']),'mni', spm_get_defaults('mat.format'));
        end
        M = mat\mni.affine/Mt;
        mat_intent = mni.code;
    else
        M = mat\eye(4);
        mat_intent = 'Aligned';
    end
    fprintf('\n');

    if isfield(job.data,'subjs')
        % Re-order data
        %------------------------------------------------------------------
        subjs = job.data.subjs;
        subj  = struct('deformation',cell(numel(subjs.deformations),1),...
                       'images',   cell(numel(subjs.deformations),1));
        for i=1:numel(subj)
            subj(i).deformation = {subjs.deformations{i}};
            subj(i).images   = cell(numel(subjs.images),1);
            for j=1:numel(subjs.images)
                subj(i).images{j} = subjs.images{j}{i};
            end
        end
    else
        subj = job.data.subj;
    end

    % Loop over subjects
    %----------------------------------------------------------------------
    out = cell(1,numel(subj));
    for i=1:numel(subj)
        % Spatially normalise data from this subject
        [pth,nam,ext] = fileparts(subj(i).deformation{1}); %#ok
        fprintf('** "%s" **\n', nam);
        out{i} = deal_with_subject(subj(i).deformation,subj(i).images, mat,dim,M,job.preserve,job.fwhm,mat_intent);
    end

    if isfield(job.data,'subjs')
        out1 = out;
        out  = cell(numel(subj),numel(subjs.images));
        for i=1:numel(subj)
            for j=1:numel(subjs.images)
                out{i,j} = out1{i}{j};
            end
        end
    end
end
%__________________________________________________________________________

%__________________________________________________________________________
function out = deal_with_subject(Py,PI,mat,dim,M,jactransf,fwhm,mat_intent)
% Py         - Filename of shoot deformation.
% PI         - Filenames of images to spatially normalise.
% mat        - Voxel-to-world matrix for header of warped images.
% dim        - Dimensions of warped images.
% M          - Matrix for adjusting the deformation fields (eg when using
%              different bounding boxes or origins, or when an additional
%              affine transform is included.
% jactransf  - Whether or not to "modulate".
% fwhm       - FWHM for smoothing.
% mat_intent - Something for the header to indicate average space or MNI
%              space.
%
% Generate deformation, which is the inverse of the usual one (it is for "pushing"
% rather than the usual "pulling"). This deformation is affine transformed to
% allow for different voxel sizes and bounding boxes, and also to incorporate
% the affine mapping between MNI space and the population average shape.
%--------------------------------------------------------------------------
NY  = nifti(Py{1});
M   = M*NY.mat;
y   = single(squeeze(NY.dat(:,:,:,:,:)));
d   = size(y);

if norm(M-eye(4))>1e-3 || any(d(1:3) ~= dim(1:3))
    y0  = zeros([dim(1:3),3],'single');
    for d=1:3
        yd = y(:,:,:,d);
        for x3=1:dim(3)
            y0(:,:,x3,d) = single(spm_slice_vol(yd,M\spm_matrix([0 0 x3]),dim(1:2),[1 NaN]));
        end
    end
else
    y0 = y;
end

odm = zeros(1,3);
oM  = zeros(4,4);
out = cell(1,numel(PI));
for m=1:numel(PI)

    % Generate headers etc for output images
    %----------------------------------------------------------------------
    [pth,nam,ext,num] = spm_fileparts(PI{m}); %#ok
    NI = nifti(fullfile(pth,[nam ext]));
    NO = NI;
    if jactransf
        if fwhm==0
            NO.dat.fname=fullfile(pth,['mw' nam ext]);
        else
            NO.dat.fname=fullfile(pth,['smw' nam ext]);
        end
        NO.dat.scl_slope = 1.0;
        NO.dat.scl_inter = 0.0;
        NO.dat.dtype     = 'float32-le';
    else
        if fwhm==0
            NO.dat.fname=fullfile(pth,['w' nam ext]);
        else
            NO.dat.fname=fullfile(pth,['sw' nam ext]);
        end
    end
    NO.dat.dim = [dim NI.dat.dim(4:end)];
    NO.mat  = mat;
    NO.mat0 = mat;
    NO.mat_intent  = mat_intent;
    NO.mat0_intent = mat_intent;

    if strcmp(mat_intent,'Aligned'), to_what = 'Average';  else
                                     to_what = mat_intent; end
    if fwhm==0
        NO.descrip = ['Shoot-normed (' to_what ')'];
    else
        NO.descrip = sprintf('Smoothed (%g) Shoot normed (%s)',fwhm, to_what);
    end
    out{m} = NO.dat.fname;
    NO.extras = [];
    create(NO);

    % Smoothing settings
    vx  = sqrt(sum(mat(1:3,1:3).^2));
    krn = max(fwhm./vx,0.1);

    % Loop over volumes within the file
    %----------------------------------------------------------------------
    fprintf('%s',nam); drawnow;
    for j=1:size(NI.dat,4)

        % Check if it is an "imported" image to normalise
        if sum(sum((NI.mat - NY.mat ).^2)) < 0.0001 && ...
           sum(sum((NI.mat - NI.mat0).^2)) > 0.0001
            % No affine transform necessary
            M0 = NI.mat0;
        else
            % Need to resample the mapping by an affine transform
            % so that it maps from voxels in the native space image
            % to voxels in the spatially normalised image.
            %--------------------------------------------------------------
            M0 = NI.mat;
            if ~isempty(NI.extras) && isstruct(NI.extras) && isfield(NI.extras,'mat')
                M1 = NI.extras.mat;
                if size(M1,3) >= j && sum(sum(M1(:,:,j).^2)) ~=0
                    M0 = M1(:,:,j);
                end
            end
        end
        M   = inv(M0);
        dm  = [size(NI.dat),1,1,1,1];
        if ~all(dm(1:3)==odm) || ~all(M(:)==oM(:))
            % Generate new deformation (if needed)
            y   = zeros(size(y0),'single');
            y(:,:,:,1) = M(1,1)*y0(:,:,:,1) + M(1,2)*y0(:,:,:,2) + M(1,3)*y0(:,:,:,3) + M(1,4);
            y(:,:,:,2) = M(2,1)*y0(:,:,:,1) + M(2,2)*y0(:,:,:,2) + M(2,3)*y0(:,:,:,3) + M(2,4);
            y(:,:,:,3) = M(3,1)*y0(:,:,:,1) + M(3,2)*y0(:,:,:,2) + M(3,3)*y0(:,:,:,3) + M(3,4);

            % Generate Jacobian determinants.
            c          = abs(spm_diffeo('jacdet',y)*(det(NI.mat(1:3,1:3))/det(NO.mat(1:3,1:3))));
            c(:,:,[1 end]) = NaN; % Boundary voxels are not handled well - so remove
            c(:,[1 end],:) = NaN;
            c([1 end],:,:) = NaN;
        end
        odm = dm(1:3);
        oM  = M;

        % Write the warped data for this time point.
        %------------------------------------------------------------------
        for k=1:size(NI.dat,5)
            for l=1:size(NI.dat,6)
                f  = single(NI.dat(:,:,:,j,k,l));
                if ~jactransf
                    % Unmodulated
                    f = spm_diffeo('pull',f,y).*c;
                    if fwhm>0
                        spm_smooth(f,f,krn); % Side effects
                        cs = zeros(size(c),'like',c);
                        spm_smooth(c,cs,krn); % Side effects
                        f  = f./max(cs,0.000001);
                    else
                        f  = f./max(c ,0.000001);
                    end
                else
                    % Modulated, by pushing
                    f = spm_diffeo('pull',f,y).*c;
                    spm_smooth(f,f,krn); % Side effects
                end
                NO.dat(:,:,:,j,k,l) = f;
                fprintf('\t%d,%d,%d', j,k,l); drawnow;
            end
        end
    end
    fprintf('\n'); drawnow;
end
%__________________________________________________________________________

author: spm


text:
function [Yl,Yh,Yscale] = dtwavexfm(X,nlevels,biort,qshift);

% Function to perform a n-level DTCWT decompostion on a 1-D column vector X
% (or on the columns of a matrix X).
%
% [Yl,Yh,Yscale] = dtwavexfm(X,nlevels,biort,qshift);
%
%     X -> real 1-D signal column vector (or matrix of vectors)
%
%     nlevels -> No. of levels of wavelet decomposition
%
%     biort ->  'antonini'   => Antonini 9,7 tap filters.
%               'legall'     => LeGall 5,3 tap filters.
%               'near_sym_a' => Near-Symmetric 5,7 tap filters.
%               'near_sym_b' => Near-Symmetric 13,19 tap filters.
%
%     qshift -> 'qshift_06' => Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters,
%                              (only 6,6 non-zero taps).
%               'qshift_a' =>  Q-shift 10,10 tap filters,
%                              (with 10,10 non-zero taps, unlike qshift_06).
%               'qshift_b' => Q-Shift 14,14 tap filters.
%               'qshift_c' => Q-Shift 16,16 tap filters.
%               'qshift_d' => Q-Shift 18,18 tap filters.
%
%
%     Yl     -> The real lowpass subband from the final level.
%     Yh     -> A cell array containing the complex highpass subband for each level.
%     Yscale -> This is an OPTIONAL output argument, that is a cell array containing
%               the real lowpass coefficients at every scale.
%
%
% Example: [Yl,Yh] = dtwavexfm(X,5,'near_sym_b','qshift_b');
% performs a 5-level transform on the real image X using the 13,19-tap filters
% for level 1 and the Q-shift 14-tap filters for levels >= 2.
%
% Nick Kingsbury and Cian Shaffrey
% Cambridge University, May 2002

if isstr(biort) & isstr(qshift)		%Check if the inputs are strings
   biort_exist = exist([biort '.mat']);
   qshift_exist = exist([qshift '.mat']);
   if biort_exist == 2 & qshift_exist == 2;  %Check to see if the filters exist as .mat files
      load (biort);
      load (qshift);
   else
      error('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTWAVEXFM for details.');
   end
else
   error('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTWAVEXFM.');
end

L = size(X);

if any(rem(L(1),2)),	 % ensure that X is an even length, thus enabling it to be extended if needs be.
   error('Size of X must be a multiple of 2');
end

if nlevels == 0, return; end

%initialise
Yh=cell(nlevels,1);
if nargout == 3
   Yscale=cell(nlevels,1);   % This is only required if the user specifies a third output component.
end

j = sqrt(-1);

% Level 1.
Hi = colfilter(X, h1o);
Lo = colfilter(X, h0o);
t = 1:2:size(Hi,1);
Yh{1} = Hi(t,:) + j*Hi(t+1,:); % Convert Hi to complex form.
if nargout == 3
   Yscale{1} = Lo;
end

if nlevels >= 2;  % Levels 2 and above.
   for level = 2:nlevels;
      if rem(size(Lo,1),4),	% Check to see if height of Lo is divisable by 4, if not extend.
         Lo = [Lo(1,:); Lo; Lo(end,:)];
      end
      Hi = coldfilt(Lo,h1b,h1a);
      Lo = coldfilt(Lo,h0b,h0a);
	   t = 1:2:size(Hi,1);
   	Yh{level} = Hi(t,:) + j*Hi(t+1,:); % Convert Hi to complex form.
   	if nargout == 3
      	Yscale{level} = Lo;
   	end
   end
end

Yl = Lo;

return

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function M = makeFslXfmMatrix(T,R,S,filename)
%MAKEFSLXFMMATRIX Make FSL-compatible transformation matrix.
%   M = MAKEFSLXFMMATRIX(T,R,S,FILENAME) outputs a 4x4 transformation
%   matrix performing the translations in T, the rotations in R and the
%   scalings in S and writes this matrix to the file specified by the
%   string FILENAME. This file is compatible with FSL's flirt and thus can
%   be directly used without further modification.
%
%   This function is actually the inverse function of FSL's avscale in that
%   avscale reads a transformation matrix file and lists the corresponding
%   translation, rotation and so on.
%
%   MAKEFSLXFMMATRIX can be especially useful if you would like to run
%   simulations.
%
%   Please note that, just as in avscale, T is in millimeters, R is in
%   radians and S is unitless.
%
%   Usage example:
%   ==============
%   T = [1.433440 19.715600 0.786690]; % in mm.
%   R = [0.102735 -0.155470 0.121564]; % in rad.
%   S = [1 1 1]; % unitless
%   filename = 'thisIsWhatIWant.mat'; % :))
%
%   M = makeFslXfmMatrix(T,R,S,filename);
%
%   Now you can use avscale to see if things are working correctly.
%
%   [status,result] = system(['avscale --allparams thisIsWhatIWant.mat']);
%   disp(result)
%
%   Best.

% Check inputs
if ~isnumeric(T) || ~isnumeric(R) || ~isnumeric(S)
    error('T, R and S have to be numeric...');
end

if ~ischar(filename)
    error('FILENAME has to be a string...');
end

% Get rotation angles about each axis.
thetaX = R(1);
thetaY = R(2);
thetaZ = R(3);

% Compute the rotation matrices.
Rx = [1 0 0; 0 cos(thetaX) sin(thetaX);0 -sin(thetaX) cos(thetaX)];
Ry = [cos(thetaY) 0 -sin(thetaY);0 1 0;sin(thetaY) 0 cos(thetaY)];
Rz = [cos(thetaZ) sin(thetaZ) 0;-sin(thetaZ) cos(thetaZ) 0;0 0 1];

% Concatenate rotations.
R_3x3 = Rx*Ry*Rz;

% Put scalings on the diagonal.
S_3x3 = diag(S);

% Make the 4x4 transformation matrix.
M = [R_3x3*S_3x3 T(:); 0 0 0 1];

% Try to open the file
fid = fopen(filename,'w');

% Check if it worked
if fid == -1
    error(['Could not open for writing: ' filename]);
end

for k = 1:4
    fprintf(fid,[num2str(M(k,:)) '\n']);
end

fclose(fid);


text:
report_this_filefun(mfilename('fullpath'));

m = [];
m =0:0.1:max(newt2.Magnitude);
m2 =min(newt2.Magnitude):0.1:max(newt2.Magnitude);
m3 = 10.^(-(m2-min(newt2.Magnitude)));
m3 = m3*newt2.Count;

k = 0:0.1:min(newt2.Magnitude)-0.1;
m4 = k*8+newt2.Count;
m = [m4 m3];
clf
%plot(0:0.1:max(newt2.Magnitude),log10(m))
%grid
%hold on
newcat = newt2;

lepo = length(m) -1;
mm = [newt2.Count m(1:lepo) ];
bval = mm-m;
bvalfl = bval(length(bval):-1:1);
maxmag = max(newcat.Magnitude);
mima = min(newcat.Magnitude);
if mima > 0 ; mima = 0 ; end

% number of mag units
nmagu = (maxmag*10)+1;

bvalsum = zeros(1,nmagu);
bvalsum3 = zeros(1,nmagu);

bvalsum = m;
bvalsum3 = m(length(m):-1:1);
xt3 = (maxmag:-0.1:mima);

backg_be = log10(bvalsum);
backg_ab = log10(bvalsum3);
orient tall
rect = [0.2,  0.3, 0.70, 0.6];           % plot Freq-Mag curves
axes('position',rect);

semilogy(xt3,bvalsum3,'-.m')
hold on
semilogy(xt3,bvalsum3,'om')
difb = [0 diff(bvalsum3) ];
semilogy(xt3,difb,'xg')
semilogy(xt3,difb,'g')
grid

% Marks the point of maximum curvature
%
i = find(difb == max(difb));
te = semilogy(xt3(i),difb(i),'xk');
set(te,'LineWidth',2,'MarkerSize',ms10)
te = semilogy(xt3(i),bvalsum3(i),'xk');
set(te,'LineWidth',2,'MarkerSize',ms10)

% Estimate the b-value
%
i2 = round(i/3);
te = semilogy(xt3(i2),difb(i2),'xk');
set(te,'LineWidth',2,'MarkerSize',ms10)
te = semilogy(xt3(i2),bvalsum3(i2),'xk');
set(te,'LineWidth',2,'MarkerSize',ms10)

xlabel('Magnitude','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
ylabel('Cumulative Number','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
set(gca,'Color',[1 1 0.6])
set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on')

par2 = 0.1 * max(bvalsum3);
par3 = 0.12 * max(bvalsum3);
M1b = [];
M1b = [xt3(i) bvalsum3(i)];
tt3=num2str(fix(100*M1b(1))/100);
text( M1b(1),M1b(2),['|: M1=',tt3] )

M2b = [];
M2b =  [xt3(i2) bvalsum3(i2)];
tt4=num2str(fix(100*M2b(1))/100);
text( M2b(1),M2b(2),['|: M2=',tt4] )

ll = xt3 > M1b(1) & xt3 < M2b(1);
x = xt3(ll);

l = newcat.Magnitude > M1b(1) & newcat.Magnitude < M2b(1);
me = 0.4343/(sum(bval.*xt3)/(sum(bval))-M1b(1));
mer = 1.96*me/(sqrt(length(newcat(l,6))));

so = log10(bval(10*M1b(1)+2)) - log10(bval(10*M2b(1)));
me= so/( M2b(1)-0.2- M1b(1));

pause(0.1)

y = backg_ab(ll);
[p,s] = polyfit(x,y,1);                   % fit a line to background
f = polyval(p,x);
f = 10.^f;
hold on
ttm= semilogy(x,f,'b');                         % plot linear fit to backg
set(ttm,'LineWidth',2)
r = corrcoef(x,y);
r = r(1,2);
std_backg = std(y - polyval(p,x));      % standard deviation of fit

p=-p(1,1);
p=fix(100*p)/100;
std_backg=fix(100*std_backg)/100;
tt2=num2str(std_backg);
tt1=num2str(p);

n = length(x)+3;
l = b(:,6) > M1b(1) & b(:,6) <= M2b(1);
les = (mean(b(l,6)) - (M1b(1)+0.05))/0.1;

mysofu = @(x)sofu(x,n,les);
so = fzero(mysofu,1.0);
me2 = log(so)/(-2.3026*0.1)

rect=[0 0 1 1];
h2=axes('position',rect);
set(h2,'visible','off');

txt1=text(.16, .18,['B-Value(L2): ',tt1,'  B(eff): ',num2str(me), '  B(mean2)= ',num2str(me2)]);
set(txt1,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
txt1=text(.16, .1,['Standard Deviation: ',tt2]);
set(txt1,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
set(gcf,'visible','on');
zmap_message_center.set_info('  ','Done')
done

author: CelsoReyes


text:
function [minC, maxC] = coor_extremals(C)
% usage
%   [minC, maxC] = coor_extremals(C)
%
% input
%   C = cell matrix of nested cell/numeric matrices, with numeric or empty
%       matrices at the bottom nesting level, of coordinates. These belong
%       to different graphics objects.
%
% output
%   minC = minimum coordinate over all graphics objects
%   maxC = maximum coordinate over all graphics objects
%
% See also PLOT_SCALINGS, MIN_CELL, MAX_CELL.

% depends
%   min_cell, max_cell

maxC = max_cell(C);
minC = min_cell(C);


text:
function [varargout] = isequaln(varargin)

% ISEQUALN is a drop-in replacement for the same function that was
% introduced in MATLAB R2012a.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% see https://github.com/fieldtrip/fieldtrip/issues/899

alternatives = which(mfilename, '-all');
if ~iscell(alternatives)
  % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
  alternatives = {alternatives};
end

keep = true(size(alternatives));
for i=1:numel(alternatives)
  keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
  keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
end
alternatives = alternatives(keep);

if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
  % remove this directory from the path
  p = fileparts(mfilename('fullpath'));
  sw = warning('off','backtrace');
  warning('Removing "%s" from your path.\nSee http://www.fieldtriptoolbox.org/faq/should_i_add_fieldtrip_with_all_subdirectories_to_my_matlab_path/', p);
  warning(sw);
  rmpath(p);
  % call the original MATLAB function
  if exist(mfilename, 'builtin')
    [varargout{1:nargout}] = builtin(mfilename, varargin{:});
  else
    [varargout{1:nargout}] = feval(mfilename, varargin{:});
  end
  return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = isequaln(varargin)

tf = isequalwithequalnans(varargin{:});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};


text:
function phi = rad(phi)
%RAD returns the given angle in radians
%
%   Usage: phi = rad(phi)
%
%   Input parameters:
%       phi     - angle, can be a scalar or matrix / degree
%
%   Output parameters:
%       phi     - angle / rad
%
%   See also: deg

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Computation =====================================================
phi = phi./180*pi;


text:
%==============================================================================
% Copyright (C) 2005, Jan Modersitzki and Nils Papenberg, see copyright.m,
% this file is part of the FLIRT Package, all rights reserved;
% http://www.math.uni-luebeck.de/SAFIR/FLIRT-MATLAB.html
%==============================================================================
% function [B,Bstr] = getElasticMatrixStg(Omega,m,mu,lambda);
% generates the elastic matrix B for the domain Omega with resolution m,
% by default, mu = 1, lambda = 0;
%==============================================================================
function [B,Bstr] = getElasticMatrixStg(Omega,m,mu,lambda)

Bstr = 'elastic-stg';
if ~exist('mu','var'),     mu     = 1;  end;
if ~exist('lambda','var'), lambda = 0;  end;
h = Omega./m;

dim = length(Omega);

switch dim
  case 2
    d11 = spdiags(ones(m(1),1)*[-1,1],[0,1],m(1),m(1)+1)/h(1);
    d12 = spdiags(ones(m(2)-1,1)*[-1,1],[0,1],m(2)-1,m(2))/h(2);
    d21 = spdiags(ones(m(1)-1,1)*[-1,1],[0,1],m(1)-1,m(1))/h(1);
    d22 = spdiags(ones(m(2),1)*[-1,1],[0,1],m(2),m(2)+1)/h(2);

    D11 = sparse(kron(speye(m(2)),d11));
    D12 = sparse(kron(d12,speye(m(1)+1)));
    D21 = sparse(kron(speye(m(2)+1),d21));
    D22 = sparse(kron(d22,speye(m(1))));

    % build the elastic operator
    %
    %      | \nabla 0      |
    %  B = |               |
    %      | 0      \nabla |
    %      |               |
    %      | Div1   Div2   |
    % B[U] = [\partial_1 U1,\partial_2 U1, \partial_1 U2,\partial_2 U2,div(U)]

    a  = sqrt(mu);
    b  = sqrt(mu+lambda);
    n1 = (m(1)+1)*m(2);
    n2 = m(1)*(m(2)+1);
    j1 = size(D11,1);
    j2 = size(D12,1);
    j3 = size(D21,1);
    j4 = size(D22,1);

    B = [  a*D11,sparse(j1,n2);
      a*D12,sparse(j2,n2);
      sparse(j3,n1),a*D21;
      sparse(j4,n1),a*D22;
      b*D11,b*D22];
    return;

  otherwise
    error('nyi');
end;
%==============================================================================

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ivec = s_to_i4vec ( s, n )

%*****************************************************************************80
%
%% S_TO_I4VEC reads a integer vector from a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character S(*), the string to be read.
%
%    Input, integer N, the number of values expected.
%
%    Output, integer IVEC(N), the values read from the string.
%
  i = 0;
  ivec = [];

  ilo = 1;

  while ( i < n )

    i = i + 1;

    [ ivec(i), lchar ] = s_to_i4 ( s(ilo:end) );

    ilo = ilo + lchar;

  end

  return
end


text:
function [cfg, M] = ft_sourcemovie(cfg, source, source2)

% FT_SOURCEMOVIE displays the source reconstruction on a cortical mesh
% and allows the user to scroll through time with a movie.
%
% Use as
%   ft_sourcemovie(cfg, source)
% where the input source data is obtained from FT_SOURCEANALYSIS, or a
% a parcellated source structure (i.e. contains a brainordinate field) and
% cfg is a configuration structure that should contain
%
%   cfg.funparameter    = string, functional parameter that is color coded (default = 'avg.pow')
%   cfg.maskparameter   = string, functional parameter that is used for opacity (default = [])
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_SOURCEPLOT, FT_SOURCEINTERPOLATE, FT_SOURCEPARCELLATE

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the initial part deals with parsing the input options and data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
if ft_abort
  return
end
ft_preamble debug
ft_preamble loadvar source
ft_preamble provenance source

% ensure that the input data is valid for this function, this will also do
% backward-compatibility conversions of old data that for example was
% read from an old *.mat file
source = ft_checkdata(source, 'datatype', 'source', 'feedback', 'yes');

% get the options
xlim          = ft_getopt(cfg, 'time');
ylim          = ft_getopt(cfg, 'frequency');
zlim          = ft_getopt(cfg, 'funcolorlim');
olim          = ft_getopt(cfg, 'opacitylim');                     % don't use alim as variable name
cfg.xparam    = ft_getopt(cfg, 'xparam', 'time');                 % use time as default
cfg.yparam    = ft_getopt(cfg, 'yparam');                         % default is dealt with below
cfg.funparameter  = ft_getopt(cfg, 'funparameter', 'avg.pow');        % use power as default
cfg.funcolormap   = ft_getopt(cfg, 'funcolormap',  'jet');
cfg.maskparameter = ft_getopt(cfg, 'maskparameter');

if isempty(cfg.yparam) && isfield(source, 'freq')
  % the default is freq (if present)
  cfg.yparam = 'freq';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the actual computation is done in the middle part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==2
  fun = getsubfield(source, cfg.funparameter); % might be avg.pow
elseif nargin>2 && isfield(source2, 'pos')
  % in this case two conditions of data seem to have been inputted
  fun  = getsubfield(source, cfg.funparameter); % might be avg.pow
  fun2 = getsubfield(source2, cfg.funparameter);
elseif nargin>2
  % assume the first data argument to be a parcellation, and the second a
  % parcellated structure

  % THIS DOES NOT EXIST ANYMORE: IF PARCELLATED DATA IS IN THE INPUT, IT
  % WILL BE UNPARCELLATED AUTOMATICALLY BY FT_CHECKDATA
  ft_error('bogus');
end
if size(source.pos)~=size(fun,1)
  ft_error('inconsistent number of vertices in the cortical mesh');
end

if ~isfield(source, 'tri')
  ft_error('source.tri missing, this function requires a triangulated cortical sheet as source model');
end

if ~isempty(cfg.maskparameter) && ischar(cfg.maskparameter)
  mask = double(getsubfield(source, cfg.maskparameter));
else
  mask = 0.5*ones(size(fun));
end

xparam = source.(cfg.xparam);
if length(xparam)~=size(fun,2)
  error('inconsistent size of "%s" compared to "%s"', cfg.funparameter, cfg.xparam);
end

if ~isempty(cfg.yparam)
  yparam = source.(cfg.yparam);
  if length(cfg.yparam)~=size(fun,3)
    error('inconsistent size of "%s" compared to "%s"', cfg.funparameter, cfg.yparam);
  end
end

if isempty(xlim)
  xlim(1) = min(xparam);
  xlim(2) = max(xparam);
end

xbeg = nearest(xparam, xlim(1));
xend = nearest(xparam, xlim(2));
% update the configuration
cfg.xlim = xparam([xbeg xend]);

if ~isempty(cfg.yparam)
  if isempty(ylim)
    ylim(1) = min(yparam);
    ylim(2) = max(yparam);
  end
  ybeg = nearest(yparam, ylim(1));
  yend = nearest(yparam, ylim(2));
  % update the configuration
  cfg.ylim = yparam([ybeg yend]);
  hasyparam = true;
else
  % this allows us not to worry about the yparam any more
  yparam = nan;
  ybeg = 1;
  yend = 1;
  cfg.ylim = [];
  hasyparam = false;
end

% make a subselection of the data
xparam  = xparam(xbeg:xend);
yparam  = yparam(ybeg:yend);
fun     = fun(:,xbeg:xend,ybeg:yend);
if exist('fun2', 'var')
  fun2 = fun2(:,xbeg:xend,ybeg:yend);
end
mask    = mask(:,xbeg:xend,ybeg:yend);
clear xbeg xend ybeg yend

if isempty(zlim)
  zlim(1) = min(fun(:));
  zlim(2) = max(fun(:));
  % update the configuration
  cfg.funcolorlim = zlim;
end

if isempty(olim)
  olim(1) = min(mask(:));
  olim(2) = max(mask(:));
  if olim(1)==olim(2)
    olim(1) = 0;
    olim(2) = 1;
  end
  % update the configuration
  cfg.opacitylim = olim;
end

% collect the data and the options to be used in the figure
opt.cfg     = cfg;
opt.xparam  = xparam;
opt.yparam  = yparam;
opt.xval    = 0;
opt.yval    = 0;
opt.dat     = fun;
opt.mask    = mask;
opt.pos     = source.pos;
opt.tri     = source.tri;
opt.vindx   = source.inside(:);
opt.speed   = 1;
opt.record  = 0;
opt.threshold = 0;
opt.frame   = 0;
opt.cleanup = false;
if isfield(source, 'parcellation')
  opt.parcellation = source.parcellation;
  opt.parcellationlabel = source.parcellationlabel;
end

% add functional data of optional third input to the opt structure
% FIXME here we should first check whether the meshes correspond!
if nargin>2 && isfield(source2, 'pos')
  opt.dat2 = fun2;
end

% get a handle to a figure
h  = figure;
set(h, 'color', [1 1 1]);
set(h, 'toolbar', 'figure');
set(h, 'visible', 'on');
set(h, 'CloseRequestFcn', @cb_quitbutton);
set(h, 'position', [100 50 1000 600]);
set(h, 'windowbuttondownfcn', @cb_getposition);

% get timer object
t = timer;
set(t, 'timerfcn', {@cb_timer, h}, 'period', 0.1, 'executionmode', 'fixedSpacing');

% make the user interface elements
cambutton    = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'light', 'userdata', 'C');
playbutton   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'play',   'userdata', 'p');
recordbutton = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'record', 'userdata', 'r');
quitbutton   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'quit',   'userdata', 'q');

thrminmin    = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'downarrow');
thrminplus   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'shift+downarrow');
thr          = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'threshold', 'userdata', 't');
thrplusmin   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'shift+uparrow');
thrplusplus  = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'uparrow');

spdmin       = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '<', 'userdata', 'shift+downarrow');
spd          = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'speed','userdata', 's');
spdplus      = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '>', 'userdata', 'shift+uparrow');
climminmin   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'leftarrow');
climminplus  = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'shift+leftarrow');
clim         = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'colorlim', 'userdata', 'z');
climplusmin  = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'shift+rightarrow');
climplusplus = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'rightarrow');
sliderx      = uicontrol('parent', h, 'units', 'normalized', 'style', 'slider',     'string', sprintf('%s = ', cfg.xparam));
stringx      = uicontrol('parent', h, 'units', 'normalized', 'style', 'text');
slidery      = uicontrol('parent', h, 'units', 'normalized', 'style', 'slider',     'string', sprintf('%s = ', cfg.yparam));
stringy      = uicontrol('parent', h, 'units', 'normalized', 'style', 'text');
stringz      = uicontrol('parent', h, 'units', 'normalized', 'style', 'text');

set(cambutton,    'position', [0.095 0.28  0.09 0.05 ], 'callback', @cb_keyboard);
set(quitbutton,   'position', [0.005 0.28  0.09 0.05 ], 'callback', @cb_keyboard);
set(playbutton,   'position', [0.005 0.22  0.09 0.05 ], 'callback', @cb_keyboard);
set(recordbutton, 'position', [0.095 0.22  0.09 0.05 ], 'callback', @cb_keyboard);
set(thrminmin,    'position', [0.005 0.16  0.03 0.025], 'callback', @cb_keyboard);
set(thrminplus,   'position', [0.005 0.185 0.03 0.025], 'callback', @cb_keyboard);
set(thr,          'position', [0.035 0.16  0.12 0.05 ], 'callback', @cb_keyboard);
set(thrplusmin,   'position', [0.155 0.16  0.03 0.025], 'callback', @cb_keyboard);
set(thrplusplus,  'position', [0.155 0.185 0.03 0.025], 'callback', @cb_keyboard);
set(climminmin,   'position', [0.005 0.10  0.03 0.025], 'callback', @cb_keyboard);
set(climminplus,  'position', [0.005 0.125 0.03 0.025], 'callback', @cb_keyboard);
set(clim,         'position', [0.035 0.10  0.12 0.05 ], 'callback', @cb_keyboard);
set(climplusmin,  'position', [0.155 0.10  0.03 0.025], 'callback', @cb_keyboard);
set(climplusplus, 'position', [0.155 0.125 0.03 0.025], 'callback', @cb_keyboard);
set(spdmin,       'position', [0.005 0.04  0.03 0.05 ], 'callback', @cb_keyboard);
set(spd,          'position', [0.035 0.04  0.12 0.05 ], 'callback', @cb_keyboard);
set(spdplus,      'position', [0.155 0.04  0.03 0.05 ], 'callback', @cb_keyboard);
set(sliderx,      'position', [0.02  0.35  0.38 0.03 ], 'callback', @cb_slider);
set(slidery,      'position', [0.200 0.005 0.78 0.03 ], 'callback', @cb_slider);
set(stringx,      'position', [0.750 0.90  0.18 0.03 ]);
set(stringy,      'position', [0.750 0.85  0.18 0.03 ]);
set(stringz,      'position', [0.60  0.95  0.33 0.03 ]);

set(stringx, 'string', sprintf('%s = ', cfg.xparam));
set(stringy, 'string', sprintf('%s = ', cfg.yparam));
set(stringz, 'string', sprintf('location = '));
set(stringx, 'horizontalalignment', 'right', 'backgroundcolor', [1 1 1]);
set(stringy, 'horizontalalignment', 'right', 'backgroundcolor', [1 1 1]);
set(stringz, 'horizontalalignment', 'right', 'backgroundcolor', [1 1 1]);

% create axes object to contain the mesh
hx = axes;
set(hx, 'position', [0.4 0.08 0.6 0.8]);
set(hx, 'tag', 'mesh');
hs = ft_plot_mesh(source, 'edgecolor', 'none', 'facecolor', [0.5 0.5 0.5], 'vertexcolor', 0.*opt.dat(:,1,1), 'facealpha', 0.*opt.mask(:,1,1), 'clim', [0 1], 'alphalim', [0 1], 'alphamap', 'rampup', 'colormap', cfg.funcolormap, 'maskstyle', 'colormix');

lighting gouraud
material dull
cam1 = camlight('left');
cam2 = camlight('right');
caxis(cfg.funcolorlim);

% create axis object to contain a time course
hy = axes;
set(hy, 'position', [0.02 0.45 0.38 0.5]);
set(hy, 'yaxislocation', 'right');

if ~hasyparam
  tline = plot(opt.xparam, mean(opt.dat(opt.vindx,:))); hold on;
  abc = axis;
  axis([opt.xparam(1) opt.xparam(end) abc(3:4)]);
  vline = plot(opt.xparam(1)*[1 1], abc(3:4), 'r');
  hline1 = plot(opt.xparam([1 end]), [0 0], 'k');
  hline2 = plot(opt.xparam([1 end]), [0 0], 'k');

  if nargin>2 && isfield(source2, 'pos')
    tline2 = plot(opt.xparam, mean(opt.dat2(opt.vindx,:)), 'r'); hold on;
  end

else
  error('not yet implemented');
end
set(hy, 'tag', 'timecourse');

% remember the various handles
opt.h   = h;  % handle to the figure
opt.hs  = hs; % handle to the mesh
opt.hx  = hx; % handle to the axes containing the mesh
opt.hy  = hy; % handle to the axes containing the timecourse
opt.cam = [cam1 cam2]; % handles to the light objects
opt.vline = vline; % handle to the vertical line in the ERF plot
opt.tline = tline; % handle to the ERF
if nargin>2 && isfield(source2, 'pos')
  opt.tline2 = tline2;
end
opt.hline1 = hline1; % handle for the horizontal line for upper threshold
opt.hline2 = hline2; % handle for hte horizontal line for lower threshold
opt.playbutton   = playbutton; % handle to the playbutton
opt.recordbutton = recordbutton; % handle to the recordbutton
opt.quitbutton   = quitbutton; % handle to the quitbutton
opt.threshold    = [0 0];

opt.t   = t;
opt.sliderx  = sliderx;
opt.slidery  = slidery;
opt.stringx  = stringx;
opt.stringy  = stringy;
opt.stringz  = stringz;

if ~hasyparam
  set(opt.slidery, 'visible', 'off');
  set(opt.stringy, 'visible', 'off');
end

setappdata(h, 'opt', opt);

while opt.cleanup==0
  uiwait(h);
  opt = getappdata(h, 'opt');
end
stop(opt.t);

if nargout
  M = opt.movie;
end

delete(h);

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

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

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

persistent previous_valx previous_valy previous_vindx previous_clim previous_thr

if isempty(previous_valx)
  previous_valx = 0;
end
if isempty(previous_valy)
  previous_valy = 0;
end
if isempty(previous_clim)
  previous_clim = [0 1];
end
if isempty(previous_thr)
  previous_thr = [0 0];
end

h    = getparent(h);
opt  = getappdata(h, 'opt');
valx = get(opt.sliderx, 'value');
valx = round(valx*(size(opt.dat,2)-1))+1;
valx = min(valx, size(opt.dat,2));
valx = max(valx, 1);

valy = get(opt.slidery, 'value');
valy = round(valy*(size(opt.dat,3)-1))+1;
valy = min(valy, size(opt.dat,3));
valy = max(valy, 1);

mask = squeeze(opt.mask(:,valx,valy));

% update stuff
if previous_valx~=valx || previous_valy~=valy || ~isequal(previous_clim, opt.cfg.funcolorlim) || ~isequal(previous_thr, opt.threshold)
  % update strings
  set(opt.stringx, 'string', sprintf('%s = %3.3f\n', opt.cfg.xparam, opt.xparam(valx)));
  set(opt.stringy, 'string', sprintf('%s = %3.3f\n', opt.cfg.yparam, opt.yparam(valy)));

  dat = opt.dat(:,valx,valy);
  mask(dat>opt.threshold(1)&dat<opt.threshold(2)) = 0;

  % convert the color-data + opacity into rgb for robust rendering
  bgcolor = repmat([0.5 0.5 0.5], [numel(mask) 1]);
  rgb     = bg_rgba2rgb(bgcolor, dat, opt.cfg.funcolormap, opt.cfg.funcolorlim, mask, 'rampup', opt.cfg.opacitylim);

  % update data in mesh
  set(opt.hs, 'FaceVertexCData', rgb, 'facecolor', 'interp');

  set(opt.vline, 'xdata', [1 1]*opt.xparam(valx));
end

% update ERF-plot
set(opt.hy,    'ylim',   opt.cfg.funcolorlim);
set(opt.vline, 'ydata',  opt.cfg.funcolorlim);
if ~(numel(previous_vindx)==numel(opt.vindx) && all(previous_vindx==opt.vindx))
  tmp = mean(opt.dat(opt.vindx,:,valy),1);
  set(opt.tline, 'ydata', tmp);
  %set(opt.hy,    'ylim',  [min(tmp(:)) max(tmp(:))]);
  %set(opt.vline, 'ydata', [min(tmp(:)) max(tmp(:))]);

  if isfield(opt, 'dat2')
    tmp = mean(opt.dat2(opt.vindx,:,valy),1);
    set(opt.tline2, 'ydata', tmp);
  end

  set(opt.hy,    'yaxislocation', 'right');
  if isfield(opt, 'parcellation')
    set(opt.stringz, 'string', sprintf('location = %s', opt.parcellationlabel{opt.parcellation(opt.vindx)}));
  else
    set(opt.stringz, 'string', sprintf('location = [%2.1f, %2.1f, %2.1f]', opt.pos(opt.vindx,:)));
  end
end

if opt.record
  tmp = get(opt.h, 'position');
  opt.frame = opt.frame + 1;
  opt.movie(opt.frame) = getframe(opt.h,[1 1 tmp(3:4)-1]);
end
setappdata(h, 'opt', opt);

previous_valx = valx;
previous_valy = valy;
previous_vindx = opt.vindx;

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

opt = getappdata(h, 'opt');
if strcmp(get(opt.playbutton, 'string'), 'pause')
  stop(opt.t);
  set(opt.playbutton, 'string', 'play');
else
  start(opt.t);
  set(opt.playbutton, 'string', 'pause');
end
setappdata(h, 'opt', opt);

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

opt = getappdata(h, 'opt');
opt.cleanup = 1;
setappdata(h, 'opt', opt);
uiresume(h);

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

opt = getappdata(h, 'opt');
if strcmp(get(opt.recordbutton, 'string'), 'stop')
  opt.record = 0;
  set(opt.recordbutton, 'string', 'record');
else
  opt.record = 1;
  set(opt.recordbutton, 'string', 'stop');
end
setappdata(h, 'opt', opt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_timer(obj, info, h)

opt   = getappdata(h, 'opt');
delta = opt.speed/size(opt.dat,2);
val = get(opt.sliderx, 'value');
val = val + delta;
if val>1
  val = val-1;
end
set(opt.sliderx, 'value', val);
setappdata(h, 'opt', opt);
cb_slider(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_alim(h, eventdata)
if ~ishandle(h)
  return
end
opt = guidata(h);
switch get(h, 'String')
  case '+'
    alim(alim*sqrt(2));
  case '-'
    alim(alim/sqrt(2));
end % switch
guidata(h, opt);
uiresume;

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

h   = getparent(h);
opt = getappdata(h, 'opt');
if strcmp(get(get(h, 'currentaxes'), 'tag'), 'timecourse')
  % get the current point
  %pos = get(opt.hy, 'currentpoint');
  %set(opt.sliderx, 'value', nearest(opt.xparam, pos(1)));
elseif strcmp(get(get(h, 'currentaxes'), 'tag'), 'mesh')
  % get the current point, which is defined as the intersection through the
  % axis-box (in 3D)
  pos       = get(opt.hx, 'currentpoint');

  % get the intersection with the mesh
  [ipos, d] = intersect_line(opt.pos, opt.tri, pos(1,:), pos(2,:));
  [md, ix]  = min(abs(d));

  if ~isempty(ix)
    dpos      = opt.pos - ipos(ix*ones(size(opt.pos,1),1),:);
    opt.vindx = nearest(sum(dpos.^2,2),0);
  end
end
setappdata(h, 'opt', opt);
cb_slider(h);
uiresume;

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

if (isempty(eventdata) && ft_platform_supports('matlabversion', -Inf, '2014a')) || isa(eventdata, 'matlab.ui.eventdata.ActionData')
  % determine the key that corresponds to the uicontrol element that was activated
  key = get(h, '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

h = getparent(h);
opt = getappdata(h, 'opt');

switch key
  case 'leftarrow' % change colorlim
    cval = opt.cfg.funcolorlim(1);
    if cval<0
      cval = cval.*1.1;
    else
      cval = cval./1.1;
    end
    opt.cfg.funcolorlim(1) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'rightarrow' % change colorlim
    cval = opt.cfg.funcolorlim(2);
    if cval<0
      cval = cval.*1.1;
    else
      cval = cval./1.1;
    end
    opt.cfg.funcolorlim(2) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'shift+leftarrow'
    cval = opt.cfg.funcolorlim(1);
    if cval<0
      cval = cval./1.1;
    else
      cval = cval.*1.1;
    end
    opt.cfg.funcolorlim(1) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'shift+rightarrow'
    cval = opt.cfg.funcolorlim(2);
    if cval<0
      cval = cval./1.1;
    else
      cval = cval.*1.1;
    end
    opt.cfg.funcolorlim(2) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'uparrow' % enhance threshold
    thrval = opt.threshold(2);
    if thrval==0 && opt.cfg.funcolorlim(2)>0
      opt.threshold(2) = opt.cfg.funcolorlim(2).*0.1;
    elseif thrval==0
      opt.threshold(2) = max(opt.dat(:)).*0.1;
    else
      opt.threshold(2) = thrval.*1.1;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline1, 'YData', [1 1].*opt.threshold(2));
  case 'shift+uparrow' % enhance threshold
    thrval = opt.threshold(2);
    if thrval>0
      opt.threshold(2) = thrval.*0.9;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline1, 'YData', [1 1].*opt.threshold(2));
  case 'downarrow' % lower threshold
    thrval = opt.threshold(1);
    if thrval==0 && opt.cfg.funcolorlim(1)<0
      opt.threshold(1) = opt.cfg.funcolorlim(1).*0.1;
    elseif thrval==0
      opt.threshold(1) = min(opt.dat(:)).*0.1;
    else
      opt.threshold(1) = thrval.*1.1;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline2, 'YData', [1 1].*opt.threshold(1));
  case 'shift+downarrow' % lower threshold
    thrval = opt.threshold(1);
    if thrval<0
      opt.threshold(1) = thrval.*0.9;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline2, 'YData', [1 1].*opt.threshold(1));
  case 'ctrl+uparrow' % change speed
    opt.speed = opt.speed*sqrt(2);
    setappdata(h, 'opt', opt);
  case 'ctrl+downarrow'
    opt.speed = opt.speed/sqrt(2);
    opt.speed = max(opt.speed, 1); % should not be smaller than 1
    setappdata(h, 'opt', opt);
  case 'ctrl+rightarrow' % change channel
  case 'C' % update camera position1.373e-14
    camlight(opt.cam(1), 'left');
    camlight(opt.cam(2), 'right');
  case 'p'
    cb_playbutton(h);
  case 'q'
    cb_quitbutton(h);
  case 'r'
    cb_recordbutton(h);
  case 's'
    % select the speed
    response = inputdlg('speed', 'specify', 1, {num2str(opt.speed)});
    if ~isempty(response)
      opt.speed = str2double(response);
      setappdata(h, 'opt', opt);
    end
  case 't'
    % select the threshold
    response = inputdlg('threshold', 'specify', 1, {num2str(opt.threshold)});
    if ~isempty(response)
      tok = tokenize(response{1}, ' ');
      tok(cellfun(@isempty,tok)) = [];
      for k = 1:numel(tok)
        opt.threshold(1,k) = str2double(tok{k});
      end
      setappdata(h, 'opt', opt);
    end
    set(opt.hline1, 'YData', [1 1].*opt.threshold(2));
    set(opt.hline2, 'YData', [1 1].*opt.threshold(1));
  case 'z'
    % select the colorlim
    response = inputdlg('colorlim', 'specify', 1, {num2str(opt.cfg.funcolorlim)});
    if ~isempty(response)
      tok = tokenize(response{1}, ' ');
      tok(cellfun(@isempty,tok)) = [];
      for k = 1:2
        opt.cfg.funcolorlim(1,k) = str2double(tok{k});
      end
      setappdata(h, 'opt', opt);
      cb_slider(h);
    end
  case 'control+control'
    % do nothing
  case 'shift+shift'
    % do nothing
  case 'alt+alt'
    % do nothing
  otherwise
    setappdata(h, 'opt', opt);
    %cb_help(h);
end
cb_slider(h);
uiresume(h);

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


text:
import client.EEGLABChat;
import client.VisualToolbar;
import java.awt.*;
import javax.swing.*;

tb = VisualToolbar();
F = gcf;
tb.setPreferredSize(Dimension(0, 75));

javacomponent(tb,'South',F);

refresh(F);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function[]=makefigs_periodize
%MAKEFIGS_PERIODIZE  Makes a sample figure for PERIODIZE.

load qgsnapshot
[xp,yp,fp]=periodize(100,200,qgsnapshot.x,qgsnapshot.y,qgsnapshot.psi);

figure
pcolor(xp,yp,fp),set(gca,'dataaspectratio',[1 1 1]),axis tight, shading interp
vlines(xp([100 end-100]),'w')
hlines(yp([200 end-200]),'w')
xtick([-5:1:5]*1000),ytick([-5:1:5]*1000)
title('Periodization of 1024x1024 QG turbulence with N=100 and M=200.')

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab