text:
function [sCodebook, settings] = grlvq(sCodebook, sData, varargin)

%GRLVQ trains a codebook using the Generalized Relevance LVQ (GRLVQ)
%algorithm.
%
% sM = grlvq(sM, sD, [argID, value, ...])
%
%  sM = grlvq(sM, sD)
%  sM = grlvq(sM, sD, 'learningRatePrototypes', [0.01 0.001], 'relevanceStart', 10);
%
% Required Input Arguments:
%   sM                     (struct) map struct containing the labels as the
%                                   first column of .labels
%   sD                     (struct) data struct containing the labels as
%                                   the first column of .labels
%
% Optional Input Arguments:
%   PrototypesPerClass     (vector) number of prototypes per class used;
%                                   a number or a vector with a number
%                                   for each class (default=1)
%   initialRelevances      (vector) initial Relevances
%   regularization         (scalar) regularization for relevances
%   testSet                (matrix) a test set to compute the test error;
%                                   the last column should be the labels
%   comparable             (scalar) a flag which resets the random
%                                   generator to produce comparable results
%                                   if set to 1
%   optimization           (string) choice for the optimization technique:
%                                   sgd or fminlbfgs (default=fminlbfgs)
%  Parameter for the stochastic gradient descent sgd:
%   nb_epochs              (scalar) the number of epochs (default=100)
%   learningRatePrototypes (vector) learning rate for the prototypes; could
%                                   be the start and end value or a vector
%                                   of length nb_epochs
%   learningRateRelevances (vector) learning rate for the relevances; could
%                                   be the start and end value or a vector
%                                   of length nb_epochs
%   relevanceStart         (scalar) epoch to start the relevance learning
%                                   (default=1)
%  Parameter for the build-in function fminlbfgs:
%   threshstop             (scalar) the training error for early stopping
%                                   (default=0)
%   useEarlyStopping       (scalar) use early stopping based on threshstop
%                                   (default=1)
%   Display                (string) the optimization output 'iter' or 'off'
%                                   (default='off')
%   GradObj                (string) turn the usage of gradient information
%                                   on or off (default='on')
%   HessUpdate             (string) the update can be 'lbfgs', 'bfgs' or
%                                   'steepdesc' (default='lbfgs')
%   TolFun                 (scalar) the tolerance (default=1e-6)
%   MaxIter                (scalar) the maximal number of iterations
%                                   (default=2500)
%   MaxFunEvals                     the maximal number of function
%                                   evaluations (default=1000000)
%   TolX                   (scalar) tolerance on the minimum
%                                   (default=1e-10)
%   DiffMinChange          (scalar) minimal change (default=1e-10)
%
% Output Arguments:
%   sM                     (struct) map struct containing the optimized
%                                   prototypes and the relevances lambda
%   settings               (struct) information on the settings of the
%                                   algorithm
%
% NOTE: does not take the vector mask into account but trains a relevance
%       for each dimension.
%
% For more help, try 'type grlvq', or check out the online documentation.
% See also GMLVQ, GRLVQ_CORE, LVQ3, LVQ1.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% grlvq
%
% PURPOSE
%
% Trains a codebook and a global metric with the Generalized Relevance LVQ
% (GRLVQ) algorithm.
%
% SYNTAX
%
%   sM = grlvq(sM, sD)
%   sM = grlvq(sM, sD, 'learningRatePrototypes', [0.01 0.001], 'relevanceStart', 10);
%   [sM, settings] = grlvq(sM, sD, argID, value, ...);
%
% DESCRIPTION
%
% GRLVQ constitutes an enhancement of the LVQ algorithm by minimizing an
% explicit error function and optimizing the metric of the space where the
% computation is performed. Distances are calculated by
%   d(x,y) = (x - y) * diag(lambda) * (x - y)'.
% The vector lambda is optimized by the algorithm.
%
% Two optimization techniques are implemented:
%     - stochastic gradient descent (sgd)
%     - limited memory Quasi Newton Broyden-Fletcher-Goldfarb-Shanno
%       (L-BFGS)
%
% Classification is performed similarly to standard LVQ, with the only
% difference that for the computation of distances the above specified
% formula is applied.
%
% This function provides an interface in the style of the som-toolbox and
% is basically a wrapper for the method grlvq_core. It deals with structs
% rather then directly with data matrices.
%
% REFERENCES
%
%   Barbara Hammer, Thomas Villmann: Generalized relevance learning vector
%   quantization. Neural Networks 15(8-9): 1059-1068 (2002).
%
%   P. Schneider, K. Bunte, B. Hammer and M. Biehl: Regularization in
%   Matrix Relevance Learning, IEEE Transactions on Neural Networks, vol.
%   21, nb. 5, pp. 831-840, 2010.
%
% REQUIRED INPUT ARGUMENTS
%   sM                     (struct) map struct containing the initialized
%                                   prototype vectors as rows of the
%                                   .codebooks matrix and the labels as the
%                                   first column of .labels
%   sD                     (struct) data struct containing the data vectors
%                                   as rows of .data and the labels as the
%                                   first column of .labels
%
% OPTIONAL INPUT ARGUMENTS
%   PrototypesPerClass     (vector) number of prototypes per class used;
%                                   a number or a vector with a number
%                                   for each class (default=1)
%   initialRelevances      (vector) initial Relevances
%                                   (default: vector of ones)
%   regularization         (scalar) regularization for relevances
%                                   (detault=0)
%   comparable             (scalar) a flag which resets the random
%                                   generator to produce comparable results
%                                   if set to 1 (default=0)
%   optimization           (string) choice for the optimization technique:
%                                   sgd (stochastic gradient descent) or
%                                   fminlbfgs (Limited memory Quasi Newton
%                                   Broyden-Fletcher-Goldfarb-Shanno)
%                                   (default=fminlbfgs)
%  Parameter for the stochastic gradient descent sgd:
%   nb_epochs              (scalar) the number of epochs (default=100)
%   learningRatePrototypes (vector) learning rate for the prototypes; could
%                                   be the start and end value or a vector
%                                   of length nb_epochs
%   learningRateRelevances (vector) learning rate for the relevances; could
%                                   be the start and end value or a vector
%                                   of length nb_epochs
%   relevanceStart         (scalar) epoch to start the relevance learning
%                                   (default=1)
%  Parameter for the build-in function fminlbfgs:
%   threshstop             (scalar) the training error for early stopping
%                                   (default=0)
%   useEarlyStopping       (scalar) use early stopping based on threshstop
%                                   (default=1)
%   Display                (string) the optimization output 'iter' or 'off'
%                                   (default='off')
%   GradObj                (string) turn the usage of gradient information
%                                   on or off (default='on')
%   HessUpdate             (string) the update can be 'lbfgs', 'bfgs' or
%                                   'steepdesc' (default=lbfgs)
%   TolFun                 (scalar) the tolerance (default=1e-6)
%   MaxIter                (scalar) the maximal number of iterations
%                                   (default=2500)
%   MaxFunEvals                     the maximal number of function
%                                   evaluations (default=1000000)
%   TolX                   (scalar) tolerance on the minimum
%                                   (default=1e-10)
%   DiffMinChange          (scalar) minimal change (default=1e-10)
%
% OUTPUT ARGUMENTS
%   sM                     (struct) map struct containing the optimized
%                                   prototypes as well as the relevances
%                                   for each dimensionlambda
%   settings               (struct) information on the settings of the
%                                   algorithm
%
% EXAMPLES
%
%   lab   = unique(sD.labels(:,1));           % different classes
%   nProt = length(lab)*5;                    % 5 prototypes for each
%   sM = som_randinit(sD,'msize',[nProt 1]);  % initial prototypes
%   sM.labels = [lab;lab;lab;lab;lab];        % and their classes
%   sM = grlvq(sM,sD);                        % use GRLVQ to adjust the
%                                             % prototypes and the metric
%
% SEE ALSO
%
%   gmlvq         Use the GMLVQ algorithm for training.
%   grlvq_core    Access the GRLVQ functionality without using structs.
%   lvq3          Use LVQ3 algorithm for training.
%   lvq1          Use LVQ1 algorithm for training.

% Contributed to SOM Toolbox vs2, December 3rd, 2012 by Alexander Schulz
% Copyright (c) Alexander Schulz
% http://www.cis.hut.fi/projects/somtoolbox/

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

cods      = sCodebook.codebook;
codLabels = sCodebook.labels;
data      = sData.data;
labels    = sData.labels;

% convert labels from cell to a vector
labelsV    = som_label2num(labels);
codLabelsV = som_label2num(codLabels);

% calculate the number of prototypes per class
n_prot_per_class = zeros(1, max(codLabelsV));
for i=1:max(codLabelsV)
    n_prot_per_class(i) = sum(codLabelsV == i);
end

% call the main function
[model, settings]  = grlvq_core(data, labelsV, 'initialPrototypes', ...
    [cods, codLabelsV], 'PrototypesPerClass', n_prot_per_class, varargin{:});

% write the results in the output struct
sCodebook.codebook = model.w;
sCodebook.lambda   = model.lambda;

if strcmp(settings.optimization, 'sgd')
    trainlen   = settings.nb_epochs;
    alpha_ini  = settings.learningRatePrototypes(1);
    alpha_type = 'power';
else
    trainlen   = NaN;
    alpha_ini  = NaN;
    alpha_type = '';
end

sTrain = som_set('som_train','algorithm','GRLVQ',...
		 'data_name',sData.name,...
		 'neigh','',...
		 'mask',ones(size(model.w,2),1),...
		 'radius_ini',NaN,...
		 'radius_fin',NaN,...
		 'alpha_ini',alpha_ini,...
		 'alpha_type',alpha_type,...
		 'trainlen',trainlen,...
		 'time',datestr(now,0));
sCodebook.trainhist(end+1) = sTrain;

author: ilarinieminen

repo: SOM-Toolbox


text:
function [ y, m ] = month_carry_alexandrian ( y, m )

%*****************************************************************************80
%
%% MONTH_CARRY_ALEXANDRIAN carries a year of months on the Alexandrian calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, the YM date.
%
  months = year_length_months_alexandrian ( y );

  while ( 1 )

    if ( m <= months )
      break
    end

    m = m - months;
    y = y + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [TorF, vstr, rdate] = have_feature_optimoptions()
%HAVE_FEATURE_OPTIMOPTIONS  Detect availability/version info for OPTIMOPTIONS
%
%   Feature detection function implementing 'optimoptions' tag for HAVE_FEATURE
%   to detect support for OPTIMOPTIONS, option setting function for MATLAB
%   Optimization Toolbox 6.3 and later.
%
%   See also HAVE_FEATURE, HAVE_FEATURE_OPTIM, OPTIMOPTIONS.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

TorF = 0;
vstr = '';
rdate = '';
if have_feature('matlab')
    v = have_feature('optim', 'all');
    if v.av && v.vnum >= 6.003      %% Opt Tbx 6.3+ (R2013a+, MATLAB 8.1+)
        TorF = 1;
        vstr = v.vstr;
        rdate = v.date;
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function b = c8mat_fss ( n, a, nb, b )

%*****************************************************************************80
%
%% C8MAT_FSS factors and solves a system with multiple right hand sides.
%
%  Discussion:
%
%    This routine uses partial pivoting, but no pivot vector is required.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, complex A(N,N), the coefficient matrix of the linear system.
%
%    Input, integer NB, the number of right hand sides.
%
%    Input, complex B(N,NB), the right hand side of the linear system.
%
%    Output, complex B(N,NB), the solution of the linear system.
%
  info = 0;

  for jcol = 1 : n
%
%  Find the maximum element in column I.
%
    piv = abs ( a(jcol,jcol) );
    ipiv = jcol;
    for i = jcol+1 : n
      if ( piv < abs ( a(i,jcol) ) )
        piv = abs ( a(i,jcol) );
        ipiv = i;
      end
    end

    if ( piv == 0.0 )
      info = jcol;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'C8MAT_FSS - Fatal error!\n' );
      fprintf ( 1, '  Zero pivot on step %d\n', info );
      return
    end
%
%  Switch rows JCOL and IPIV, and B.
%
    if ( jcol ~= ipiv )

      for j = 1 : n
        t         = a(jcol,j);
        a(jcol,j) = a(ipiv,j);
        a(ipiv,j) = t;
      end

      t(1,1:nb)    = b(jcol,1:nb);
      b(jcol,1:nb) = b(ipiv,1:nb);
      b(ipiv,1:nb) = t(1,1:nb);

    end
%
%  Scale the pivot row.
%
    temp = a(jcol,jcol);
    a(jcol,jcol) = 1.0;
    a(jcol,jcol+1:n) = a(jcol,jcol+1:n) / temp;
    b(jcol,1:nb) = b(jcol,1:nb) / temp;
%
%  Use the pivot row to eliminate lower entries in that column.
%
    for i = jcol+1 : n
      if ( a(i,jcol) ~= 0.0 )
        temp = - a(i,jcol);
        a(i,jcol) = 0.0;
        a(i,jcol+1:n) = a(i,jcol+1:n) + temp * a(jcol,jcol+1:n);
        b(i,1:nb) = b(i,1:nb) + temp * b(jcol,1:nb);
      end
    end

  end
%
%  Back solve.
%
  for j = 1 : nb
    for jcol = n : -1 : 2
      b(1:jcol-1,j) = b(1:jcol-1,j) - a(1:jcol-1,jcol) * b(jcol,j);
    end
  end

  return
end


text:
function A = linop_explicit( op )
%LINOP_EXPLICIT Outputs the explicit matrix representation
%   of a (implicitly defined) linear operator.
%   Useful for checking correctness of code.
% A = LINOP_EXPLICIT( OP )
%   returns the matrix A such that A*X = OP(X)
%
% Note: may not play well with linear operators defined on complex numbers
%   (but should be OK if the output inclides complex numbers)

% Introduced June 2016

if nargin == 0,
    error( 'Not enough input arguments.' );
end
sz = op([],0);
if isnumeric(sz),
    sz = { [sz(2),1], [sz(1),1] };
end
% convert [n1;n2] to [n1,n2] if necessary:
for kk = 1:2
    sz{kk} = sz{kk}(:).';
end
% If inputs and outputs are not vectors, then we cannot represent
%   with matrix multiplication
if sz{1}(2) ~= 1 || sz{2}(2) ~= 1
    error('Cannot represent this operator as matrix since input/output is not a vector');
end
m   = sz{2}(1);
n   = sz{1}(1);
e   = zeros(n,1);
A   = zeros(m,n);
for i = 1:n
    e(i) = 1;
    A(:,i)  = op( e, 1 );
    e(i) = 0;
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function V = VBA_spm_vol_nifti(P,n)
% Get header information for a NIFTI-1 image
% FORMAT V = spm_vol_nifti(P,n)
% P   - filename or NIfTI object
% n   - volume id (a 1x2 array, e.g. [3,1])
%
% V   - a structure containing the image volume information
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_vol_nifti.m 6079 2014-06-30 18:25:37Z spm $

if nargin<2,  n = [1 1];      end
if ischar(n), n = str2num(n); end

if ischar(P)
    N = nifti(P);
else
    N = P;
end
n  = [n 1 1];
n  = n(1:2);
dm = [N.dat.dim 1 1 1 1];
if any(n>dm(4:5)), V = struct([]); return; end

dt = struct(N.dat);
dt = double([dt.dtype dt.be]);

if isfield(N.extras,'mat') && ...
        size(N.extras.mat,3) >= n(1) && ...
        sum(sum(N.extras.mat(:,:,n(1)))) ~= 0
    mat = N.extras.mat(:,:,n(1));
else
    mat = N.mat;
end

off = (n(1)-1+dm(4)*(n(2)-1))*ceil(spm_type(dt(1),'bits')*dm(1)*dm(2)/8)*dm(3) + ...
    N.dat.offset;

V   = struct('fname',   N.dat.fname,...
             'dim',     dm(1:3),...
             'dt',      dt,...
             'pinfo',   [N.dat.scl_slope N.dat.scl_inter off]',...
             'mat',     mat,...
             'n',       n,...
             'descrip', N.descrip,...
             'private', N);

author: MBB-team

repo: VBA-toolbox


text:
function plotConfMat(varargin)
%PLOTCONFMAT plots the confusion matrix with colorscale, absolute numbers
%   and precision normalized percentages
%
%   usage:
%   PLOTCONFMAT(confmat) plots the confmat with integers 1 to n as class labels
%   PLOTCONFMAT(confmat, labels) plots the confmat with the specified labels
%
%   Vahe Tshitoyan
%   20/08/2017
%
%   Arguments
%   confmat:            a square confusion matrix
%   labels (optional):  vector of class labels

% number of arguments
switch (nargin)
    case 0
       confmat = 1;
       labels = {'1'};
    case 1
       confmat = varargin{1};
       labels = 1:size(confmat, 1);
    otherwise
       confmat = varargin{1};
       labels = varargin{2};
end

confmat(isnan(confmat))=0; % in case there are NaN elements
numlabels = size(confmat, 1); % number of labels

% calculate the percentage accuracies
confpercent = 100*confmat./repmat(sum(confmat, 1),numlabels,1);

% plotting the colors
imagesc(confpercent);
title(sprintf('Accuracy: %.2f%%', 100*trace(confmat)/sum(confmat(:))));
%ylabel('Predicted Class'); xlabel('Target Class');
ylabel('Predicted Class'); xlabel('Target Class');

% set the colormap
colormap(flipud(gray));

% Create strings from the matrix values and remove spaces
textStrings = num2str([confpercent(:), confmat(:)], '%.1f%%\n%d\n');
textStrings = strtrim(cellstr(textStrings));

% Create x and y coordinates for the strings and plot them
[x,y] = meshgrid(1:numlabels);
hStrings = text(x(:),y(:),textStrings(:), ...
    'HorizontalAlignment','center');

% Get the middle value of the color range
midValue = mean(get(gca,'CLim'));

% Choose white or black for the text color of the strings so
% they can be easily seen over the background color
textColors = repmat(confpercent(:) > midValue,1,3);
set(hStrings,{'Color'},num2cell(textColors,2));

% Setting the axis labels
set(gca,'XTick',1:numlabels,...
    'XTickLabel',labels,...
    'YTick',1:numlabels,...
    'YTickLabel',labels,...
    'TickLength',[0 0]);
end

author: opencobra

repo: cobratoolbox


text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 9
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%% Figure 9.1a

% load EEG data
load sampleEEGdata.mat

% plot a few trials from one channel...

% specify the label of the channel to plot
which_channel_to_plot = 'fcz';
% and find the index (channel number) of that label
channel_index = strcmpi(which_channel_to_plot,{EEG.chanlocs.labels});

x_axis_limit = [-200 1000]; % in ms

num_trials2plot = 12;

figure
set(gcf,'Name',[ num2str(num_trials2plot) ' random trials from channel ' which_channel_to_plot ],'Number','off')
for i=1:num_trials2plot

    % figure out how many subplots we need
    subplot(ceil(num_trials2plot/ceil(sqrt(num_trials2plot))),ceil(sqrt(num_trials2plot)),i)

    % pick a random trial (using randsample, which is in the stats toolbox)
    random_trial_to_plot = randsample(EEG.trials,1);
    % if you don't have the stats toolbox, use the following two lines:
    %random_trial_to_plot = randperm(EEG.trials);
    %random_trial_to_plot = random_trial_to_plot(1);

    % plot trial and specify x-axis and title
    plot(EEG.times,squeeze(EEG.data(channel_index,:,random_trial_to_plot)));
    set(gca,'xlim',x_axis_limit,'ytick',[])
    title([ 'Trial ' num2str(random_trial_to_plot) ])
end

%% Figure 9.1b

figure
% plot all trials
plot(EEG.times,squeeze(EEG.data(channel_index,:,:)),'y')

hold on
% plot ERP (simply the average time-domain signal)
plot(EEG.times,squeeze(mean(EEG.data(channel_index,:,:),3)),'k','linew',2)
set(gca,'xlim',[-300 1000],'ylim',[-60 60],'ydir','reverse')

% now plot only the ERP
figure
plot(EEG.times,squeeze(mean(EEG.data(channel_index,:,:),3))) % Note the "3" as second input to "mean"; this takes the average of the 3rd dimension.

% plot lines indicating baseline activity and stim onset
hold on
plot(get(gca,'xlim'),[0 0],'k')
plot([0 0],get(gca,'ylim'),'k:')

% add axis labels and title
xlabel('Time (ms)')
ylabel('\muV') % note that matlab interprets "\mu" as the Greek character for micro
title([ 'ERP (average of ' num2str(EEG.trials) ' trials) from electrode ' EEG.chanlocs(channel_index).labels ])

% plot upside down, following ERP convention
set(gca,'ydir','reverse')
% axis ij % this does the same as the previous trial

% below is some advanced but flexible code to change the x-axis label
set(gca,'xlim',[-300 1000])
xticklabel=cellstr(get(gca,'xticklabel'));
xticklabel{str2double(xticklabel)==0}='stim';
set(gca,'xticklabel',xticklabel)

%% Figure 9.2

% pick a channel
chan2plot = 'p7';

% compute ERP
erp = double(squeeze(mean(EEG.data(strcmpi(chan2plot,{EEG.chanlocs.labels}),:,:),3)));

% low-pass filter data (requires signal processing toolbox)
% you'll learn about how filtering works, and what this code means, in chapter 14.
nyquist          = EEG.srate/2;
transition_width = 0.15; % percent

% first, filter from 0-40
filter_cutoff = 40; % Hz
ffrequencies  = [ 0 filter_cutoff filter_cutoff*(1+transition_width) nyquist ]/nyquist;
idealresponse = [ 1 1 0 0 ];
filterweights = firls(100,ffrequencies,idealresponse);
erp_0to40     = filtfilt(filterweights,1,double(erp));

% next, filter from 0-10
filter_cutoff = 10; % Hz
ffrequencies  = [ 0 filter_cutoff filter_cutoff*(1+transition_width) nyquist ]/nyquist;
idealresponse = [ 1 1 0 0 ];
filterweights = firls(100,ffrequencies,idealresponse);
erp_0to10     = filtfilt(filterweights,1,double(erp));

% finally, filter from 5-15
filter_low    =  5; % Hz
filter_high   = 15; % Hz
ffrequencies  = [ 0 filter_low*(1-transition_width) filter_low filter_high filter_high*(1+transition_width) nyquist ]/nyquist;
idealresponse = [ 0 0 1 1 0 0 ];
filterweights = firls(round(3*(EEG.srate/filter_low)),ffrequencies,idealresponse);
erp_5to15     = filtfilt(filterweights,1,double(erp));

% now plot all filtered ERPs
figure
plot(EEG.times,erp,'k')
hold on
plot(EEG.times,erp_0to40,'b','linew',2)
plot(EEG.times,erp_0to10,'r')
plot(EEG.times,erp_5to15,'m')

set(gca,'xlim',[-200 1200],'ydir','r')
xlabel('Time (ms)'), ylabel('Voltage (\muV)')
title([ 'ERP from electrode ' chan2plot ])
legend({'None';'0-40';'0-10';'5-15'})

%% Figure 9.3

figure

% Butterfly plot
subplot(211)
plot(EEG.times,squeeze(mean(EEG.data,3)))
set(gca,'xlim',[-200 1000],'ydir','reverse')
xlabel('Time (ms)'), ylabel('\muV')
title('ERP from all sensors')

% topographical variance plot
subplot(212)
plot(EEG.times,squeeze(var( mean(EEG.data,3) ))) % note: change var to std for global field power
set(gca,'xlim',[-200 1000])
xlabel('Time (ms)'), ylabel('var(\muV)')
title('Topographical variance')

%% figure 9.4

% topoplot with colored dots vs. interpolated surface
% topoplot is a function in the eeglab toolbox, which you can download for
% free from the internet.

figure
subplot(121)
% To get the following line to work, you need to modify the eeglab function 'topoplot'
% Replace the first line of code (the function definition) with:
% function [handle,Zi,grid,Xi,Yi,x,y,pltchans] = topoplot(Values,loc_file,varargin)
[~,~,~,~,~,xi,yi,pltchans] = topoplot(zeros(64,1),EEG.chanlocs,'electrodes','off','plotrad',.53);

hold on
c = -squeeze(mean(EEG.data(:,300,:),3));
c = c-min(c); % these lines scale the variable c
c = c./max(c);
c = repmat(c,1,3);
for i=1:length(xi)
    hold on
    plot3(yi(i),xi(i),.1,'o','markerfacecolor',c(i,:),'markersize',5,'markeredge','none');
    % here, plot3 is used to plot the electrode on top of the topomap
    % (compare with the plot function) (a trick I learned from topoplot.m)
end

subplot(122)
c = -squeeze(mean(EEG.data(:,300,:),3));
c = c-min(c);
c = c./max(c);
c = repmat(c,1,3);
topoplot(double(c(:,1)),EEG.chanlocs,'plotrad',.53,'electrodes','off','numcontour',0);

colormap gray

%% introduction to topographical plotting

timepoint2plot = 100; % in ms
trial2plot = randsample(EEG.trials,1);

color_limit = 20; % more-or-less arbitrary, but this is a good value

% convert time point from ms to index
[junk,timepointidx] = min(abs(EEG.times-timepoint2plot));

% First step is to get X and Y coordinates of electrodes
% These must be converted to polar coordinates
th=pi/180*[EEG.chanlocs.theta];
[electrode_locs_X,electrode_locs_Y] = pol2cart(th,[EEG.chanlocs.radius]);

% interpolate to get nice surface
interpolation_level = 100; % you can try changing this number, but 100 is probably good
interpX = linspace(min(electrode_locs_X),max(electrode_locs_X),interpolation_level);
interpY = linspace(min(electrode_locs_Y),max(electrode_locs_Y),interpolation_level);

% meshgrid is a function that creates 2D grid locations based on 1D inputs
[gridX,gridY] = meshgrid(interpX,interpY);

% let's look at these matrices
figure
subplot(121)
imagesc(gridX)

subplot(122)
imagesc(gridY)

% now interpolate the data on a 2D grid
interpolated_EEG_data = griddata(electrode_locs_Y,electrode_locs_X,double(squeeze(EEG.data(:,timepointidx,trial2plot))),gridX,gridY);

figure
set(gcf,'number','off','name',[ 'Topographical data from trial ' num2str(trial2plot) ', time=' num2str(round(EEG.times(timepointidx))) ]);
subplot(121)
contourf(interpY,interpX,interpolated_EEG_data,100,'linecolor','none');
axis square
set(gca,'clim',[-color_limit color_limit],'xlim',[min(interpY) max(interpY)]*1.1,'ylim',[min(interpX) max(interpX)]*1.1)
title('Interpolated data in matlab')

subplot(122)
topoplot(double(squeeze(EEG.data(:,timepointidx,trial2plot))),EEG.chanlocs,'maplimits',[-color_limit color_limit]); % eeglab's topoplot function
title('eeglab ''topoplot'' function')

figure
set(gcf,'name','a landscape of cortical electrophysiological dynamics')
surf(interpY,interpX,interpolated_EEG_data);
xlabel('left-right of scalp'), ylabel('anterior-posterior of scalp'), zlabel('\muV')
shading interp, axis tight
set(gca,'clim',[-color_limit color_limit],'xlim',[min(interpY) max(interpY)]*1.1,'ylim',[min(interpX) max(interpX)]*1.1)
rotate3d on, view(0,90)

%% Figure 9.5

% find indices for requested plot times
times2plot = dsearchn(EEG.times',(-100:50:600)');

figure
for i=1:length(times2plot)
    subplot(3,5,i)

    % extract EEG data and replace FC4 data with noise
    eegdata2plot = double(squeeze(mean(EEG.data(:,times2plot(i),:),3)));
    eegdata2plot(strcmpi('fc4',{EEG.chanlocs.labels})) = randn*10;

    topoplot(eegdata2plot,EEG.chanlocs,'maplimits',[-8 8]);
    title([ num2str(round(EEG.times(times2plot(i)))) ' ms' ])
end

%% Figure 9.6

useRTs = true; % or false

% get RTs from each trial, to use for sorting trials. In this experiment,
% the RT was always the first event after the stimulus (the time=0 event).
% Normally, you should build in exceptions in case there was no response or
% another event occured between the stimulus and response. This was already
% done for the current dataset.

rts = zeros(size(EEG.epoch));
for ei=1:length(EEG.epoch)

    % first, find the index at which the time=0 event occurs
    time0event = find(cell2mat(EEG.epoch(ei).eventlatency)==0);

    % then get reaction time
    rts(ei) = EEG.epoch(ei).eventlatency{time0event+1};
end

% now find the sorting indices for RTs
if useRTs
    [dontneed,rts_idx] = sort(rts);
else
    [dontneed,rts_idx] = sort(squeeze(EEG.data(47,334,:)));
end

% now plot
figure
imagesc(EEG.times,1:EEG.trials,squeeze(EEG.data(47,:,rts_idx))');
set(gca,'clim',[-30 30],'xlim',[-200 1200],'ydir','n')

% also plot the RTs on each trial
if useRTs
    hold on
    plot(rts(rts_idx),1:EEG.trials,'k','linew',3)
end

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
function test_tutorial_beamformer(datadir)

% MEM 8gb
% WALLTIME 03:30:00
% DEPENDENCY ft_redefinetrial ft_freqanalysis ft_volumesegment ft_prepare_singleshell ft_sourceanalysis ft_prepare_leadfield ft_sourceinterpolate ft_sourceplot ft_volumenormalise

if nargin==0
  % this is where the data should be located
  datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer');
end

%load(fullfile(datadir, 'dataFIC.mat'));
load(fullfile(datadir, 'data_all.mat'));
load(fullfile(datadir, 'segmentedmri.mat'));
mri = ft_read_mri(fullfile(datadir, 'Subject01.mri'));

%% Preprocess time windows of interest

cfg = [];
cfg.toilim = [-0.5 0];
dataPre = ft_redefinetrial(cfg, data_all);

cfg.toilim = [0.8 1.3];
dataPost = ft_redefinetrial(cfg, data_all);

%% Cross-spectral density

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqPre = ft_freqanalysis(cfg, dataPre);

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqPost = ft_freqanalysis(cfg, dataPost);

%% Compute (or load) the forward model)

%try
  %if ~exist(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer/segmentedmri.mat'), 'file')
  cfg = [];
  cfg.write        = 'no';
  [segmentedmri] = ft_volumesegment(cfg, mri);
%catch
%  mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));
%  load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer/segmentedmri.mat'));
%end

%% Prepare head model
cfg = [];
cfg.method = 'singleshell';
vol = ft_prepare_headmodel(cfg, segmentedmri);

%% Prepare leadfield
cfg                 = [];
cfg.grad            = freqPost.grad;
cfg.headmodel       = vol;
cfg.reducerank      = 2;
cfg.channel         = {'MEG','-MLP31', '-MLO12'};
cfg.sourcemodel.resolution = 1;   % use a 3-D grid with a 1 cm resolution
cfg.sourcemodel.unit       = 'cm';
[grid] = ft_prepare_leadfield(cfg);

%% Source analysis without contrasting condition

cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.headmodel    = vol;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda       = 0;

sourcePost_nocon = ft_sourceanalysis(cfg, freqPost);

%save sourcePost_nocon.mat sourcePost_nocon;

% call ft_volumereslice so that figures appear correct side up
mri = ft_volumereslice([], mri);

% Plot the result
cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
sourcePostInt_nocon  = ft_sourceinterpolate(cfg, sourcePost_nocon , mri);

cfg              = [];
cfg.method       = 'slice';
cfg.funparameter = 'avg.pow';
figure
ft_sourceplot(cfg,sourcePostInt_nocon);

%% Compute and plot Neural Activity Index
sourceNAI = sourcePost_nocon;
sourceNAI.avg.pow = sourcePost_nocon.avg.pow ./ sourcePost_nocon.avg.noise;

cfg = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
sourceNAIInt = ft_sourceinterpolate(cfg, sourceNAI , mri);

cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [4.0 6.2];
cfg.opacitylim    = [4.0 6.2];
cfg.opacitymap    = 'rampup';
figure
ft_sourceplot(cfg, sourceNAIInt);

%% Exercise 4: lead field normalization
cfg                 = [];
cfg.grad            = freqPost.grad;
cfg.headmodel       = vol;
cfg.reducerank      = 2;
cfg.channel         = {'MEG','-MLP31', '-MLO12'};
cfg.sourcemodel.resolution = 1;   % use a 3-D grid with a 1 cm resolution
cfg.sourcemodel.unit       = 'cm';
cfg.normalize       = 'yes';
[gridn] = ft_prepare_leadfield(cfg);

cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = gridn;
cfg.headmodel    = vol;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda       = 0;
sourcePostn = ft_sourceanalysis(cfg, freqPost);

cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
sourcePostIntn  = ft_sourceinterpolate(cfg, sourcePostn , mri);
cfg              = [];
cfg.method       = 'slice';
cfg.funparameter = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [.6 1.2]*1e-26;
cfg.opacitylim    = [.6 1.2]*1e-26;
cfg.opacitymap    = 'rampup';
figure
ft_sourceplot(cfg,sourcePostIntn);

%% Source analysis with constrasting condition
dataAll = ft_appenddata([], dataPre, dataPost);

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqAll = ft_freqanalysis(cfg, dataAll);

cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.headmodel    = vol;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda       = '5%';
cfg.dics.keepfilter   = 'yes';
cfg.dics.realfilter   = 'yes';
sourceAll = ft_sourceanalysis(cfg, freqAll);

cfg.sourcemodel.filter = sourceAll.avg.filter;
sourcePre_con  = ft_sourceanalysis(cfg, freqPre );
sourcePost_con = ft_sourceanalysis(cfg, freqPost);

%save sourcePre_con sourcePre_con
%save sourcePost_con sourcePost_con

% %% Plot results
% cfg            = [];
% cfg.downsample = 2;
% cfg.parameter  = 'avg.pow';
% sourcePostInt  = ft_sourceinterpolate(cfg, sourcePost , mri);
%
% cfg              = [];
% cfg.method       = 'slice';
% cfg.funparameter = 'avg.pow';
% figure
% ft_sourceplot(cfg,sourcePostInt);
%
% %% Plot Condition contrast

sourceDiff = sourcePost_con;
sourceDiff.avg.pow = (sourcePost_con.avg.pow - sourcePre_con.avg.pow) ./ sourcePre_con.avg.pow;

cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
sourceDiffInt  = ft_sourceinterpolate(cfg, sourceDiff , mri);

cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0 1.2];
cfg.opacitylim    = [0 1.2];
cfg.opacitymap    = 'rampup';
figure
ft_sourceplot(cfg, sourceDiffInt);

%% Exercise 6: regularization
cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.headmodel    = vol;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda       = '0%';
cfg.dics.keepfilter   = 'yes';
cfg.dics.realfilter   = 'yes';
sourceAll = ft_sourceanalysis(cfg, freqAll);
cfg.sourcemodel.filter = sourceAll.avg.filter;
source0Pre  = ft_sourceanalysis(cfg, freqPre );
source0Post = ft_sourceanalysis(cfg, freqPost);
cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.headmodel    = vol;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda       = '10%';
cfg.dics.keepfilter   = 'yes';
cfg.dics.realfilter   = 'yes';
sourceAll = ft_sourceanalysis(cfg, freqAll);
cfg.sourcemodel.filter = sourceAll.avg.filter;
source10Pre  = ft_sourceanalysis(cfg, freqPre );
source10Post = ft_sourceanalysis(cfg, freqPost);

source0Diff = source0Post;
source0Diff.avg.pow = (source0Post.avg.pow - source0Pre.avg.pow) ./ source0Pre.avg.pow;
cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
source0DiffInt  = ft_sourceinterpolate(cfg, source0Diff , mri);
cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0 1.2];
cfg.opacitylim    = [0 1.2];
cfg.opacitymap    = 'rampup';
figure
ft_sourceplot(cfg, source0DiffInt);

source10Diff = source10Post;
source10Diff.avg.pow = (source10Post.avg.pow - source10Pre.avg.pow) ./ source10Pre.avg.pow;
cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
source10DiffInt  = ft_sourceinterpolate(cfg, source10Diff , mri);
cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0 1.2];
cfg.opacitylim    = [0 1.2];
cfg.opacitymap    = 'rampup';
figure
ft_sourceplot(cfg, source10DiffInt);

%% Orthogonal cut

cfg = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
figure
ft_sourceplot(cfg, sourceDiffInt);

%% Template MRI

cfg = [];
cfg.nonlinear     = 'no';
sourceDiffIntNorm = ft_volumenormalise(cfg, sourceDiffInt);

cfg = [];
cfg.method        = 'ortho';
%cfg.interactive   = 'yes';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
figure
ft_sourceplot(cfg, sourceDiffIntNorm);

%% Project to a surface
% cfg = [];
% cfg.nonlinear     = 'no';
% sourceDiffIntNorm = ft_volumenormalise(cfg, sourceDiffInt);

cfg = [];
cfg.method         = 'surface';
cfg.funparameter   = 'avg.pow';
cfg.maskparameter  = cfg.funparameter;
cfg.funcolorlim    = [0.0 1.2];
cfg.funcolormap    = 'jet';
cfg.opacitylim     = [0.0 1.2];
cfg.opacitymap     = 'rampup';
cfg.projmethod     = 'nearest';
cfg.surffile       = 'surface_white_both.mat';
cfg.surfdownsample = 10;
figure
ft_sourceplot(cfg, sourceDiffIntNorm);
view ([90 0])

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%DEMGMM1 Demonstrate EM for Gaussian mixtures.
%
%	Description
%	This script demonstrates the use of the EM algorithm to fit a mixture
%	of Gaussians to a set of data using maximum likelihood. A colour
%	coding scheme is used to illustrate the evaluation of the posterior
%	probabilities in the E-step of the EM algorithm.
%
%	See also
%	DEMGMM2, DEMGMM3, DEMGMM4, GMM, GMMEM, GMMPOST
%

%	Copyright (c) Ian T Nabney (1996-2001)

clc;
disp('This demonstration illustrates the use of the EM (expectation-')
disp('maximization) algorithm for fitting of a mixture of Gaussians to a')
disp('data set by maximum likelihood.')
disp(' ')
disp('The data set consists of 40 data points in a 2-dimensional')
disp('space, generated by sampling from a mixture of 2 Gaussian')
disp('distributions.')
disp(' ')
disp('Press any key to see a plot of the data.')
pause;

% Generate the data
randn('state', 0); rand('state', 0);
gmix = gmm(2, 2, 'spherical');
ndat1 = 20; ndat2 = 20; ndata = ndat1+ndat2;
gmix.centres =  [0.3 0.3; 0.7 0.7];
gmix.covars = [0.01 0.01];
x = gmmsamp(gmix, ndata);

h = figure;
hd = plot(x(:, 1), x(:, 2), '.g', 'markersize', 30);
hold on; axis([0 1 0 1]); axis square; set(gca, 'box', 'on');
ht = text(0.5, 1.05, 'Data', 'horizontalalignment', 'center');
disp(' ');
disp('Press any key to continue.')
pause; clc;

disp('We next create and initialize a mixture model consisting of a mixture')
disp('of 2 Gaussians having ''spherical'' covariance matrices, using the')
disp('function GMM. The Gaussian components can be displayed on the same')
disp('plot as the data by drawing a contour of constant probability density')
disp('for each component having radius equal to the corresponding standard')
disp('deviation. Component 1 is coloured red and component 2 is coloured')
disp('blue.')
disp(' ')
disp('Note that a particulary poor choice of initial parameters has been')
disp('made in order to illustrate more effectively the operation of the')
disp('EM algorithm.')
disp(' ')
disp('Press any key to see the initial configuration of the mixture model.')
pause;

% Set up mixture model
ncentres = 2; input_dim = 2;
mix = gmm(input_dim, ncentres, 'spherical');

% Initialise the mixture model
mix.centres = [0.2 0.8; 0.8, 0.2];
mix.covars = [0.01 0.01];

% Plot the initial model
ncirc = 30; theta = linspace(0, 2*pi, ncirc);
xs = cos(theta); ys = sin(theta);
xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
set(ht, 'string', 'Initial Configuration');
figure(h);
disp(' ')
disp('Press any key to continue');
pause; clc;

disp('Now we adapt the parameters of the mixture model iteratively using the')
disp('EM algorithm. Each cycle of the EM algorithm consists of an E-step')
disp('followed by an M-step.  We start with the E-step, which involves the')
disp('evaluation of the posterior probabilities (responsibilities) which the')
disp('two components have for each of the data points.')
disp(' ')
disp('Since we have labelled the two components using the colours red and')
disp('blue, a convenient way to indicate the value of a posterior')
disp('probability for a given data point is to colour the point using a')
disp('scale ranging from pure red (corresponding to a posterior probability')
disp('of 1.0 for the red component and 0.0 for the blue component) through')
disp('to pure blue.')
disp(' ')
disp('Press any key to see the result of applying the first E-step.')
pause;

% Initial E-step.
set(ht, 'string', 'E-step');
post = gmmpost(mix, x);
dcols = [post(:,1), zeros(ndata, 1), post(:,2)];
delete(hd);
for i = 1 : ndata
  hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ...
          'marker', '.', 'markersize', 30);
end
figure(h);

disp(' ');
disp('Press any key to continue')
pause; clc;

disp('Next we perform the corresponding M-step. This involves replacing the')
disp('centres of the component Gaussians by the corresponding weighted means')
disp('of the data. Thus the centre of the red component is replaced by the')
disp('mean of the data set, in which each data point is weighted according to')
disp('the amount of red ink (corresponding to the responsibility of')
disp('component 1 for explaining that data point). The variances and mixing')
disp('proportions of the two components are similarly re-estimated.')
disp(' ')
disp('Press any key to see the result of applying the first M-step.')
pause;

% M-step.
set(ht, 'string', 'M-step');
options = foptions;
options(14) = 1; % A single iteration
options(1) = -1; % Switch off all messages, including warning
mix = gmmem(mix, x, options);
delete(hc);
xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
figure(h);
disp(' ')
disp('Press any key to continue')
pause; clc;

disp('We can continue making alternate E and M steps until the changes in')
disp('the log likelihood at each cycle become sufficiently small.')
disp(' ')
disp('Press any key to see an animation of a further 9 EM cycles.')
pause;
figure(h);

% Loop over EM iterations.
numiters = 9;
for n = 1 : numiters

  set(ht, 'string', 'E-step');
  post = gmmpost(mix, x);
  dcols = [post(:,1), zeros(ndata, 1), post(:,2)];
  delete(hd);
  for i = 1 : ndata
    hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ...
                 'marker', '.', 'markersize', 30);
  end
  pause(1)

  set(ht, 'string', 'M-step');
  [mix, options] = gmmem(mix, x, options);
  fprintf(1, 'Cycle %4d  Error %11.6f\n', n, options(8));
  delete(hc);
  xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs;
  yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys;
  hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r');
  hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b');
  pause(1)

end

disp(' ');
disp('Press any key to end.')
pause; clc; close(h); clear all

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
clear all
close all
clc

i = 1;
switch i
    case 1
        load Dat_1.mat

    case 2
        load Dat_1.mat
end

N = 10000;
f = 1/N;
iterations = 12;
bits = 3;
levels = 2^bits;
m_max = 3.9;
delta = (2*m_max)/levels;
m = linspace(-m_max, m_max, levels+1);

for i = 1:iterations
    m
    for k = 1:levels

        sum=0;
        count=0;
        for j = 1:N
            if m(k) <= X(j) && X(j) < m(k+1)
                count = count + 1;
                sum = sum + X(j);
            end
        end

        if count == 0 && k <= levels/2
            v(k) = m(k);
        else
        if count == 0 && k > levels/2
            v(k) = m(k+1);
        else v(k) = sum/count;
        end
        end

    end

    for k = 2:levels-1
        m(k) = (v(k) + v(k+1))/2;
    end
   v
end

author: Sable

repo: mcbench-benchmarks


text:
classdef meanshift_Y_is_XP < dagnn.ElementWise
    properties (Transient)
        numInputs
        SIZE_
    end

    % [TODO]: current version only supports batchSize=1; need to extend to
    % multiple input images
    methods
        function outputs = forward(obj, inputs, params)
            obj.numInputs = numel(inputs);
            X = inputs{1};
            P = inputs{2};
            gpuMode = isa(X, 'gpuArray');
            % [height width channels batchsize]
            [h, w, ch, bs] = size(X);
            obj.SIZE_ = [h, w, ch, bs];
            if gpuMode
                Y = gpuArray(zeros(h, w, ch, bs, 'single'));
            else
                Y = zeros(h, w, ch, bs, 'single');
            end
            for i = 1:bs
                cur_X = X(:,:,:,i);
                cur_P = P(:,:,:,i);
                cur_X = reshape(cur_X, [h*w, ch]);
                cur_X = cur_X'*cur_P;
                cur_X = reshape(cur_X', [h, w, ch]);
                Y(:,:,:,i) = cur_X;
            end
            outputs{1} = Y;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            derInputs = cell(1, numel(inputs));
            h = obj.SIZE_(1);
            w = obj.SIZE_(2);
            ch = obj.SIZE_(3);
            bs = obj.SIZE_(4);

            dzdy = derOutputs{1};
            X = inputs{1};
            P = inputs{2};
            gpuMode = isa(X, 'gpuArray');
            if gpuMode
                dzdx_X = gpuArray(zeros(size(X,1), size(X,2), size(X,3), size(X,4), 'single'));
                dzdx_P = gpuArray(zeros(size(P,1), size(P,2), size(P,3), size(P,4), 'single'));
            else
                dzdx_X = zeros(size(X,1), size(X,2), size(X,3), size(X,4), 'single');
                dzdx_P = zeros(size(P,1), size(P,2), size(P,3), size(P,4), 'single');
            end
            for i = 1:size(X,4)
                cur_X = X(:,:,:,i);
                cur_P = P(:,:,:,i);

                cur_X = reshape(cur_X, [h*w, ch]);
                cur_X = cur_X';

                cur_dzdy = dzdy(:,:,:,i);
                cur_dzdy = reshape(cur_dzdy, [h*w, ch])';

                cur_dzdx_X = cur_dzdy*cur_P';
                cur_dzdx_P = cur_X'*cur_dzdy;

                cur_dzdx_X = reshape(cur_dzdx_X', [h, w, ch]);

                dzdx_X(:,:,:,i) = cur_dzdx_X;
                dzdx_P(:,:,:,i) = cur_dzdx_P;
            end

            derInputs{1} = dzdx_X;
            derInputs{2} = dzdx_P;
            derParams = {} ;
        end

        function obj = meanshift_Y_is_XP(varargin)
            obj.load(varargin) ;
        end
    end
end

author: aimerykong


text:
function Fitting_from_bb_vis( Image, DepthImage, bounding_box, PDM, patchExperts, clmParams, out_dir, varargin)
%FITTING Summary of this function goes here
%   Detailed explanation goes here

    % the bounding box format is [minX, minY, maxX, maxY];

    % the mean model shape
    M = PDM.M;

    num_points = numel(M) / 3;

    if(any(strcmp(varargin,'orientation')))
        orientation = varargin{find(strcmp(varargin, 'orientation'))+1};
        rot = Euler2Rot(orientation);
    else
        rot = eye(3);
        orientation = [0;0;0];
    end

    rot_m = rot * reshape(M, num_points, 3)';
    width_model = max(rot_m(1,:)) - min(rot_m(1,:));
    height_model = max(rot_m(2,:)) - min(rot_m(2,:));

    a = (((bounding_box(3) - bounding_box(1)) / width_model) + ((bounding_box(4) - bounding_box(2))/ height_model)) / 2;

    tx = (bounding_box(3) + bounding_box(1))/2;
    ty = (bounding_box(4) + bounding_box(2))/2;

    % correct it so that the bounding box is just around the minimum
    % and maximum point in the initialised face
    tx = tx - a*(min(rot_m(1,:)) + max(rot_m(1,:)))/2;
    ty = ty - a*(min(rot_m(2,:)) + max(rot_m(2,:)))/2;

    % visualisation of the initial state
    %hold off;imshow(Image);hold on;plot(a*rot_m(1,:)+tx, a*rot_m(2,:)+ty,'.r');hold on;rectangle('Position', [bounding_box(1), bounding_box(2), bounding_box(3)-bounding_box(1), bounding_box(4)-bounding_box(2)]);
    global_params = [a, 0, 0, 0, tx, ty]';
    global_params(2:4) = orientation;

    local_params = zeros(numel(PDM.E), 1);

    if(any(strcmp(varargin,'gparam')))
        global_params = varargin{find(strcmp(varargin, 'gparam'))+1};
    end

    if(any(strcmp(varargin,'lparam')))
        local_params = varargin{find(strcmp(varargin, 'lparam'))+1};
    end

    scale = clmParams.startScale;

    if(size(Image, 3) == 1)
        GrayImage = Image;
    else
        GrayImage = rgb2gray(Image);
    end

    [heightImg, widthImg] = size(GrayImage);

    % Some predefinitions for faster patch extraction
    [xi, yi] = meshgrid(0:widthImg-1,0:heightImg-1);
    xi = double(xi);
    yi = double(yi);

    GrayImageDb = double(GrayImage);

    % multi iteration refinement using NU-RLMS in each one
    i=1;

    current_patch_scaling = patchExperts(scale).trainingScale;
    visibilities = patchExperts(scale).visibilities;

    view = GetView(patchExperts(scale).centers, global_params(2:4));

    % The shape fitting is performed in the reference frame of the
    % patch training scale
    refGlobal = [current_patch_scaling, 0, 0, 0, 0, 0]';

    % the reference shape
    refShape = GetShapeOrtho(M, PDM.V, local_params, refGlobal);

    % shape around which the patch experts will be evaluated in the original image
    [shape2D] = GetShapeOrtho(M, PDM.V, local_params, global_params);
    shape2D_img = shape2D(:,1:2);

    % Create transform using a slightly modified version of Kabsch that
    % takes scaling into account as well, in essence we get a
    % similarity transform from current estimate to reference shape
    [A_img2ref, T_img2ref, ~, ~] = AlignShapesWithScale(shape2D_img(:,1:2),refShape(:,1:2));

    % Create a transform, from shape in image to reference shape
    T = maketform('affine', [A_img2ref;T_img2ref]);

    shape_2D_ref = tformfwd(T, shape2D_img);

    % transform the current shape to the reference one, so we can
    % interpolate
    shape2D_in_ref = (A_img2ref * shape2D_img')';

    sideSizeX = (clmParams.window_size(i,1) - 1)/2;
    sideSizeY = (clmParams.window_size(i,2) - 1)/2;

    patches = zeros(size(shape2D_in_ref,1), clmParams.window_size(i,1) * clmParams.window_size(i,2));

    Ainv = inv(A_img2ref);

    % extract patches on which patch experts will be evaluted
    for l=1:size(shape2D_in_ref,1)
        if(visibilities(view,l))

            xs = (shape2D_in_ref(l,1)-sideSizeX):(shape2D_in_ref(l,1)+sideSizeX);
            ys = (shape2D_in_ref(l,2)-sideSizeY):(shape2D_in_ref(l,2)+sideSizeY);

            [xs, ys] = meshgrid(xs, ys);

            pairs = [xs(:), ys(:)];

            actualLocs = (Ainv * pairs')';

            actualLocs(actualLocs(:,1) < 0,1) = 0;
            actualLocs(actualLocs(:,2) < 0,2) = 0;
            actualLocs(actualLocs(:,1) > widthImg - 1,1) = widthImg - 1;
            actualLocs(actualLocs(:,2) > heightImg - 1,2) = heightImg - 1;

            [t_patch] = interp2_mine(xi, yi, GrayImageDb, actualLocs(:,1), actualLocs(:,2), 'bilinear');
            t_patch = reshape(t_patch, size(xs));

            patches(l,:) = t_patch(:);

        end
    end

    % Calculate patch responses, either SVR or CCNF
    if(strcmp(patchExperts(scale).type, 'SVR'))
        responses = PatchResponseSVM_multi_modal( patches, patchExperts(scale).patch_experts(view,:), visibilities(view,:), patchExperts(scale).normalisationOptionsCol, clmParams, clmParams.window_size(i,:));

        for r=1:numel(responses)
            out_patch = reshape(patches(r,:)/255, size(xs));
            imwrite(out_patch, [out_dir, '/', num2str(r), '_a.png']);
           imwrite(responses{r}/max(responses{r}(:)), [out_dir, '/', num2str(r), '_svr.png']);
        end

    elseif(strcmp(patchExperts(scale).type, 'CCNF'))
        responses = PatchResponseCCNF( patches, patchExperts(scale).patch_experts(view,:), visibilities(view,:), patchExperts(scale), clmParams.window_size(i,:));
        for r=1:numel(responses)
           imwrite(responses{r}/max(responses{r}(:)), [out_dir, '/', num2str(r), '_lnf.png']);
        end
    elseif(strcmp(patchExperts(scale).type, 'DNN'))
        responses = PatchResponseDNN( patches, patchExperts(scale).patch_experts(view,:), visibilities(view,:), patchExperts(scale), clmParams.window_size(i,:));
        for r=1:numel(responses)
           imwrite(responses{r}/max(responses{r}(:)), [out_dir, '/', num2str(r), '_dnn.png']);
        end
    end

end

function [id] = GetView(centers, rotation)

    [~,id] = min(sum((centers * pi/180 - repmat(rotation', size(centers,1), 1)).^2,2));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
num_pes: 256
l1_size_cstr: 100
l2_size_cstr: 3000
noc_bw_cstr: 1000
offchip_bw_cstr: 50

author: maestro-project

repo: maestro


text:
function p = addInequality(p,row)
if ~anyCones(p.K)
    % Just append in the end
    p.F_struc = [p.F_struc;row];
else
    % Insert before conics
    if p.K.f==0 % a bit faster, append on top
        p.F_struc = [row;
                     p.F_struc];
    else
        p.F_struc = [p.F_struc(1:(p.K.f+p.K.l),:);
                     row;
                     p.F_struc(1+p.K.f+p.K.l:end,:)];
    end
end
p.K.l = p.K.l + size(row,1);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function grains = subSet(grains,ind)
%
%
% Input
%  grains - @grain2d
%  ind    - indices
%
% Output
%  grains - @grain2d
%

% restrict boundary
if islogical(ind)
  % the problem is grainId is with respect to grain.id
  % but ind is with respect to the order of the grains
  % therefore we have to enlarge ind
  indLarge = false(max(grains.boundary.grainId(:)),1);
  indLarge(grains.id) = ind;

  grId = grains.boundary.grainId;
  grId(grId>0) = indLarge(grId(grId>0));
  indBd = any(grId,2);
  %repeat for inner Boundary
  indinLarge = false(max(grains.innerBoundary.grainId(:)),1);
  indinLarge(grains.id) = ind;

  grinId = grains.innerBoundary.grainId;
  grinId(grinId>0) = indinLarge(grinId(grinId>0));
  indinnerBd = any(grinId,2);

else
  indBd = any(ismember(grains.boundary.grainId,grains.id(ind)),2);
  indinnerBd = any(ismember(grains.innerBoundary.grainId,grains.id(ind)),2);
end

grains = subSet@dynProp(grains,ind);

grains.poly = grains.poly(ind);
grains.inclusionId = grains.inclusionId(ind);
grains.id = grains.id(ind);
grains.phaseId = reshape(grains.phaseId(ind),[],1);
grains.grainSize = grains.grainSize(ind);

if ~islogical(ind)
  grains.prop.meanRotation = reshape(grains.prop.meanRotation, size(ind));
end

grains.boundary = subSet(grains.boundary,indBd);
grains.innerBoundary = subSet(grains.innerBoundary,indinnerBd);

% if we have only one grain - sort boundary segments
if length(grains) == 1
  FNew = [grains.poly{1}(1:end-1).',grains.poly{1}(2:end).'];

  % remove inclusion embeddings
  if grains.inclusionId > 0
    ie = sum(FNew == FNew(1),2)==1;
    ie([1,end-grains.inclusionId]) = false;
    FNew(ie,:) = [];
  end

  % sort minimum entry first
  FNew = sort(FNew,2);

  % sort such the order of F follows the boundary
  [~,ind1] = sortrows(grains.boundary.F);
  [~,ind2] = sortrows(FNew);
  inverseorder(ind2) = ind1;

  grains.boundary = grains.boundary(inverseorder);

end

author: mtex-toolbox

repo: mtex

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


text:
function voxels = carveall( voxels, cameras )
%CARVEALL  carve away voxels using all cameras
%
%   VOXELS = CARVEALL(VOXELS, CAMERAS) simple calls CARVE for each of the
%   cameras specified

%   Copyright 2005-2009 The MathWorks, Inc.
%   $Revision: 1.0 $    $Date: 2006/06/30 00:00:00 $

for ii=1:numel(cameras);
    voxels = carve(voxels,cameras(ii));
end


text:
%
% [Sin,SNin] = Sin_vsh_ctf_fast(r_sphere,R,EX,EY,EZ,Lin)
%
% Calculate the internal SSS basis Sin for a CTF system
% using vector spherical harmonics
%
function [Sin,SNin] = Sin_vsh_ctf_fast(r_sphere,R,EX,EY,EZ,Lin)

% Copyright (c) 2016, Elekta Oy
% ---------------------------------------
%
% Redistribution and use of the Software in source and binary forms, with or without
% modification, are permitted for non-commercial use.
%
% The Software is provided "as is" without warranties of any kind, either express or
% implied including, without limitation, warranties that the Software is free of defects,
% merchantable, fit for a particular purpose. Developer/user agrees to bear the entire risk
% in connection with its use and distribution of any and all parts of the Software under this license.
%

mu0 = 1.25664e-6; % Permeability of vacuum
%
% For numerical surface integration:
%
%baseline = 50e-3;
dx = 4.5e-3;
dy = 4.5e-3;
dz1 = 0;
dz2 = 50e-3;
D = [dx dy dz1; dx -dy dz1; -dx dy dz1; -dx -dy dz1; dx dy dz2; dx -dy dz2; -dx dy dz2; -dx -dy dz2]';
for j = 1:8
   if j <= 4
      weights(j) = 1/(4*1);
   else
      weights(j) = -1/(4*1);
   end
end
weights = weights';

for ch = 1:size(R,2)
   count = 1;
   R(:,ch) = R(:,ch) - r_sphere;
   Sin(ch,:) = -mu0*vsh_response(R(:,ch),EX(:,ch),EY(:,ch),EZ(:,ch),D,weights,Lin);
end
for j = 1:size(Sin,2)
   SNin(:,j) = Sin(:,j)/norm(Sin(:,j));
end

function Sin_elements = vsh_response(r,ex,ey,ez,D,weights,Lin)

for j = 1:length(weights)
    r_this = r + D(1,j)*ex + D(2,j)*ey;
    rn(j) = norm(r_this);
    theta(j) = acos(r_this(3)/rn(j));
    phi(j) = atan2(r_this(2),r_this(1));
    sint(j) = sin(theta(j));
    sinp(j) = sin(phi(j));
    cost(j) = cos(theta(j));
    cosp(j) = cos(phi(j));
    for l = 1:Lin
       p0{l}(:,j) = legendre(l,cos(theta(j)));
       rnv(j,l) = rn(j)^(l+2);
   end
end
Sin_elements = [];
for l = 1:Lin
  for m = -l:l
    for j = 1:length(weights)
      %vs = vsh_modified_in_fast(theta(j),phi(j),l,m,p0{l}(:,j))'/rn(j)^(l+2);
      vs = vsh_modified_in_fast(theta(j),phi(j),l,m,p0{l}(:,j))'/rnv(j,l);
      V(1,j) = vs(1)*sint(j)*cosp(j) + vs(2)*cost(j)*cosp(j) - vs(3)*sinp(j);
      V(2,j) = vs(1)*sint(j)*sinp(j) + vs(2)*cost(j)*sinp(j) + vs(3)*cosp(j);
      V(3,j) = vs(1)*cost(j) - vs(2)*sint(j);
    end
    Sin_elements = [Sin_elements dot(V*weights,ez)]; % Cartesian coordinates
                                                     %Sin_element = Sin_element/sqrt((l+1)*(2*l+1));  % Back to orthonormal presentation
  end
end
%Sin_element = dot(V*weights,ez); % Cartesian coordinates
%Sin_element = Sin_element/sqrt((l+1)*(2*l+1));  % Back to orthonormal presentation

author: spm


text:
% plot of sampling location on warped grid

rep = 'results/';

if ~exist(rep)
    mkdir(rep);
end

save_warped = 0;
save_dich = 1;

s = 50;
c = 0.45;
n = 16;
t = linspace(0,1,n);
[Y,X] = meshgrid(t,t);

x = X(:); y = Y(:);

xw = x;
yw = y + c * x.^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot of the warped locations
clf;
hold on;
scatter(x,y, s, 'filled');
plot(t, 1-c*t.^2);
hold off;

clf;
scatter(x,y, s, 'filled');
axis equal;
axis off;
if save_warped
saveas(gcf, [rep 'samples.eps'], 'epsc');
end

clf;
scatter(xw,yw, 25, 'filled');
axis equal; axis off;
if save_warped
saveas(gcf, [rep 'samples_warped.eps'], 'epsc');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot of the grouping process
posw = [xw, yw]';
pos = [x, y]';
options.part_type = '1axis';
for i=0:5
    options.depthmax = i;
    [part,B,G] = dichotomic_grouping([yw, xw]',options);
    options.point_size = 80;
    clf;
    plot_dichotomic_partition(posw, part, options);
    axis([0 1.5 0 1.5]); axis equal; axis off;
    if save_dich
        saveas(gcf, [rep '_dich_' num2str(i) '.eps'], 'epsc');
    end
    clf;
    plot_dichotomic_partition(pos, part, options);
    axis([0 1.5 0 1.5]); axis equal; axis off;
    if save_dich
        saveas(gcf, [rep '_dich_' num2str(i) 'w.eps'], 'epsc');
    end
end

author: gpeyre

repo: matlab-toolboxes


text:
function result = query(BIDS,query,varargin)
% Query a directory structure formated according to the BIDS standard
% FORMAT result = bids.query(BIDS,query,...)
% BIDS   - BIDS directory name or BIDS structure (from bids.layout)
% query  - type of query: {'data', 'metadata', 'sessions', 'subjects',
%          'runs', 'tasks', 'runs', 'types', 'modalities'}
% result - outcome of query
%__________________________________________________________________________
%
% BIDS (Brain Imaging Data Structure): https://bids.neuroimaging.io/
%   The brain imaging data structure, a format for organizing and
%   describing outputs of neuroimaging experiments.
%   K. J. Gorgolewski et al, Scientific Data, 2016.
%__________________________________________________________________________

% Copyright (C) 2016-2018, Guillaume Flandin, Wellcome Centre for Human Neuroimaging
% Copyright (C) 2018--, BIDS-MATLAB developers

if nargin < 2
    error('Not enough input arguments.');
end

BIDS = bids.layout(BIDS);

opts = parse_query(varargin);

switch query
%   case 'subjects'
%       result = regexprep(unique({BIDS.subjects.name}),'^[a-zA-Z0-9]+-','');
    case 'modalities'
        hasmod = arrayfun(@(y) structfun(@(x) isstruct(x) & ~isempty(x),y),...
            BIDS.subjects,'UniformOutput',false);
        hasmod = any([hasmod{:}],2);
        mods   = fieldnames(BIDS.subjects)';
        result = mods(hasmod);
    case {'sessions','subjects', 'tasks', 'runs', 'types', 'data', 'metadata'}
        %-Initialise output variable
        result = {};
        %-Filter according to subjects
        if any(ismember(opts(:,1),'sub'))
            subs = opts{ismember(opts(:,1),'sub'),2};
            opts(ismember(opts(:,1),'sub'),:) = [];
        else
            subs = unique({BIDS.subjects.name});
            subs = regexprep(subs,'^[a-zA-Z0-9]+-','');
        end
        %-Filter according to modality
        if any(ismember(opts(:,1),'modality'))
            mods = opts{ismember(opts(:,1),'modality'),2};
            opts(ismember(opts(:,1),'modality'),:) = [];
        else
            mods = bids.query(BIDS,'modalities');
        end
        %-Get optional target option for metadata query
        if strcmp(query,'metadata') && any(ismember(opts(:,1),'target'))
            target = opts{ismember(opts(:,1),'target'),2};
            opts(ismember(opts(:,1),'target'),:) = [];
            if iscellstr(target)
                target = substruct('.',target{1});
            end
        else
            target = [];
        end
        %-Perform query
        for i=1:numel(BIDS.subjects)
            if ~ismember(BIDS.subjects(i).name(5:end),subs), continue; end
            for j=1:numel(mods)
                d = BIDS.subjects(i).(mods{j});
                for k=1:numel(d)
                    sts = true;
                    for l=1:size(opts,1)
                        if ~isfield(d(k),opts{l,1}) || ~ismember(d(k).(opts{l,1}),opts{l,2})
                            sts = false;
                        end
                    end
                    switch query
                        case 'subjects'
                            if sts
                                result{end+1} = BIDS.subjects(i).name;
                            end
                        case 'sessions'
                            if sts
                                result{end+1} = BIDS.subjects(i).session;
                            end
                        case 'data'
                            if sts && isfield(d(k),'filename')
                                if strcmp(mods{j},'other')
                                    result{end+1} = fullfile(BIDS.subjects(i).path,d(k).filename);
                                else
                                    result{end+1} = fullfile(BIDS.subjects(i).path,mods{j},d(k).filename);
                                end
                            end
                        case 'metadata'
                            if sts && isfield(d(k),'filename')
                                f = fullfile(BIDS.subjects(i).path,mods{j},d(k).filename);
                                result{end+1} = get_metadata(f);
                                if ~isempty(target)
                                    try
                                        result{end} = subsref(result{end},target);
                                    catch
                                        warning('Non-existent field for metadata.');
                                        result{end} = [];
                                    end
                                end
                            end
%                             if sts && isfield(d(k),'meta')
%                                 result{end+1} = d(k).meta;
%                             end
                        case 'runs'
                            if sts && isfield(d(k),'run')
                                result{end+1} = d(k).run;
                            end
                        case 'tasks'
                            if sts && isfield(d(k),'task')
                                result{end+1} = d(k).task;
                            end
                        case 'types'
                            if sts && isfield(d(k),'type')
                                result{end+1} = d(k).type;
                            end
                    end
                end
            end
        end
        %-Postprocessing output variable
        switch query
            case 'subjects'
                result = unique(result);
                result = regexprep(result,'^[a-zA-Z0-9]+-','');
            case 'sessions'
                result = unique(result);
                result = regexprep(result,'^[a-zA-Z0-9]+-','');
                result(cellfun('isempty',result)) = [];
            case 'data'
                result = result';
            case 'metadata'
                if numel(result) == 1
                    result = result{1};
                end
            case {'tasks','runs','types'}
                result = unique(result);
                result(cellfun('isempty',result)) = [];
        end
    otherwise
        error('Unable to perform BIDS query.');
end

%==========================================================================
%-Parse BIDS query
%==========================================================================
function query = parse_query(query)
if numel(query) == 1 && isstruct(query{1})
    query = [fieldnames(query{1}), struct2cell(query{1})];
else
    if mod(numel(query),2)
        error('Invalid input syntax.');
    end
    query = reshape(query,2,[])';
end
for i=1:size(query,1)
    if ischar(query{i,2})
        query{i,2} = cellstr(query{i,2});
    end
    for j=1:numel(query{i,2})
        if iscellstr(query{i,2})
            query{i,2}{j} = regexprep(query{i,2}{j},sprintf('^%s-',query{i,1}),'');
        end
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% TEST_THICK_RING_G_NMNN: data function for Neumann boundary condition.

function g = test_thick_ring_g_nmnn (x, y, z, ind)
  [theta, r] = cart2pol (x, y);
  switch ind
    case 1
      g = -cos (theta) .* exp (x) .* cos (z) .* (sin (x.*y) + y .* cos (x.*y)) -...
            sin (theta) .* x .* exp (x) .* cos (x.*y) .* cos (z);
    case 2
      g = cos (theta) .* exp (x) .* cos (z) .* (sin (x.*y) + y .* cos (x.*y)) +...
            sin (theta) .* x .* exp (x) .* cos (x.*y) .* cos (z);
    case 3
      g = -x .* exp (x) .* cos (x.*y) .* cos (z);
    case 4
      g = -exp (x) .* cos (z) .* (sin (x.*y) + y .* cos (x.*y));
    case 5
      g = exp (x) .* sin (x.*y) .* sin (z);
    case 6
      g = -exp (x) .* sin (x.*y) .* sin (z);
    otherwise
      error ('g_nmnn: unknown reference number');
  end
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
fprintf('\nHere we train an RBM with Continuous inputs (2D Gaussian Mixture).\n');

load('gaussianData.mat');
trainData = data; clear data;
testData = testdata; clear testdata;

arch = struct('size', [2,25], 'classifier',true, 'inputType','gaussian');

opts = {'eta',.01, ...
		'batchSz',64, ...
		'nEpoch',50, ...
		'sparse',0.002, ...
		'displayEvery',10};
%  		'visFun',@visGaussianClassLearning};
arch.opts = opts;

clear r;
r = rbm(arch);
r = r.train(trainData,labels);
[pred,classError,misClass] = r.classify(testData,testlabels);

clf;
scatter(testData(:,1),testData(:,2),[],pred,'.'); colormap lines(3)
hold on
m = plot(testData(misClass,1),testData(misClass,2),'ko');
hold off;
title(sprintf('Predicted classes -- Error=%1.2f %%',classError*100));
legend(m,'Misclassified');

author: dustinstansbury

repo: medal


text:
function [hdr,otherendian] = spm_read_hdr(fname)
% Read (SPM customised) Analyze header
% FORMAT [hdr,otherendian] = spm_read_hdr(fname)
% fname       - .hdr filename
% hdr         - structure containing Analyze header
% otherendian - byte swapping necessary flag
%_______________________________________________________________________
% @(#)spm_read_hdr.m	2.2 John Ashburner 03/07/17

fid         = fopen(fname,'r','native');
otherendian = 0;
if (fid > 0)
	dime = read_dime(fid);
	if dime.dim(1)<0 | dime.dim(1)>15, % Appears to be other-endian
		% Re-open other-endian
		fclose(fid);
		if spm_platform('bigend'), fid = fopen(fname,'r','ieee-le');
		else,                      fid = fopen(fname,'r','ieee-be'); end;
		otherendian = 1;
		dime = read_dime(fid);
	end;
	hk       = read_hk(fid);
	hist     = read_hist(fid);
	hdr.hk   = hk;
	hdr.dime = dime;
	hdr.hist = hist;

	% SPM specific bit - unused
	%if hdr.hk.sizeof_hdr > 348,
	%	spmf = read_spmf(fid,dime.dim(5));
	%	if ~isempty(spmf),
	%		hdr.spmf = spmf;
	%	end;
	%end;

	fclose(fid);
else,
	hdr = [];
	otherendian = NaN;
	%error(['Problem opening header file (' fopen(fid) ').']);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function hk = read_hk(fid)
% read (struct) header_key
%-----------------------------------------------------------------------
fseek(fid,0,'bof');
hk.sizeof_hdr 		= fread(fid,1,'int32');
hk.data_type  		= mysetstr(fread(fid,10,'uchar'))';
hk.db_name    		= mysetstr(fread(fid,18,'uchar'))';
hk.extents    		= fread(fid,1,'int32');
hk.session_error	= fread(fid,1,'int16');
hk.regular			= mysetstr(fread(fid,1,'uchar'))';
hk.hkey_un0			= mysetstr(fread(fid,1,'uchar'))';
if isempty(hk.hkey_un0), error(['Problem reading "hk" of header file (' fopen(fid) ').']); end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function dime = read_dime(fid)
% read (struct) image_dimension
%-----------------------------------------------------------------------
fseek(fid,40,'bof');
dime.dim		= fread(fid,8,'int16')';
dime.vox_units	= mysetstr(fread(fid,4,'uchar'))';
dime.cal_units	= mysetstr(fread(fid,8,'uchar'))';
dime.unused1	= fread(fid,1,'int16');
dime.datatype	= fread(fid,1,'int16');
dime.bitpix		= fread(fid,1,'int16');
dime.dim_un0	= fread(fid,1,'int16');
dime.pixdim		= fread(fid,8,'float')';
dime.vox_offset	= fread(fid,1,'float');
dime.funused1	= fread(fid,1,'float');
dime.funused2	= fread(fid,1,'float');
dime.funused3	= fread(fid,1,'float');
dime.cal_max	= fread(fid,1,'float');
dime.cal_min	= fread(fid,1,'float');
dime.compressed	= fread(fid,1,'int32');
dime.verified	= fread(fid,1,'int32');
dime.glmax		= fread(fid,1,'int32');
dime.glmin		= fread(fid,1,'int32');
if isempty(dime.glmin), error(['Problem reading "dime" of header file (' fopen(fid) ').']); end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function hist = read_hist(fid)
% read (struct) data_history
%-----------------------------------------------------------------------
fseek(fid,148,'bof');
hist.descrip	= mysetstr(fread(fid,80,'uchar'))';
hist.aux_file	= mysetstr(fread(fid,24,'uchar'))';
hist.orient		= fread(fid,1,'uchar');
hist.origin		= fread(fid,5,'int16')';
hist.generated	= mysetstr(fread(fid,10,'uchar'))';
hist.scannum	= mysetstr(fread(fid,10,'uchar'))';
hist.patient_id	= mysetstr(fread(fid,10,'uchar'))';
hist.exp_date	= mysetstr(fread(fid,10,'uchar'))';
hist.exp_time	= mysetstr(fread(fid,10,'uchar'))';
hist.hist_un0	= mysetstr(fread(fid,3,'uchar'))';
hist.views		= fread(fid,1,'int32');
hist.vols_added	= fread(fid,1,'int32');
hist.start_field= fread(fid,1,'int32');
hist.field_skip	= fread(fid,1,'int32');
hist.omax		= fread(fid,1,'int32');
hist.omin		= fread(fid,1,'int32');
hist.smax		= fread(fid,1,'int32');
hist.smin		= fread(fid,1,'int32');
if isempty(hist.smin), error(['Problem reading "hist" of header file (' fopen(fid) ').']); end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function spmf = read_spmf(fid,n)
% Read SPM specific fields
% This bit may be used in the future for extending the Analyze header.

fseek(fid,348,'bof');
mgc = fread(fid,1,'int32');    % Magic number
if mgc ~= 20020417, spmf = []; return; end;

for j=1:n,
	spmf(j).mat    = fread(fid,16,'double'); % Orientation information
	spmf(j).unused = fread(fid,384,'uchar'); % Extra unused stuff
	if length(spmf(j).unused)<384,
		error(['Problem reading "spmf" of header file (' fopen(fid) ').']);
	end;
 	spmf(j).mat = reshape(spmf(j).mat,[4 4]);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function out = mysetstr(in)
tmp = find(in == 0);
tmp = min([min(tmp) length(in)]);
out = setstr([in(1:tmp)' zeros(1,length(in)-(tmp))])';
return;
%_______________________________________________________________________
%_______________________________________________________________________


text:
function img = flowToColor(flow, varargin)

%  flowToColor(flow, maxFlow) flowToColor color codes flow field, normalize
%  based on specified value,
%
%  flowToColor(flow) flowToColor color codes flow field, normalize
%  based on maximum flow present otherwise

%   According to the c++ source code of Daniel Scharstein
%   Contact: schar@middlebury.edu

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-31 18:33:30 (Wed, 31 Oct 2006) $

% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%
%                         All Rights Reserved
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for any purpose other than its incorporation into a
% commercial product is hereby granted without fee, provided that the
% above copyright notice appear in all copies and that both that
% copyright notice and this permission notice appear in supporting
% documentation, and that the name of the author and Brown University not be used in
% advertising or publicity pertaining to distribution of the software
% without specific, written prior permission.
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
% INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
% PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR BROWN UNIVERSITY BE LIABLE FOR
% ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

UNKNOWN_FLOW_THRESH = 1e9;
UNKNOWN_FLOW = 1e10;            %

[height widht nBands] = size(flow);

if nBands ~= 2
    error('flowToColor: image must have two bands');
end;

u = flow(:,:,1);
v = flow(:,:,2);

maxu = -999;
maxv = -999;

minu = 999;
minv = 999;
maxrad = -1;

% fix unknown flow
idxUnknown = (abs(u)> UNKNOWN_FLOW_THRESH) | (abs(v)> UNKNOWN_FLOW_THRESH) ;
u(idxUnknown) = 0;
v(idxUnknown) = 0;

maxu = max(maxu, max(u(:)));
minu = min(minu, min(u(:)));

maxv = max(maxv, max(v(:)));
minv = min(minv, min(v(:)));

rad = sqrt(u.^2+v.^2);
maxrad = max(maxrad, max(rad(:)));

fprintf('max flow: %.4f flow range: u = %.3f .. %.3f; v = %.3f .. %.3f\n', maxrad, minu, maxu, minv, maxv);

if isempty(varargin) ==0
    maxFlow = varargin{1};
    if maxFlow > 0
        maxrad = maxFlow;
    end;
end;

u = u/(maxrad+eps);
v = v/(maxrad+eps);

% compute color

img = computeColor(u, v);

% unknown flow
IDX = repmat(idxUnknown, [1 1 3]);
img(IDX) = 0;

author: kristinbranson

repo: JAABA


text:
function element_node = grid_t6_element ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_T6_ELEMENT produces a grid of pairs of 6 node triangles.
%
%  Example:
%
%    Input:
%
%      NELEMX = 3, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NODE =
%         1,  3, 15,  2,  9,  8;
%        17, 15,  3, 16,  9, 10;
%         3,  5, 17,  4, 11, 10;
%        19, 17,  5, 18, 11, 12;
%         5,  7, 19,  6, 13, 12;
%        21, 19,  7, 20, 13, 14;
%        15, 17, 29, 16, 23, 22;
%        31, 29, 17, 30, 23, 24;
%        17, 19, 31, 18, 25, 24;
%        33, 31, 19, 32, 25, 26;
%        19, 21, 33, 20, 27, 26;
%        35, 33, 21, 34, 27, 28.
%
%  Grid:
%
%   29-30-31-32-33-34-35
%    |\ 8  |\10  |\12  |
%    | \   | \   | \   |
%   22 23 24 25 26 27 28
%    |   \ |   \ |   \ |
%    |  7 \|  9 \| 11 \|
%   15-16-17-18-19-20-21
%    |\ 2  |\ 4  |\ 6  |
%    | \   | \   | \   |
%    8  9 10 11 12 13 14
%    |   \ |   \ |   \ |
%    |  1 \|  3 \|  5 \|
%    1--2--3--4--5--6--7
%
%  Reference Element T6:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  6  5
%    |  |   \
%    |  |    \
%    0  1--4--2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    2 * NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NODE(6,2*NELEMX*NELEMY), the nodes that form
%    each element.
%

%
%  Node labeling:
%
%    NW---N--NE
%     | \     |
%     W   C   E
%     |    \  |
%    SW---S--SE
%
  element = 0;

  for j = 1 : nelemy
    for i = 1 : nelemx

      sw = 2 * ( j - 1 )  * ( 2 * nelemx + 1 ) + 2 * ( i - 1 ) + 1;
      w  = sw +               2 * nelemx + 1;
      nw = sw +         2 * ( 2 * nelemx + 1 );

      s  = sw + 1;
      c  = sw + 1 +               2 * nelemx + 1;
      n  = sw + 1 +         2 * ( 2 * nelemx + 1 );

      se = sw + 2;
      e  = sw + 2 +               2 * nelemx + 1;
      ne = sw + 2 +         2 * ( 2 * nelemx + 1 );

      element = element + 1;

      element_node(1,element) = sw;
      element_node(2,element) = se;
      element_node(3,element) = nw;
      element_node(4,element) = s;
      element_node(5,element) = c;
      element_node(6,element) = w;

      element = element + 1;

      element_node(1,element) = ne;
      element_node(2,element) = nw;
      element_node(3,element) = se;
      element_node(4,element) = n;
      element_node(5,element) = c;
      element_node(6,element) = e;

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Load AOI from file
function ApToolLoadImage(filenames,hObject,handles,varargin)

if isfield(handles,'H')
    delete(handles.H);
end

if ((nargin<1)||isempty(filenames)) % If no filename was give, use dialog box to ask for one
    % Recall last interactively selected path used
    if ispref('matlab_sar_toolbox','last_used_directory')
        pathname = getpref('matlab_sar_toolbox','last_used_directory');
        if ~ischar(pathname)||~exist(pathname,'dir')
            pathname = pwd;
        end
    else
        pathname = pwd;
    end
    [filenames,pathname]=uigetfile(sar_file_extensions('amplitude'),...
        'Open SAR Data File',pathname,'MultiSelect','on');
    if(iscell(filenames)), filenames=sort(filenames); end
    setpref('matlab_sar_toolbox','last_used_directory',pathname);
    if isfield(handles,'phdhandle')
        handles = rmfield(handles,'phdhandle');
    end
else % Path was already passed in with filenames
    pathname='';
end
fullfilename = {}; % Reset if files were previously open
if(iscell(filenames)) % Multiple files requested
    for j=1:length(filenames)
        fullfilename{j}=[pathname filenames{j}];
    end
elseif(filenames)
    fullfilename{1}=[pathname filenames];
else % filename=0.  Cancel was pressed, instead of a file being chosen.
    return;
end
if ~isempty(pathname)
    setpref('matlab_sar_toolbox','last_used_directory',pathname); %store path
end

% Open readers for all files, and build cell array of all readers
reader_obj ={};
sf_ind = []; fn_ind = []; % Indices to filenames and subimages for each reader object
for i = 1:length(fullfilename)
    newreader = open_reader(fullfilename{i});
    fn_ind = [fn_ind i*ones(size(newreader))];
    sf_ind = [sf_ind 1:numel(newreader)];
    reader_obj = [reader_obj newreader];
end
clear newreader;

[readerobj_indices,reader_obj] = group_reader_objs_by_pol(reader_obj); % Consolidate reader objects into polarimetric sets

%get metadata
if (iscell(reader_obj))
    meta=reader_obj{1}.get_meta(); % Assume first frame for getting datasize
    % Adjust later if frame changes
else
    meta=reader_obj.get_meta();
end
nx = double(meta.ImageData.NumCols);
ny = double(meta.ImageData.NumRows);

%set filename or corename as the title of the GUI
if ~iscell(filenames)
    temp = filenames;
    clear filenames;
    filenames{1} = temp;
end
[~,name,ext] = fileparts(filenames{1});
fname = strcat(name,ext);
set(handles.figure1,'Name',['ApertureTool, Filename: ', fname]);

% Select AOI and frame number if necessary
if isempty(varargin) || isempty(varargin{1}) || varargin{1}(1) == 0
    %determine if we'll just load the image or if we'll load an AOI
    %selection MITM_Viewer
    if (nx <= handles.MaxLoadSize && ny <= handles.MaxLoadSize)
        %load image directly into tool
        aoi = [1 1 nx ny];
        FrameNum = 1;
    else
        %launch a MITM viewer with an AOI selection box
         ButtonName = questdlg(sprintf('%s %i x %i\n%s %i x %i\n%s',...
            'Suggested AOI dimensions are less than',handles.MaxLoadSize,handles.MaxLoadSize,...
            'Current AOI dimensions are',nx,ny,...
            'Do you want to redefine the AOI?'), ...
            'AptertureTool AOI', ...
            'Yes', 'No', 'Yes');
        if isempty(ButtonName), return; end %window was closed CANCEL
        if strcmpi(ButtonName,'No')
            aoi = [1 1 nx ny]; %use AOI anyway
            FrameNum = 1;
        else
            [aoi, FrameNum] = mitm_viewer(fullfilename,'mode','aoi','closeAfterSelect',true);
            FrameNum = find(cellfun(@(x) isequal(FrameNum,sf_ind(x)), readerobj_indices));
        end
        if isempty(aoi), return; end % Cancel was selected
    end
else
    aoi = varargin{1};
    if nargin>4 && ~isempty(varargin{2}) && all(varargin{2} ~= 0)
        FrameNum = find(cellfun(@(x) isequal(varargin{2},sf_ind(x)), readerobj_indices));
    else
        FrameNum = 1;
    end
end
if (iscell(reader_obj))
    meta=reader_obj{FrameNum}.get_meta(); % Assume first frame for getting datasize
end

set(handles.nx,'String',nx);
set(handles.ny,'String',ny);

%Draw metaicon
MetaIcon(meta,'handle',handles.metaicon);

% Handle complex data normalization parameters from metadata
% Set data normalization controls based on metadata
if isfield(meta,'Grid') && all(isfield(meta.Grid,{'Row','Col'})) && ...
        isfield(meta.Grid.Row,'WgtType') && ...
        isfield(meta.Grid.Col,'WgtType') && ...
        ((strcmpi(meta.Grid.Col.WgtType,'UNIFORM') && ...
        strcmpi(meta.Grid.Row.WgtType,'UNIFORM')) || ...
        (isfield(meta.Grid.Row.WgtType,'WindowName') && ...
        isfield(meta.Grid.Col.WgtType,'WindowName') && ...
        strcmpi(meta.Grid.Col.WgtType.WindowName,'UNIFORM') && ...
        strcmpi(meta.Grid.Row.WgtType.WindowName,'UNIFORM')))
    % If we already know that data is uniform weighted, this is not an option.
    set(handles.UniformWeighting,'Value',true);
    set(handles.UniformWeighting,'Enable','off');
else
    set(handles.UniformWeighting,'Value',false);
    set(handles.UniformWeighting,'Enable','on');
end
% Set default values for required frequency support deskew values
handles.isvalid_row_deltakcoapoly = isfield(meta,'Grid') && ...
    isfield(meta.Grid,'Row') && isfield(meta.Grid.Row,'DeltaKCOAPoly') && ...
    ~isempty(meta.Grid.Row.DeltaKCOAPoly) && ...
    all(isfinite(meta.Grid.Row.DeltaKCOAPoly(:)));
if handles.isvalid_row_deltakcoapoly && any(meta.Grid.Row.DeltaKCOAPoly(:)~=0)
    set(handles.DeskewFast,'Enable','on');
else
    set(handles.DeskewFast,'Enable','off');
    meta.Grid.Row.DeltaKCOAPoly = 0; % Give default value (no shift)
end
handles.isvalid_col_deltakcoapoly = isfield(meta,'Grid') && ...
    isfield(meta.Grid,'Col') && isfield(meta.Grid.Col,'DeltaKCOAPoly') && ...
    ~isempty(meta.Grid.Col.DeltaKCOAPoly) && ...
    all(isfinite(meta.Grid.Col.DeltaKCOAPoly(:)));
if handles.isvalid_col_deltakcoapoly && any(meta.Grid.Col.DeltaKCOAPoly(:)~=0)
    set(handles.DeskewSlow,'Enable','on');
else
    set(handles.DeskewSlow,'Enable','off');
    meta.Grid.Col.DeltaKCOAPoly = 0; % Give default value (no shift)
end
% Is it possible to deskew frequency support in both dimensions at once?
if isscalar(meta.Grid.Col.DeltaKCOAPoly) && ...
        isscalar(meta.Grid.Row.DeltaKCOAPoly)
    set(handles.DeskewSlow,'Value',handles.isvalid_col_deltakcoapoly); % Center both dimensions
    set(handles.DeskewFast,'Value',handles.isvalid_row_deltakcoapoly);
else % Make user select dimension to deskew
    set(handles.DeskewSlow,'Value',handles.isvalid_col_deltakcoapoly && ...
        all(meta.Grid.Col.DeltaKCOAPoly(:) == 0));
    set(handles.DeskewFast,'Value',handles.isvalid_row_deltakcoapoly && ...
        all(meta.Grid.Row.DeltaKCOAPoly(:) == 0));
end
handles.manual_offset = [0 0];
% Determine weighting function
try
    handles.weight_fun_az = sicdweight2fun(meta.Grid.Col);
    handles.weight_fun_rng = sicdweight2fun(meta.Grid.Row);
    handles.weighting_computed_from_data = false;
catch % Unable to determine weighting functions from metadata.
    % Estimate weighting from complex data instead
    if iscell(reader_obj)
        handles.weight_fun_az = estimate_weighting_file(reader_obj{FrameNum}, 1);
        handles.weight_fun_rng = estimate_weighting_file(reader_obj{FrameNum}, 2);
    else
        handles.weight_fun_az = estimate_weighting_file(reader_obj, 1);
        handles.weight_fun_rng = estimate_weighting_file(reader_obj, 2);
    end
    handles.weighting_computed_from_data = true;
end
% Determine if we have metadata to compute zeropad accurately
if isfield(meta,'Grid') && all(isfield(meta.Grid,{'Row','Col'})) && ...
        all(isfield(meta.Grid.Row,{'SS','ImpRespBW'})) && ...
        all(isfield(meta.Grid.Col,{'SS','ImpRespBW'}))
    handles.AzPad = max(1,1/(meta.Grid.Col.SS*meta.Grid.Col.ImpRespBW));
    handles.RnPad = max(1,1/(meta.Grid.Row.SS*meta.Grid.Row.ImpRespBW));
else % Otherwise, just default to full data.
    handles.AzPad = 1;
    handles.RnPad = 1;
end % User can always manual adjust zeropad area

% Use correct FFT function based on data
% Sgn of -1 means that a forward FFT is the correct transform to move from
% image domain to spatial frequency.  However, we want to flip the
% direction, so that high RF frequencies are at the top (low indices) when
% displayed in MATLAB, and time/angle increases toward the right (higher
% indices).  So we swap the FFT/IFFTs which, in effect, flips left/right
% and up/down in spatial frequency.
if isfield(meta,'Grid') && all(isfield(meta.Grid,{'Row','Col'})) && ...
        isfield(meta.Grid.Row,'Sgn') && isfield(meta.Grid.Col,'Sgn') && ...
        meta.Grid.Row.Sgn > 0 && meta.Grid.Col.Sgn > 0
    handles.fft_sp = @fft2; % Function to transform from image to spatial frequency
    handles.fft_im = @ifft2; % Function to transfrom from spatial frequency to image
elseif isfield(meta,'Grid') && all(isfield(meta.Grid,{'Row','Col'})) && ...
        isfield(meta.Grid.Row,'Sgn') && isfield(meta.Grid.Col,'Sgn') && ...
        meta.Grid.Row.Sgn ~= meta.Grid.Col.Sgn
    error('APERTURETOOL:MIXED_FFT_SGN','Mixed FFT signs not handled.');
else
    handles.fft_sp = @ifft2; % Function to transform from image to spatial frequency
    handles.fft_im = @fft2; % Function to transfrom from spatial frequency to image
end

%get polarization information
if isfield(meta,'ImageFormation') && isfield(meta.ImageFormation,'TxRcvPolarizationProc')
    handles.TxRcvPolarizationProc = meta.ImageFormation.TxRcvPolarizationProc;
    if ~iscell(handles.TxRcvPolarizationProc)
        handles.TxRcvPolarizationProc = {handles.TxRcvPolarizationProc};
        set(handles.image_rep_combo,'value',1);
    else
        set(handles.image_rep_combo,'value',...
            find(strcmpi('Pauli',get(handles.image_rep_combo,'string'))));
    end
    handles.polar_PHD_idx = 1;
    if numel(handles.TxRcvPolarizationProc) == 1
        set(handles.image_rep_combo,'value',1,'visible','off');
        set(handles.PHDPolCombo,'visible','off');
    else
        set(handles.image_rep_combo,'value',2,'visible','on');
        set(handles.PHDPolCombo,'visible','on');
    end
    %set strings in PHDPolCombo
    set(handles.PHDPolCombo,'String',meta.ImageFormation.TxRcvPolarizationProc);
else
    set(handles.image_rep_combo,'value',1,'visible','off');
    set(handles.PHDPolCombo,'visible','off');
end

% Read data
handles.complex_data = double(reader_obj{FrameNum}.read_chip(...
    [aoi(1) aoi(1)+aoi(3)-1],[aoi(2) aoi(2)+aoi(4)-1]));
reader_obj{FrameNum}.close();
% Condition data so that imcontrast control works better
handles.complex_data = 10*handles.complex_data./mean(handles.complex_data(:));
handles.aoi = aoi; % store aoi for later
handles.segment = FrameNum;

% This is derived redundant information, but we compute it once, so we
% don't need to do this in multiple places:
handles.ZpLimsAz = round((aoi(3)/2) * (1 - 1/handles.AzPad)) + 1;
handles.ZpLimsAz(2) = aoi(3) - handles.ZpLimsAz + 1;
handles.ZpWidthAz = diff(handles.ZpLimsAz)+1;
handles.ZpLimsRn = round((aoi(4)/2) * (1 - 1/handles.RnPad)) + 1;
handles.ZpLimsRn(2) = aoi(4) - handles.ZpLimsRn + 1;
handles.ZpWidthRn = diff(handles.ZpLimsRn)+1;

% update metadata
% chip info
meta.ImageData.FirstRow = meta.ImageData.FirstRow + aoi(2) - 1;
meta.ImageData.FirstCol = meta.ImageData.FirstCol + aoi(1) - 1;
meta.ImageData.NumRows = aoi(4);
meta.ImageData.NumCols = aoi(3);
try
    meta = add_sicd_corners(meta);
end
% If ImpRespWid is not defined, but ImpRespBW is, we just use what it would
% be for uniform weighting (since we're probably going to want to apply a
% uniform weighting to the data anyway.)
set(handles.PHDZoomCombo,'value',1); % Reset to "Full"
resolution_defined = true;
if ~isfield(meta.Grid.Col,'ImpRespWid')
    if isfield(meta.Grid.Col,'ImpRespBW')
        meta.Grid.Col.ImpRespWid = .886/meta.Grid.Col.ImpRespBW;
    else % Just set to zero if we don't know it
        meta.Grid.Col.ImpRespWid = 0;
        meta.Grid.Col.SS = 0;
        resolution_defined = false;
    end
end
if ~isfield(meta.Grid.Row,'ImpRespWid')
    if isfield(meta.Grid.Row,'ImpRespBW')
        meta.Grid.Row.ImpRespWid = .886/meta.Grid.Row.ImpRespBW;
    else % Just set to zero if we don't know it
        meta.Grid.Row.ImpRespWid = 0;
        meta.Grid.Row.SS = 0;
        resolution_defined = false;
    end
end

%set initial SS, Ap Time and Ap Bandwidth

handles.meta = meta; % store metadata for later

ApToolprocess_phd(hObject,handles); % Apply PHD options
handles = guidata(hObject); % process_phd change handles (We may fix this later.)

%set initial magnifications to fit
mag = handles.apiSP1.findFitMag();
handles.apiSP1.setMagnification(mag);
set(handles.ImageZoomCombo,'Value',1);
strings = get(handles.ImageZoomCombo,'String');
strings{1} = sprintf('Fit (%d%%)',round(mag*100));
set(handles.ImageZoomCombo,'String',strings);

% Draw imrect on phd
if isfield(handles,'H') && ishandle(handles.H)
    delete(handles.H); % Delete any previous imrects
end
handles.H = imrect(handles.phd, [handles.ZpLimsAz(1), handles.ZpLimsRn(1),...
    handles.ZpWidthAz handles.ZpWidthRn]);
if verLessThan('matlab', '7.6')
    disp('Warning. Matlab version does not support imrect constraint function');
else
    setPositionConstraintFcn(handles.H,...
        makeConstrainToRectFcn('imrect',handles.ZpLimsAz,handles.ZpLimsRn));
end
api_handle = iptgetapi(handles.H);
api_handle.addNewPositionCallback(@(x) ApToolnewpos(x,hObject));

% Now enable relevant controls
set(handles.Play,'enable','on');
set(handles.Last,'enable','on');
set(handles.First,'enable','on');
set(handles.Next,'enable','off'); % Full aperture is not a frame, so
set(handles.Prev,'enable','off'); % no next or previous
set(handles.Recording,'Enable','on');
set(handles.SaveImage,'Enable','on');
set(handles.ManualAdjust,'enable','on');

% Update handles structure
guidata(hObject, handles);

% Update resolution fields
ApToolnewpos(getPosition(handles.H),hObject);

author: ngageoint