text:
% crossf() - Returns estimates and plots event-related coherence (ERCOH)
%        between two input data time series (X,Y). A lower panel (optionally)
%        shows the coherence phase difference between the processes.
%        In this panel, output by   > crossf(X,Y,...);
%            90 degrees (orange) means X leads Y by a quarter cycle.
%           -90 degrees (blue)   means Y leads X by a quarter cycle.
%        Coherence phase units may be radians, degrees, or msec.
%        Click on any subplot to view separately and zoom in/out.
%
% Function description:
%        Uses EITHER fixed-window, zero-padded FFTs (fastest) OR constant-Q
%        0-padded wavelet DFTs (more even sensitivity across frequencies),
%        both Hanning-tapered.  Output frequency spacing is the lowest
%        frequency ('srate'/'winsize') divided by the 'padratio'.
%
%        If an 'alpha' value is given, then bootstrap statistics are
%        computed (from a distribution of 'naccu' {200} surrogate baseline
%        data epochs) for the baseline epoch, and non-significant features
%        of the output plots are zeroed (and shown in green). The baseline
%        epoch is all windows with center latencies < the given 'baseline'
%        value, or if 'baseboot' is 1, the whole epoch.
% Usage:
%        >> [coh,mcoh,timesout,freqsout,cohboot,cohangles] ...
%                       = crossf(X,Y,frames,tlimits,srate,cycles, ...
%                                        'key1', 'val1', 'key2', val2' ...);
% Required inputs:
%       X       = first single-channel data set (1,frames*nepochs)
%       Y       = second single-channel data set (1,frames*nepochs)
%       frames  = frames per epoch                                 {default: 750}
%       tlimits = [mintime maxtime] (ms) epoch latency limits {def: [-1000 2000]}
%       srate   = data sampling rate (Hz)                          {default: 250}
%       cycles  = 0  -> Use FFTs (with constant window length)
%               = >0 -> Number of cycles in each analysis wavelet
%               = [cycles expfactor] -> if 0 < expfactor < 1,  the number
%                 of wavelet cycles expands with frequency from cycles
%                 If expfactor = 1, no expansion; if = 0, constant
%                 window length (as in FFT)                          {default: 0}
% Optional Coherence Type:
%       'type'  = ['coher'|'phasecoher'] Compute either linear coherence
%                 ('coher') or phase coherence ('phasecoher') also known
%                 as phase coupling factor' {default: 'phasecoher'}.
%       'subitc' = ['on'|'off'] subtract stimulus locked Inter-Trial Coherence
%                 from X and Y. This computes the  'intrinsic' coherence
%                 X and Y not arising from common synchronization to
%                 experimental events. See notes. {default: 'off'}
%       'shuffle' = integer indicating the number of estimates to compute
%                 bootstrap coherence based on shuffled trials. This estimates
%                 the coherence arising only from time locking of X and Y
%                 to experimental events (opposite of 'subitc')      {default: 0}
% Optional Detrend:
%       'detret' = ['on'|'off'], Linearly detrend data within epochs {def: 'off'}
%       'detrep' = ['on'|'off'], Linearly detrend data across trials {def: 'off'}
%
% Optional FFT/DFT:
%       'winsize'  = If cycles==0: data subwindow length (fastest, 2^n<frames);
%                    if cycles >0: *longest* window length to use. This
%                    determines the lowest output frequency  {default: ~frames/8}
%       'timesout' = Number of output latencies (int<frames-winsize)   {def: 200}
%       'padratio' = FFTlength/winsize (2^k)                         {default: 2}
%                    Multiplies the number of output frequencies by
%                    dividing their spacing. When cycles==0, frequency
%                    spacing is (low_frequency/padratio).
%       'maxfreq'  = Maximum frequency (Hz) to plot (& output if cycles>0)
%                    If cycles==0, all FFT frequencies are output  {default: 50}
%       'baseline' = Coherence baseline end latency (ms). NaN -> No baseline
%                      {default:NaN}
%       'powbase'  = Baseline spectrum to log-subtract      {default: from data}
%
% Optional Bootstrap:
%       'alpha'    = If non-0, compute two-tailed bootstrap significance prob.
%                    level. Show non-signif output values as green. {def: 0}
%       'naccu'    = Number of bootstrap replications to compute {def: 200}
%       'boottype' = ['times'|'timestrials'] Bootstrap type: Either shuffle
%                    windows ('times') or windows and trials ('timestrials')
%                    Option 'timestrials' requires more memory  {default: 'times'}
%       'memory'   = ['low'|'high'] 'low' -> decrease memory use {default: 'high'}
%       'baseboot' = Extent of bootstrap shuffling (0=to 'baseline'; 1=whole epoch)
%                    If no baseline is given (NaN), extent of bootstrap shuffling
%                    is the whole epoch                         {default: 0}
%       'rboot'    = Input bootstrap coherence limits (e.g., from crossf())
%                    The bootstrap type should be identical to that used
%                    to obtain the input limits. {default: compute from data}
% Optional Scalp Map:
%       'topovec'  = (2,nchans) matrix, plot scalp maps to plot {default: []}
%                    ELSE (c1,c2), plot two cartoons showing channel locations.
%       'elocs'    = Electrode location structure or file for scalp map
%                    {default: none}
%       'chaninfo' = Electrode location additional information (nose position...)
%                    {default: none}
%
% Optional Plot and Compute Features:
%       'compute'   = ['matlab'|'c'] Use C subroutines to speed up the
%                     computation (currently unimplemented) {def: 'matlab'}
%       'savecoher' - [0|1] 1 --> Accumulate the individual trial coherence
%                     vectors; output them as cohangles {default: 0 = off}
%       'plotamp'   = ['on'|'off'], Plot coherence magnitude    {def: 'on'}
%       'maxamp'    = [real] Set the maximum for the amp. scale {def: auto}
%       'plotphase' = ['on'|'off'], Plot coherence phase angle  {def: 'on'}
%       'angleunit' = Phase units: 'ms' -> msec, 'deg' -> degrees,
%                     or 'rad' -> radians                  {default: 'deg'}
%       'title'     = Optional figure title                {default:  none}
%       'vert'      = Latencies to mark with a dotted vertical line
%                                                           {default: none}
%       'linewidth' = Line width for marktimes traces (thick=2, thin=1)
%                                                              {default: 2}
%       'cmax'      = Maximum amplitude for color scale  {def: data limits}
%       'axesfont'  = Axes font size                          {default: 10}
%       'titlefont' = Title font size                          {default: 8}
%
% Outputs:
%       coh         = Matrix (nfreqs,timesout) of coherence magnitudes
%       mcoh        = Vector of mean baseline coherence at each frequency
%       timesout    = Vector of output latencies (window centers) (ms).
%       freqsout    = Vector of frequency bin centers (Hz).
%       cohboot     = Matrix (nfreqs,2) of [lower;upper] coher signif. limits
%                     if 'boottype' is 'trials',  (nfreqs,timesout, 2)
%       cohangle    = (nfreqs,timesout) matrix of coherence angles (in radians)
%       cohangles   = (nfreqs,timesout,trials) matrix of single-trial coherence
%                      angles (in radians), saved and output only if 'savecoher',1
%
% Plot description:
%   Assuming both 'plotamp' and 'plotphase' options are 'on' (=default), the upper panel
%   presents the magnitude of either phase coherence or linear coherence, depending on
%   the 'type' parameter (above). The lower panel presents the coherence phase difference
%   (in degrees). Click on any plot to pop up a new window (using 'axcopy()').
%   -- The upper left marginal panel shows mean coherence during the baseline period
%      (blue), and when significance is set, the significance threshold (dotted black-green).
%   -- The horizontal panel under the coherence magnitude image indicates the maximum
%      (green) and minimum (blue) coherence values across all frequencies. When significance
%      is set (using option 'trials' for 'boottype'), an additional curve indicates the
%      significance threshold (dotted black-green).
%
% Notes: 1) When cycles==0, nfreqs is total number of FFT frequencies.
%        2) As noted above: 'blue' coherence angle -> X leads Y; 'red' -> Y leads X
%        3) The 'boottype' should be ideally 'timesframes', but this creates high
%           memory demands, so the 'times' method must be used in many cases.
%        4) If 'boottype' is 'trials', the average of the complex bootstrap
%           is subtracted from the coherence to compensate for phase differences
%           (the average is also subtracted from the bootstrap distribution).
%           For other bootstraps, this is not necessary since the phase is random.
%        5) If baseline is non-NaN, the baseline is subtracted from
%           the complex coherence. On the left hand side of the coherence
%           amplitude image, the baseline is displayed as a magenta line
%           (if no baseline is selected, this curve represents the average
%           coherence at every given frequency).
%        6) If a out-of-memory error occurs, set the 'memory' option to 'low'
%           (Makes computation time slower; Only the 'times' bootstrap method
%           can be used in this mode).
%
% Authors: Arnaud Delorme, Sigurd Enghoff & Scott Makeig
%          CNL/Salk Institute 1998-2001; SCCN/INC/UCSD, La Jolla, 2002-
%
% See also: timef()

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

% 11-20-98 defined g.linewidth constant -sm
% 04-01-99 made number of frequencies consistent -se
% 06-29-99 fixed constant-Q freq indexing -se
% 08-13-99 added cohangle plotting -sm
% 08-20-99 made bootstrap more efficient -sm
% 08-24-99 allow nan values introduced by possible eventlock() preproc. -sm
% 03-05-2007 eventlock.m deprecated to eegalign.m. -tf
% 03-16-00 added lead/lag interpretation to help msg - sm & eric visser
% 03-16-00 added axcopy() feature -sm & tpj
% 04-20-00 fixed Rangle sign for wavelets, added verts array -sm
% 01-22-01 corrected help msg when nargin<2 -sm & arno delorme
% 01-25-02 reformated help & license, added links -ad
% 03-09-02 function restructuration -ad
%  add 'key', val arguments (+ external baseboot, baseline, color axis, angleunit...)
%  add detrending (across time and trials) + 'coher' option for amplitude coherence
%  significance only if alpha is given, ploting options in 'plotamp' and 'plotphase'
% 03-16-02 timeout automatically adjusted if too high -ad
% 04-03-02 added new options for bootstrap -ad

% Note: 3 "objects" (Tf, Coher and Boot) are handled by specific functions under Matlab
%    (Tf) function Tf = tfinit(...) - create object Time Frequency (Tf) associated with some data
%    (Tf) function [Tf, itcvals] = tfitc(...) - compute itc for the selected data
%    (Tf) function [Tf, itcvals] = tfitcpost(Tf, trials) - itc normlisation
%    (Tf) function [Tf, tmpX] = tfcomp(Tf, trials, times) - compute time freq. decomposition
%    (Coher) function Coher = coherinit(...) - initialize coherence object
%    (Coher) function Coher = cohercomp(Coher, tmpX, tmpY, trial, time) - compute coherence
%    (Coher) function Coher = cohercomppost(Coher, trials) - coherence normalization
%    (Boot) function Boot = bootinit(...) - intialize bootstrap object
%    (Boot) function Boot = bootcomp(...) - compute bootstrap
%    (Boot) function [Boot, Rbootout] = bootcomppost(...) - bootstrap normalization
% and by real objects under C++ (C++ code, incomplete)

function [R,mbase,times,freqs,Rbootout,Rangle, trialcoher, Tfx, Tfy] = crossf(X, Y, frame, tlimits, Fs, varwin, varargin)

%varwin,winsize,nwin,oversmp,maxfreq,alpha,verts,caxmax)

% ------------------------
% Commandline arg defaults:
% ------------------------
DEFAULT_ANGLEUNIT = 'deg'; % angle plotting units - 'rad', 'ms', or 'deg'
DEFAULT_EPOCH	= 750;			% Frames per epoch
DEFAULT_TIMELIM = [-1000 2000];	% Time range of epochs (ms)
DEFAULT_FS		= 250;			% Sampling frequency (Hz)
DEFAULT_NWIN	= 200;			% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;			% Fixed window length or base on cycles.

% =0: fix window length to nwin
% >0: set window length equal varwin cycles
%     bounded above by winsize, also determines
%     the min. freq. to be computed.

DEFAULT_OVERSMP	= 2;			% Number of times to oversample = vertical resolution
DEFAULT_MAXFREQ = 50;			% Maximum frequency to display (Hz)
DEFAULT_TITLE	= 'Event-Related Coherence';			% Figure title
DEFAULT_ALPHA   = NaN;			% Default two-sided significance probability threshold

if (nargin < 2)
   help crossf
   return
end

if ~iscell(X)
	if (min(size(X))~=1 | length(X)<2)
		fprintf('crossf(): X must be a row or column vector.\n');
		return
	elseif (min(size(Y))~=1 | length(Y)<2)
		fprintf('crossf(): Y must be a row or column vector.\n');
		return
	elseif (length(X) ~= length(Y))
		fprintf('crossf(): X and Y must have same length.\n');
		return
	end
end;

if (nargin < 3)
   frame = DEFAULT_EPOCH;
elseif (~isnumeric(frame) | length(frame)~=1 | frame~=round(frame))
   fprintf('crossf(): Value of frames must be an integer.\n');
   return
elseif (frame <= 0)
   fprintf('crossf(): Value of frames must be positive.\n');
   return
elseif ~iscell(X) & (rem(length(X),frame) ~= 0)
   fprintf('crossf(): Length of data vectors must be divisible by frames.\n');
   return
end

if (nargin < 4)
   tlimits = DEFAULT_TIMELIM;
elseif (~isnumeric(tlimits) | sum(size(tlimits))~=3)
   error('crossf(): Value of tlimits must be a vector containing two numbers.');
elseif (tlimits(1) >= tlimits(2))
   error('crossf(): tlimits interval must be [min,max].');
end

if (nargin < 5)
   Fs = DEFAULT_FS;
elseif (~isnumeric(Fs) | length(Fs)~=1)
   error('crossf(): Value of srate must be a number.');
elseif (Fs <= 0)
   error('crossf(): Value of srate must be positive.');
end

if (nargin < 6)
   varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) | length(varwin)>2)
   error('crossf(): Value of cycles must be a number or a (1,2) vector.');
elseif (varwin < 0)
   error('crossf(): Value of cycles must be either zero or positive.');
end

% consider structure for these arguments
% --------------------------------------
vararginori = varargin;
for index=1:length(varargin)
	if iscell(varargin{index}), varargin{index} = { varargin{index} }; end;
end;
if ~isempty(varargin)
   try, g = struct(varargin{:});
   catch, error('Argument error in the {''param'', value} sequence'); end;
else
	g = [];
end;

try, g.shuffle;    catch, g.shuffle = 0; end;
try, g.title;      catch, g.title = DEFAULT_TITLE; end;
try, g.winsize;    catch, g.winsize = max(pow2(nextpow2(frame)-3),4); end;
try, g.pad;        catch, g.pad = max(pow2(nextpow2(g.winsize)),4); end;
try, g.timesout;   catch, g.timesout = DEFAULT_NWIN; end;
try, g.padratio;   catch, g.padratio = DEFAULT_OVERSMP; end;
try, g.maxfreq;    catch, g.maxfreq = DEFAULT_MAXFREQ; end;
try, g.topovec;    catch, g.topovec = []; end;
try, g.elocs;      catch, g.elocs = ''; end;
try, g.alpha;      catch, g.alpha = DEFAULT_ALPHA; end;
try, g.marktimes;  catch, g.marktimes = []; end; % default no vertical lines
try, g.marktimes = g.vert;       catch, g.vert = []; end; % default no vertical lines
try, g.powbase;    catch, g.powbase = nan; end;
try, g.rboot;      catch, g.rboot = nan; end;
try, g.plotamp;    catch, g.plotamp = 'on'; end;
try, g.plotphase;  catch, g.plotphase  = 'on'; end;
try, g.plotbootsub;  catch, g.plotbootsub  = 'on'; end;
try, g.detrep;     catch, g.detrep = 'off'; end;
try, g.detret;     catch, g.detret = 'off'; end;
try, g.baseline;   catch, g.baseline = NaN; end;
try, g.baseboot;   catch, g.baseboot = 0; end;
try, g.linewidth;  catch, g.linewidth = 2; end;
try, g.naccu;      catch, g.naccu = 200; end;
try, g.angleunit;  catch, g.angleunit = DEFAULT_ANGLEUNIT; end;
try, g.cmax;       catch, g.cmax = 0; end; % 0=use data limits
try, g.type;       catch, g.type = 'phasecoher'; end;
try, g.boottype;   catch, g.boottype = 'times'; end;
try, g.subitc;     catch, g.subitc = 'off'; end;
try, g.memory;     catch, g.memory = 'high'; end;
try, g.compute;    catch, g.compute = 'matlab'; end;
try, g.maxamp;     catch, g.maxamp = []; end;
try, g.savecoher;  catch, g.savecoher = 0; end;
try, g.noinput;    catch, g.noinput = 'no'; end;
try, g.chaninfo;   catch, g.chaninfo = []; end;

allfields = fieldnames(g);
for index = 1:length(allfields)
	switch allfields{index}
	 case { 'shuffle' 'title' 'winsize' 'pad' 'timesout' 'padratio' 'maxfreq' 'topovec' 'elocs' 'alpha' ...
		  'marktimes' 'vert' 'powbase' 'rboot' 'plotamp' 'plotphase' 'plotbootsub' 'detrep' 'detret' ...
		  'baseline' 'baseboot' 'linewidth' 'naccu' 'angleunit' 'cmax' 'type' 'boottype' 'subitc' ...
		  'memory' 'compute' 'maxamp' 'savecoher' 'noinput' 'chaninfo' };
	  case {'plotersp' 'plotitc' }, disp(['crossf warning: timef option ''' allfields{index} ''' ignored']);
	 otherwise disp(['crossf error: unrecognized option ''' allfields{index} '''']); beep; return;
	end;
end;

g.tlimits = tlimits;
g.frame   = frame;
g.srate   = Fs;
g.cycles  = varwin(1);
if length(varwin)>1
	g.cyclesfact = varwin(2);
else
	g.cyclesfact = 1;
end;
g.type       = lower(g.type);
g.boottype   = lower(g.boottype);
g.detrep     = lower(g.detrep);
g.detret     = lower(g.detret);
g.plotphase  = lower(g.plotphase);
g.plotbootsub = lower(g.plotbootsub);
g.subitc     = lower(g.subitc);
g.plotamp    = lower(g.plotamp);
g.shuffle    = lower(g.shuffle);
g.compute    = lower(g.compute);
g.AXES_FONT  = 10;
g.TITLE_FONT = 14;

% testing arguments consistency
% -----------------------------
if (~ischar(g.title))
   error('Title must be a string.');
end

if (~isnumeric(g.winsize) | length(g.winsize)~=1 | g.winsize~=round(g.winsize))
   error('Value of winsize must be an integer number.');
elseif (g.winsize <= 0)
   error('Value of winsize must be positive.');
elseif (g.cycles == 0 & pow2(nextpow2(g.winsize)) ~= g.winsize)
   error('Value of winsize must be an integer power of two [1,2,4,8,16,...]');
elseif (g.winsize > g.frame)
   error('Value of winsize must be less than frame length.');
end

if (~isnumeric(g.timesout) | length(g.timesout)~=1 | g.timesout~=round(g.timesout))
   error('Value of timesout must be an integer number.');
elseif (g.timesout <= 0)
   error('Value of timesout must be positive.');
end
if (g.timesout > g.frame-g.winsize)
   g.timesout = g.frame-g.winsize;
   disp(['Value of timesout must be <= frame-winsize, timeout adjusted to ' int2str(g.timesout) ]);
end

if (~isnumeric(g.padratio) | length(g.padratio)~=1 | g.padratio~=round(g.padratio))
   error('Value of padratio must be an integer.');
elseif (g.padratio <= 0)
   error('Value of padratio must be positive.');
elseif (pow2(nextpow2(g.padratio)) ~= g.padratio)
   error('Value of padratio must be an integer power of two [1,2,4,8,16,...]');
end

if (~isnumeric(g.maxfreq) | length(g.maxfreq)~=1)
   error('Value of g.maxfreq must be a number.');
elseif (g.maxfreq <= 0)
   error('Value of g.maxfreq must be positive.');
elseif (g.maxfreq > Fs/2)
   fprintf('Warning: input value of g.maxfreq larger that Nyquist frequency %3.4 Hz\n\n',Fs/2);
end

if isempty(g.topovec)
   g.topovec = [];
elseif min(size(g.topovec))==1
   g.topovec = g.topovec(:);
   if size(g.topovec,1)~=2
      error('topovec must be a row or column vector.');
   end
end;

if isempty(g.elocs)
   g.elocs = '';
elseif (~ischar(g.elocs)) & ~isstruct(g.elocs)
   error('Channel location file must be a valid text file.');
end

if (~isnumeric(g.alpha) | length(g.alpha)~=1)
   error('timef(): Value of g.alpha must be a number.\n');
elseif (round(g.naccu*g.alpha) < 2)
   fprintf('Value of g.alpha is out of the normal range [%g,0.5]\n',2/g.naccu);
   g.naccu = round(2/g.alpha);
   fprintf('  Increasing the number of bootstrap iterations to %d\n',g.naccu);
end
if g.alpha>0.5 | g.alpha<=0
   error('Value of g.alpha is out of the allowed range (0.00,0.5).');
end
if ~isnan(g.alpha)
   if g.baseboot > 0
      fprintf('Bootstrap analysis will use data in baseline (pre-0) subwindows only.\n')
   else
      fprintf('Bootstrap analysis will use data in all subwindows.\n')
   end
end
switch g.angleunit
   case { 'rad', 'ms', 'deg' },;
   otherwise error('Angleunit must be either ''rad'', ''deg'', or ''ms''');
end;
switch g.type
   case { 'coher', 'phasecoher' 'phasecoher2' },;
   otherwise error('Type must be either ''coher'' or ''phasecoher''');
end;
switch g.boottype
   case { 'times' 'timestrials' 'trials'},;
   otherwise error('Boot type must be either ''times'', ''trials'' or ''timestrials''');
end;
if (~isnumeric(g.shuffle))
   error('Shuffle argument type must be numeric');
end;
switch g.memory
   case { 'low', 'high' },;
   otherwise error('memory must be either ''low'' or ''high''');
end;
if strcmp(g.memory, 'low') & ~strcmp(g.boottype, 'times')
   error(['Bootstrap type ''' g.boottype ''' cannot be used in low memory mode']);
end;

switch g.compute
   case { 'matlab', 'c' },;
   otherwise error('compute must be either ''matlab'' or ''c''');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare 2 conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if iscell(X)
	if length(X) ~= 2 | length(Y) ~= 2
		error('crossf: to compare conditions, X and Y input must be 2-elements cell arrays');
	end;
	if ~strcmp(g.boottype, 'times')
		disp('crossf warning: The significance bootstrap type is irrelevant when comparing conditions');
	end;
	for index = 1:2:length(vararginori)
		if index<=length(vararginori) % needed: if elemenets are deleted
			%if strcmp(vararginori{index}, 'alpha'), vararginori(index:index+1) = [];
			if strcmp(vararginori{index}, 'title'), vararginori(index:index+1) = [];
			end;
		end;
	end;
	if iscell(g.title)
		if length(g.title) <= 2,
			g.title{3} = 'Condition 2 - condition 1';
		end;
	else
		g.title = { 'Condition 1', 'Condition 2', 'Condition 2 - condition 1' };
	end;

	fprintf('Running crossf on condition 1 *********************\n');
	fprintf('Note: If an out-of-memory error occurs, try reducing the\n');
	fprintf('      number of time points or number of frequencies\n');
	if ~strcmp(g.type, 'coher')
	   fprintf('Note: Type ''coher'' takes 3 times as much memory as other options!)\n');
        end
	figure;
	subplot(1,3,1); title(g.title{1});
	if ~strcmp(g.type, 'coher')
		[R1,mbase,times,freqs,Rbootout1,Rangle1, savecoher1] = crossf(X{1}, Y{1}, ...
								frame, tlimits, Fs, varwin, 'savecoher', 1, 'title', ' ',vararginori{:});
	else
		[R1,mbase,times,freqs,Rbootout1,Rangle1, savecoher1, Tfx1, Tfy1] = crossf(X{1}, Y{1}, ...
								frame, tlimits, Fs, varwin, 'savecoher', 1,'title', ' ',vararginori{:});
	end;
	R1 = R1.*exp(j*Rangle1); % output Rangle is in radians

	% Asking user for memory limitations
	% if ~strcmp(g.noinput, 'yes')
	%	  tmp = whos('Tfx1');
	%	  fprintf('This function will require an additional %d bytes, do you wish\n', ...
    %     tmp.bytes*6+size(savecoher1,1)*size(savecoher1,2)*g.naccu*8);
	%	  res = input('to continue (y/n) (use the ''noinput'' option to disable this message):', 's');
	%  	if res == 'n', return; end;
	% end;

	fprintf('\nRunning crossf on condition 2 *********************\n');
	subplot(1,3,2); title(g.title{2});
	if ~strcmp(g.type, 'coher')
		  [R2,mbase,times,freqs,Rbootout2,Rangle2, savecoher2] = crossf(X{2}, Y{2}, ...
								frame, tlimits, Fs, varwin,'savecoher', 1, 'title', ' ',vararginori{:});
	else
		  [R2,mbase,times,freqs,Rbootout2,Rangle2, savecoher2, Tfx2, Tfy2] = crossf(X{2}, Y{2}, ...
								frame, tlimits, Fs, varwin,'savecoher', 1, 'title', ' ',vararginori{:});
	end;
	R2 = R2.*exp(j*Rangle2); % output Rangle is in radians

	subplot(1,3,3); title(g.title{3});
	if isnan(g.alpha)
		plotall(R2-R1, [], [], times, freqs, mbase,  find(freqs <= g.maxfreq), g);
	else
		% accumulate coherence images (all arrays [nb_points * timesout * trials])
		% ---------------------------
		allsavedcoher = zeros(size(savecoher1,1), ...
                          size(savecoher1,2), ...
                          size(savecoher1,3)+size(savecoher2,3));
		allsavedcoher(:,:,1:size(savecoher1,3))     = savecoher1;
		allsavedcoher(:,:,size(savecoher1,3)+1:end) = savecoher2;
		clear savecoher1 savecoher2;

		if strcmp(g.type, 'coher')
			alltfx = zeros(size(Tfx1,1), size(Tfx2,2), size(Tfx1,3)+size(Tfx2,3));
			alltfx(:,:,1:size(Tfx1,3))     = Tfx1;
			alltfx(:,:,size(Tfx1,3)+1:end) = Tfx2;
			clear Tfx1 Tfx2;

			alltfy = zeros(size(Tfy1,1), size(Tfy2,2), size(Tfy1,3)+size(Tfy2,3));
			alltfy(:,:,1:size(Tfy1,3))   = Tfy1;
			alltfy(:,:,size(Tfy1,3)+1:end) = Tfy2;
			clear Tfy1 Tfy2;
		end;

		coherimages = zeros(size(allsavedcoher,1), size(allsavedcoher,2), g.naccu);
		cond1trials = length(X{1})/g.frame;
		cond2trials = length(X{2})/g.frame;
		alltrials = [1:cond1trials+cond2trials];
		fprintf('Accumulating bootstrap:');

		% preprocess data
		% ---------------
		switch g.type
		 case 'coher', % take the square of alltfx and alltfy
		  alltfx = alltfx.^2;
		  alltfy = alltfy.^2;
		 case 'phasecoher', % normalize
		  allsavedcoher = allsavedcoher ./ abs(allsavedcoher);
		 case 'phasecoher2', % don't do anything
		end;

		if strcmp(g.type, 'coher')
			[coherdiff coher1 coher2] = coher2conddiff( allsavedcoher, alltrials, ...
                                                        cond1trials, g.type, alltfx, alltfy);
		else
			[coherdiff coher1 coher2] = coher2conddiff( allsavedcoher, alltrials, ...
                                                        cond1trials, g.type);
		end;
		%figure; g.alpha = NaN; & to check that the new images are the same as the original
		%subplot(1,3,1); plotall(coher1, [], [], times, freqs, mbase, find(freqs <= g.maxfreq), g);
		%subplot(1,3,2); plotall(coher2, [], [], times, freqs, mbase, find(freqs <= g.maxfreq), g);
		%return;

		for index=1:g.naccu
			if rem(index,10) == 0,  fprintf(' %d',index); end
			if rem(index,120) == 0, fprintf('\n'); end

			if strcmp(g.type, 'coher')
				coherimages(:,:,index) = coher2conddiff( allsavedcoher, shuffle(alltrials), ...
                                                        cond1trials, g.type, alltfx, alltfy);
			else
				coherimages(:,:,index) = coher2conddiff( allsavedcoher, shuffle(alltrials), ...
                                                        cond1trials, g.type);
			end;
		end;
		fprintf('\n');

		% create articially a Bootstrap object to compute significance
		Boot = bootinit( [], size(allsavedcoher,1), g.timesout, g.naccu, 0, g.baseboot, ...
														'noboottype', g.alpha, g.rboot);
		Boot.Coherboot.R = coherimages;
		Boot = bootcomppost(Boot, [], [], []);
		g.title = '';
		plotall(coherdiff, Boot.Coherboot.R, Boot.Rsignif, times, freqs, mbase, ...
														find(freqs <= g.maxfreq), g);
	end;
	return; % ********************************** END PROCESSING TWO CONDITIONS
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shuffle trials if necessary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if g.shuffle ~= 0
   fprintf('x and y data trials being shuffled %d times\n',g.shuffle);
   XX = reshape(X, 1, frame, length(X)/g.frame);
   YY = Y;
   X = [];
   Y = [];
   for index = 1:g.shuffle
      XX = shuffle(XX,3);
      X = [X XX(:,:)];
      Y = [Y YY];
   end;
end;

% detrend over epochs (trials) if requested
% -----------------------------------------
switch g.detrep
case 'on'
   X = reshape(X, g.frame, length(X)/g.frame);
   X = X - mean(X,2)*ones(1, length(X(:))/g.frame);
   Y = reshape(Y, g.frame, length(Y)/g.frame);
   Y = Y - mean(Y,2)*ones(1, length(Y(:))/g.frame);
end;

% time limits
wintime = 500*g.winsize/g.srate;
times = [g.tlimits(1)+wintime:(g.tlimits(2)-g.tlimits(1)-2*wintime)/(g.timesout-1):g.tlimits(2)-wintime];

%%%%%%%%%%
% baseline
%%%%%%%%%%
if ~isnan(g.baseline)
   baseln = find(times < g.baseline); % subtract means of pre-0 (centered) windows
   if isempty(baseln)
      baseln = 1:length(times); % use all times as baseline
      disp('Bootstrap baseline empty, using the whole epoch.');
   end;
   baselength = length(baseln);
else
   baseln = 1:length(times); % use all times as baseline
   baselength = length(times); % used for bootstrap
end;

%%%%%%%%%%%%%%%%%%%%
% Initialize objects
%%%%%%%%%%%%%%%%%%%%
tmpsaveall = (~isnan(g.alpha) & isnan(g.rboot) & strcmp(g.memory, 'high')) ...
                 | (strcmp(g.subitc, 'on') & strcmp(g.memory, 'high'));
trials = length(X)/g.frame;
if ~strcmp(lower(g.compute), 'c')
   Tfx = tfinit(X, g.timesout, g.winsize, g.cycles, g.frame, g.padratio, g.detret, ...
									g.srate, g.maxfreq, g.subitc, g.type, g.cyclesfact, tmpsaveall);
   Tfy = tfinit(Y, g.timesout, g.winsize, g.cycles, g.frame, g.padratio, g.detret, ...
									g.srate, g.maxfreq, g.subitc, g.type, g.cyclesfact, tmpsaveall);
   Coher     = coherinit(Tfx.nb_points, trials, g.timesout, g.type);
   Coherboot = coherinit(Tfx.nb_points, trials, g.naccu   , g.type);
   Boot      = bootinit( Coherboot, Tfx.nb_points, g.timesout, g.naccu, baselength, ...
									g.baseboot, g.boottype, g.alpha, g.rboot);
   freqs = Tfx.freqs;
   dispf = find(freqs <= g.maxfreq);
   freqs = freqs(dispf);
else
   freqs = g.srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
end;
dispf     = find(Tfx.freqs <= g.maxfreq);

%-------------
% Reserve space
%-------------
% R  = zeros(tfx.nb_points,g.timesout);       % mean coherence
% RR = repmat(nan,tfx.nb_points,g.timesout); % initialize with nans
% Rboot = zeros(tfx.nb_points,g.naccu);  % summed bootstrap coher
% switch g.type
% case 'coher',
%    cumulXY = zeros(tfx.nb_points,g.timesout);
%    cumulXYboot = zeros(tfx.nb_points,g.naccu);
% end;
% if g.bootsub > 0
%    Rboottrial = zeros(tfx.nb_points, g.timesout, g.bootsub); % summed bootstrap coher
%    cumulXYboottrial = zeros(tfx.nb_points, g.timesout, g.bootsub);
% end;
% if ~isnan(g.alpha) & isnan(g.rboot)
%    tf.tmpalltimes = repmat(nan,tfx.nb_points,g.timesout);
% end

% --------------------
% Display text to user
% --------------------
fprintf('\nComputing Event-Related ');
switch g.type
    case 'phasecoher',  fprintf('Phase Coherence (ITC) images for %d trials.\n',length(X)/g.frame);
    case 'phasecoher2', fprintf('Phase Coherence 2 (ITC) images for %d trials.\n',length(X)/g.frame);
    case 'coher',       fprintf('Linear Coherence (ITC) images for %d trials.\n',length(X)/g.frame);
end;
fprintf('The trial latency range is from %4.5g ms before to %4.5g ms after\n     the time-locking event.\n', g.tlimits(1),g.tlimits(2));
fprintf('The frequency range displayed will be %g-%g Hz.\n',min(freqs),g.maxfreq);
if ~isnan(g.baseline)
   if length(baseln) == length(times)
      fprintf('Using the full trial latency range as baseline.\n');
   else
      fprintf('Using trial latencies from %4.5g ms to %4.5g ms as baseline.\n', g.tlimits,g.baseline);
   end;
else
   fprintf('No baseline time range was specified.\n');
end;
if g.cycles==0
   fprintf('The data window size will be %d samples (%g ms).\n',g.winsize,2*wintime);
   fprintf('The FFT length will be %d samples\n',g.winsize*g.padratio);
else
   fprintf('The window size will be %2.3g cycles.\n',g.cycles);
   fprintf('The maximum window size will be %d samples (%g ms).\n',g.winsize,2*wintime);
end
fprintf('The window will be applied %d times\n',g.timesout);
fprintf('     with an average step size of %2.2g samples (%2.4g ms).\n', Tfx.stp,1000*Tfx.stp/g.srate);
fprintf('Results will be oversampled %d times.\n',g.padratio);
if ~isnan(g.alpha)
   fprintf('Bootstrap confidence limits will be computed based on alpha = %g\n', g.alpha);
else
   fprintf('Bootstrap confidence limits will NOT be computed.\n');
end
switch g.plotphase
case 'on',
    if strcmp(g.angleunit,'deg')
       fprintf(['Coherence angles will be imaged in degrees.\n']);
    elseif strcmp(g.angleunit,'rad')
       fprintf(['Coherence angles will be imaged in radians.\n']);
    elseif strcmp(g.angleunit,'ms')
       fprintf(['Coherence angles will be imaged in ms.\n']);
    end
end;
fprintf('\nProcessing trial (of %d): ',trials);

% firstboot = 1;
% Rn=zeros(trials,g.timesout);
% X = X(:)'; % make X and Y column vectors
% Y = Y(:)';
% tfy = tfx;

if strcmp(lower(g.compute), 'c')
   % C PART
   filename = [ 'tmpcrossf' num2str(round(rand(1)*1000)) ];
   f = fopen([ filename '.in'], 'w');
   fwrite(f, tmpsaveall, 'int32');
   fwrite(f, g.detret, 'int32');
   fwrite(f, g.srate, 'int32');
   fwrite(f, g.maxfreq, 'int32');
   fwrite(f, g.padratio, 'int32');
   fwrite(f, g.cycles, 'int32');
   fwrite(f, g.winsize, 'int32');
   fwrite(f, g.timesout, 'int32');
   fwrite(f, g.subitc, 'int32');
   fwrite(f, g.type, 'int32');
   fwrite(f, trials, 'int32');
   fwrite(f, g.naccu, 'int32');
   fwrite(f, length(X), 'int32');
   fwrite(f, X, 'double');
   fwrite(f, Y, 'double');
   fclose(f);

   command = [ '!cppcrosff ' filename '.in ' filename '.out' ];
   eval(command);

   f = fopen([ filename '.out'], 'r');
   size1 = fread(f, 'int32', 1);
   size2 = fread(f, 'int32', 1);
   Rreal = fread(f, 'double', [size1 size2]);
   Rimg  = fread(f, 'double', [size1 size2]);
   Coher.R = Rreal + j*Rimg;
   Boot.Coherboot.R = [];
   Boot.Rsignif = [];
else
   % ------------------------
   % MATLAB PART
   % compute ITC if necessary
   % ------------------------
   if strcmp(g.subitc, 'on')
      for t=1:trials
         if rem(t,10) == 0,  fprintf(' %d',t); end
         if rem(t,120) == 0, fprintf('\n'); end
         Tfx = tfitc( Tfx, t, 1:g.timesout);
         Tfy = tfitc( Tfy, t, 1:g.timesout);
      end;
	  fprintf('\n');
      Tfx = tfitcpost( Tfx, trials);
      Tfy = tfitcpost( Tfy, trials);
   end;

   % ---------
   % Main loop
   % ---------
   if g.savecoher,
	   trialcoher = zeros(Tfx.nb_points, g.timesout, trials);
   else
	   trialcoher = [];
   end;
   for t=1:trials
      if rem(t,10) == 0,  fprintf(' %d',t); end
      if rem(t,120) == 0, fprintf('\n'); end

      Tfx = tfcomp( Tfx, t, 1:g.timesout);
      Tfy = tfcomp( Tfy, t, 1:g.timesout);
	  if g.savecoher
		  [Coher trialcoher(:,:,t)] = cohercomp( Coher, Tfx.tmpalltimes, ...
                                             Tfy.tmpalltimes, t, 1:g.timesout);
      else
		  Coher = cohercomp( Coher, Tfx.tmpalltimes, Tfy.tmpalltimes, t, 1:g.timesout);
	  end;

      Boot = bootcomp( Boot, Coher.Rn(t,:), Tfx.tmpalltimes, Tfy.tmpalltimes);
   end % t = trial
   [Boot Rbootout] = bootcomppost(Boot, Coher.Rn, Tfx.tmpall, Tfy.tmpall);
      % Note that the bootstrap thresholding is actually performed
      %      in the display subfunction plotall()

   Coher  = cohercomppost(Coher, trials);
end;

% ----------------------------------
% If coherence, perform the division
% ----------------------------------
% switch g.type
% case 'coher',
%   R = R ./ cumulXY;
%   if ~isnan(g.alpha) & isnan(g.rboot)
%      Rboot = Rboot ./ cumulXYboot;
%   end;
%   if g.bootsub > 0
%      Rboottrial = Rboottrial ./ cumulXYboottrial;
%   end;
% case 'phasecoher',
%   Rn = sum(Rn, 1);
%   R = R ./ (ones(size(R,1),1)*Rn);               % coherence magnitude
%   if ~isnan(g.alpha) & isnan(g.rboot)
%      Rboot = Rboot / trials;
%   end;
%   if g.bootsub > 0
%      Rboottrial = Rboottrial / trials;
%   end;
% end;

% ----------------
% Compute baseline
% ----------------
mbase = mean(abs(Coher.R(:,baseln)'));     % mean baseline coherence magnitude

% ---------------
% Plot everything
% ---------------
plotall(Coher.R, Boot.Coherboot.R, Boot.Rsignif, times, freqs, mbase, dispf, g);

% --------------------------------------
% Convert output Rangle to degrees or ms - Disabled to keep original default: radians output
% --------------------------------------
% Rangle = angle(Coher.R); % returns radians
% if strcmp(g.angleunit,'ms')  % convert to ms
%    Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times));
% elseif strcmp(g.angleunit,'deg')  % convert to deg
%    Rangle = Rangle*180/pi; % convert to degrees
% else % angleunit is 'rad'
%    % Rangle = Rangle;
% end
% Rangle(find(Rraw==0)) = 0; % mask for significance - set angle at non-signif coher points to 0

R = abs(Coher.R);
Rsignif = Boot.Rsignif;
Tfx = permute(Tfx.tmpall, [3 2 1]); % from [trials timesout nb_points]
%                                     to   [nb_points timesout trials]
Tfy = permute(Tfy.tmpall, [3 2 1]);

return; % end crossf() *************************************************

%
% crossf() plotting functions
% ----------------------------------------------------------------------
function plotall(R, Rboot, Rsignif, times, freqs, mbase, dispf, g)

switch lower(g.plotphase)
case 'on',
   switch lower(g.plotamp),
   case 'on', ordinate1 = 0.67; ordinate2 = 0.1; height = 0.33; g.plot = 1;
   case 'off', ordinate2 = 0.1; height = 0.9; g.plot = 1;
   end;
case 'off', ordinate1 = 0.1; height = 0.9;
   switch lower(g.plotamp),
   case 'on', ordinate1 = 0.1; height = 0.9;  g.plot = 1;
   case 'off', g.plot = 0;
   end;
end;

%
% Compute cross-spectral angles
% -----------------------------
Rangle = angle(R); % returns radians

%
% Optionally convert Rangle to degrees or ms
% ------------------------------------------
if strcmp(g.angleunit,'ms')  % convert to ms
   Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times));
   maxangle = max(max(abs(Rangle)));
elseif strcmp(g.angleunit,'deg')  % convert to degrees
   Rangle = Rangle*180/pi; % convert to degrees
   maxangle = 180; % use full-cycle plotting
else
   maxangle = pi;  % radians
end

R = abs(R);

% if ~isnan(g.baseline)
% 	R = R - repmat(mbase',[1 g.timesout]); % remove baseline mean
% end;

Rraw = R; % raw coherence (e.g., coherency) magnitude values output

if g.plot
   fprintf('\nNow plotting...\n');
   set(gcf,'DefaultAxesFontSize',g.AXES_FONT)
   colormap(jet(256));

   pos = get(gca,'position'); % plot relative to current axes
   q = [pos(1) pos(2) 0 0];
   s = [pos(3) pos(4) pos(3) pos(4)];
   axis('off')
end;

switch lower(g.plotamp)
case 'on'
   %
   % Image the coherence [% perturbations]
   %
   RR = R;
   if ~isnan(g.alpha) % zero out (and 'green out') nonsignif. R values
      RR(find(RR < repmat(Rboot(:),[1 g.timesout]))) = 0;
      Rraw(find(repmat(Rsignif(:),[1,size(Rraw,2)])>=Rraw))=0;
   end

   if g.cmax == 0
      coh_caxis = max(max(R(dispf,:)))*[-1 1];
   else
      coh_caxis = g.cmax*[-1 1];
   end

   h(6) = axes('Units','Normalized', 'Position',[.1 ordinate1 .8 height].*s+q);

   map=hsv(300); % install circular color map - green=0, yellow, orng, red, violet = max
   %                                         cyan, blue, violet = min
   map = flipud([map(251:end,:);map(1:250,:)]);
   map(151,:) = map(151,:)*0.9; % tone down the (0=) green!
   colormap(map);

   imagesc(times,freqs(dispf),RR(dispf,:),coh_caxis); % plot the coherence image
   if ~isempty(g.maxamp)
	   caxis([-g.maxamp g.maxamp]);
   end;
   tmpscale = caxis;

   hold on
   plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth)
   for i=1:length(g.marktimes)
      plot([g.marktimes(i) g.marktimes(i)],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
   end;
   hold off
   set(h(6),'YTickLabel',[],'YTick',[])
   set(h(6),'XTickLabel',[],'XTick',[])
   %title('Event-Related Coherence')

   h(8) = axes('Position',[.95 ordinate1 .05 height].*s+q);
   cbar(h(8),151:300, [0 tmpscale(2)]); % use only positive colors (gyorv)

   %
   % Plot delta-mean min and max coherence at each time point on bottom of image
   %
   h(10) = axes('Units','Normalized','Position',[.1 ordinate1-0.1 .8 .1].*s+q);
                                                            % plot marginal means below
   Emax = max(R(dispf,:)); % mean coherence at each time point
   Emin = min(R(dispf,:)); % mean coherence at each time point
   plot(times,Emin, times, Emax, 'LineWidth',g.linewidth); hold on;
   plot([times(1) times(length(times))],[0 0],'LineWidth',0.7);
   plot([0 0],[-500 500],'--m','LineWidth',g.linewidth);
   for i=1:length(g.marktimes)
       plot([g.marktimes(i) g.marktimes(i)],[-500 500],'--m','LineWidth',g.linewidth);
   end;
   if ~isnan(g.alpha) & strcmp(g.boottype, 'trials')
       % plot bootstrap significance limits (base mean +/-)
      plot(times,mean(Rboot(dispf,:)),'g','LineWidth',g.linewidth); hold on;
      plot(times,mean(Rsignif(dispf,:)),'k:','LineWidth',g.linewidth);
      axis([min(times) max(times) 0 max([Emax(:)' Rsignif(:)'])*1.2])
   else
      axis([min(times) max(times) 0 max(Emax)*1.2])
   end;
   tick = get(h(10),'YTick');
   set(h(10),'YTick',[tick(1) ; tick(length(tick))])
   set(h(10),'YAxisLocation','right')
   xlabel('Time (ms)')
   ylabel('coh.')

   %
   % Plot mean baseline coherence at each freq on left side of image
   %
   h(11) = axes('Units','Normalized','Position',[0 ordinate1 .1 height].*s+q);
                                                            % plot mean spectrum
   E = abs(mbase(dispf)); % baseline mean coherence at each frequency
   plot(freqs(dispf),E,'LineWidth',g.linewidth); % plot mbase
   if ~isnan(g.alpha) % plot bootstrap significance limits (base mean +/-)
      hold on
      % plot(freqs(dispf),Rboot(:,dispf)+[E;E],'g','LineWidth',g.linewidth);
      plot(freqs(dispf),mean(Rboot  (dispf,:),2),'g','LineWidth',g.linewidth);
      plot(freqs(dispf),mean(Rsignif(dispf,:),2),'k:','LineWidth',g.linewidth);
      axis([freqs(1) freqs(max(dispf)) 0 max([E Rsignif(:)'])*1.2]);
   else             % plot marginal mean coherence only
      if ~isnan(max(E))
         axis([freqs(1) freqs(max(dispf)) 0 max(E)*1.2]);
      end;
   end

   tick = get(h(11),'YTick');
   set(h(11),'YTick',[tick(1) ; tick(length(tick))])
   set(h(11),'View',[90 90])
   xlabel('Freq. (Hz)')
   ylabel('coh.')
end;

switch lower(g.plotphase)
case 'on'
   %
   % Plot coherence phase lags in bottom panel
   %
   h(13) = axes('Units','Normalized','Position',[.1 ordinate2 .8 height].*s+q);
   Rangle(find(Rraw==0)) = 0; % when plotting, mask for significance
                              % = set angle at non-signif coher points to 0

   imagesc(times,freqs(dispf),Rangle(dispf,:),[-maxangle maxangle]); % plot the
   hold on                                             % coherence phase angles
   plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth); % zero-time line
   for i=1:length(g.marktimes)
      plot([g.marktimes(i) g.marktimes(i)],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
   end;

   ylabel('Freq. (Hz)')
   xlabel('Time (ms)')

   h(14)=axes('Position',[.95 ordinate2 .05 height].*s+q);
   cbar(h(14),0,[-maxangle maxangle]); % two-sided colorbar
end

if g.plot
	try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
    if (length(g.title) > 0) % plot title
        if h(6) ~= 0, axes(h(6)); else axes(h(13)); end;
        %h = subplot('Position',[0 0  1 1].*s+q, 'Visible','Off');
        %h(13) = text(-.05,1.01,g.title);
        h(13) = title(g.title);
        %set(h(13),'VerticalAlignment','bottom')
        %set(h(13),'HorizontalAlignment','left')
        set(h(13),'FontSize',g.TITLE_FONT);
    end
   %
   %%%%%%%%%%%%%%% plot topoplot() %%%%%%%%%%%%%%%%%%%%%%%
   %
   if (~isempty(g.topovec))
      h(15) = subplot('Position',[-.1 .43 .2 .14].*s+q);
      if size(g.topovec,2) <= 2
         topoplot(g.topovec(1),g.elocs,'electrodes','off', ...
            'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);
      else
         topoplot(g.topovec(1,:),g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);
      end;
      axis('square')

      h(16) = subplot('Position',[.9 .43 .2 .14].*s+q);
      if size(g.topovec,2) <= 2
         topoplot(g.topovec(2),g.elocs,'electrodes','off', ...
            'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);
      else
         topoplot(g.topovec(2,:),g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);
      end;
      axis('square')
   end

   axcopy(gcf);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  TIME FREQUENCY   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for time freq initialisation
% -------------------------------------
function Tf = tfinit(X, timesout, winsize, ...
   cycles, frame, padratio, detret, srate, maxfreq, subitc, type, cyclesfact, saveall);
Tf.X         = X(:)'; % make X column vectors
Tf.winsize   = winsize;
Tf.cycles    = cycles;
Tf.frame     = frame;
Tf.padratio  = padratio;
Tf.detret    = detret;
Tf.stp       = (frame-winsize)/(timesout-1);
Tf.subitc    = subitc; % for ITC
Tf.type      = type; % for ITC
Tf.saveall   = saveall;
if (Tf.cycles == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%
   % Tf.freqs = srate/winsize*[1:2/padratio:winsize]/2; % incorect for padratio > 2
   Tf.freqs = linspace(0, srate/2, length([1:2/padratio:winsize])+1);
   Tf.freqs = Tf.freqs(2:end);
   Tf.win   = hanning(winsize);
   Tf.nb_points = padratio*winsize/2;
else % %%%%%%%%%%%%%%%%%% Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   Tf.freqs = srate*cycles/winsize*[2:2/padratio:winsize]/2;
   Tf.win = dftfilt(winsize,maxfreq/srate,cycles,padratio,cyclesfact);
   Tf.nb_points = size(Tf.win,2);
end;
Tf.tmpalltimes = zeros(Tf.nb_points, timesout);
trials = length(X)/frame;
if saveall
   Tf.tmpall     = repmat(nan,[trials timesout Tf.nb_points]);
else
   Tf.tmpall = [];
end
Tf.tmpallbool = zeros(trials,timesout);
Tf.ITCdone = 0;
if Tf.subitc
   Tf.ITC  = zeros(Tf.nb_points, timesout);
   switch Tf.type,
	   case { 'coher' 'phasecoher2' }
		Tf.ITCcumul  = zeros(Tf.nb_points, timesout);
   end;
end;

% function for itc
% ----------------
function [Tf, itcvals] = tfitc(Tf, trials, times);
Tf = tfcomp(Tf, trials, times);
switch Tf.type
   case 'coher',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes; % complex coher.
      Tf.ITCcumul(:,times) = Tf.ITCcumul(:,times)+abs(Tf.tmpalltimes).^2;
   case 'phasecoher2',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes; % complex coher.
      Tf.ITCcumul(:,times) = Tf.ITCcumul(:,times)+abs(Tf.tmpalltimes);
   case 'phasecoher',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes ./ abs(Tf.tmpalltimes);
                                                            % complex coher.
end % ~any(isnan())
return;

function [Tf, itcvals] = tfitcpost(Tf, trials);
switch Tf.type
   case 'coher',       Tf.ITC = Tf.ITC ./ sqrt(trials * Tf.ITCcumul);
   case 'phasecoher2', Tf.ITC = Tf.ITC ./ Tf.ITCcumul;
   case 'phasecoher',  Tf.ITC = Tf.ITC / trials; % complex coher.
end % ~any(isnan())

if Tf.saveall
  Tf.ITC = transpose(Tf.ITC); % do not use ' otherwise conjugate

	%imagesc(abs(Tf.ITC)); colorbar; figure;
	%squeeze(Tf.tmpall(1,1,1:Tf.nb_points))
	%squeeze(Tf.ITC   (1,1,1:Tf.nb_points))
	%Tf.ITC = shiftdim(Tf.ITC, -1);

	Tf.ITC = repmat(shiftdim(Tf.ITC, -1), [trials 1 1]);
	Tf.tmpall = (Tf.tmpall - abs(Tf.tmpall) .* Tf.ITC) ./ abs(Tf.tmpall);

  %	for index = 1:trials
  %		imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); drawnow; figure;
  %		Tf.tmpall(index,:,:) = (Tf.tmpall(index,:,:) - Tf.tmpall(index,:,:) .* Tf.ITC)./Tf.tmpall(index,:,:);
  %		imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); drawnow;
  %		subplot(10,10, index); imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); caxis([0 1]); drawnow;
  %	end;
  %	squeeze(Tf.tmpall(1,1,1:Tf.nb_points))
  %	figure; axcopy;

end;
Tf.ITCdone = 1;
return;

% function for time freq decomposition
% ------------------------------------
function [Tf, tmpX] = tfcomp(Tf, trials, times);
% tf is an structure containing all the information about the decomposition
for trial = trials
   for index = times
      if ~Tf.tmpallbool(trial, index) % already computed
         tmpX = Tf.X([1:Tf.winsize]+floor((index-1)*Tf.stp)+(trial-1)*Tf.frame);

         if ~any(isnan(tmpX)) % perform the decomposition
            tmpX = tmpX - mean(tmpX);
            switch Tf.detret, case 'on',
               tmpX = detrend(tmpX);
            end;

            if Tf.cycles == 0 % use FFTs
               tmpX = Tf.win .* tmpX(:);
               tmpX = fft(tmpX,Tf.padratio*Tf.winsize);
               tmpX = tmpX(2:Tf.padratio*Tf.winsize/2+1);
            else
               tmpX = transpose(Tf.win) * tmpX(:);
            end
         else
            tmpX = NaN;
         end;
         if Tf.ITCdone
            tmpX = (tmpX - abs(tmpX) .* Tf.ITC(:,index)) ./ abs(tmpX);
         end;
         Tf.tmpalltimes(:,index) = tmpX;
         if Tf.saveall
            Tf.tmpall(trial, index,:) = tmpX;
            Tf.tmpallbool(trial, index) = 1;
         end
	  else
		  Tf.tmpalltimes(:,index) = Tf.tmpall(trial, index,:);
     end;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    COHERENCE    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for coherence initialisation
% -------------------------------------
function Coher = coherinit(nb_points, trials, timesout, type);
Coher.R  = zeros(nb_points,timesout);       % mean coherence
% Coher.RR = repmat(nan,nb_points,timesout); % initialize with nans
Coher.type = type;
Coher.Rn=zeros(trials,timesout);
switch type
 case 'coher',
  Coher.cumulX = zeros(nb_points,timesout);
  Coher.cumulY = zeros(nb_points,timesout);
 case 'phasecoher2',
  Coher.cumul  = zeros(nb_points,timesout);
end;

% function for coherence calculation
% -------------------------------------
% function Coher = cohercomparray(Coher, tmpX, tmpY, trial);
% switch Coher.type
%   case 'coher',
%      Coher.R = Coher.R + tmpX.*conj(tmpY); % complex coher.
%      Coher.cumulXY = Coher.cumulXY + abs(tmpX).*abs(tmpY);
%   case 'phasecoher',
%      Coher.R = Coher.R + tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY)); % complex coher.
%      Coher.Rn(trial,:) = 1;
% end % ~any(isnan())

function [Coher,tmptrialcoh] = cohercomp(Coher, tmpX, tmpY, trial, time);
tmptrialcoh = tmpX.*conj(tmpY);
switch Coher.type
   case 'coher',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh; % complex coher.
      Coher.cumulX(:,time) = Coher.cumulX(:,time) + abs(tmpX).^2;
      Coher.cumulY(:,time) = Coher.cumulY(:,time) + abs(tmpY).^2;
 case 'phasecoher2',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh; % complex coher.
      Coher.cumul(:,time) = Coher.cumul(:,time) + abs(tmptrialcoh);
   case 'phasecoher',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh ./ abs(tmptrialcoh); % complex coher.
	  %figure; imagesc(abs(tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY))));
      Coher.Rn(trial,time) = Coher.Rn(trial,time)+1;
end % ~any(isnan())

% function for post coherence calculation
% ---------------------------------------
function Coher = cohercomppost(Coher, trials);
switch Coher.type
 case 'coher',
   Coher.R = Coher.R ./ sqrt(Coher.cumulX) ./ sqrt(Coher.cumulY);
 case 'phasecoher2',
   Coher.R = Coher.R ./ Coher.cumul;
 case 'phasecoher',
   Coher.Rn = sum(Coher.Rn, 1);
   Coher.R  = Coher.R ./ (ones(size(Coher.R,1),1)*Coher.Rn); % coherence magnitude
end;

% function for 2 conditions coherence calculation
% -----------------------------------------------
function [coherimage, coherimage1, coherimage2] = coher2conddiff( allsavedcoher, alltrials, cond1trials, type, tfx, tfy);
	t1s = alltrials(1:cond1trials);
	t2s = alltrials(cond1trials+1:end);
	switch type
	 case 'coher',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) ./ sqrt(sum(tfx(:,:,t1s),3)) ./ sqrt(sum(tfy(:,:,t1s),3));
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) ./ sqrt(sum(tfx(:,:,t2s),3)) ./ sqrt(sum(tfy(:,:,t1s),3));
	 case 'phasecoher2',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) ./ sum(abs(allsavedcoher(:,:,t1s)),3);
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) ./ sum(abs(allsavedcoher(:,:,t2s)),3);
	 case 'phasecoher',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) / cond1trials;
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) / (size(allsavedcoher,3)-cond1trials);
	end;
	coherimage = coherimage2 - coherimage1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BOOTSTRAP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for bootstrap initialisation
% -------------------------------------
function Boot = bootinit(Coherboot, nb_points, timesout, naccu, baselength, baseboot, boottype, alpha, rboot);
Boot.Rboot       = zeros(nb_points,naccu);  % summed bootstrap coher
Boot.boottype    = boottype;
Boot.baselength  = baselength;
Boot.baseboot    = baseboot;
Boot.Coherboot   = Coherboot;
Boot.naccu       = naccu;
Boot.alpha       = alpha;
Boot.rboot       = rboot;

% function for bootstrap computation
% ----------------------------------
function Boot = bootcomp(Boot, Rn, tmpalltimesx, tmpalltimesy);
if ~isnan(Boot.alpha) & isnan(Boot.rboot)
   if strcmp(Boot.boottype, 'times') % get g.naccu bootstrap estimates for each trial
      goodbasewins = find(Rn==1);
      if Boot.baseboot % use baseline windows only
         goodbasewins = find(goodbasewins<=Boot.baselength);
      end
      ngdbasewins = length(goodbasewins);
      j=1;
      tmpsX = zeros(size(tmpalltimesx,1), Boot.naccu);
      tmpsY = zeros(size(tmpalltimesx,1), Boot.naccu);
      if ngdbasewins > 1
         while j<=Boot.naccu
            s = ceil(rand([1 2])*ngdbasewins); % random ints [1,g.timesout]
            s = goodbasewins(s);
            if ~any(isnan(tmpalltimesx(:,s(1)))) & ~any(isnan(tmpalltimesy(:,s(2))))
               tmpsX(:,j) = tmpalltimesx(:,s(1));
               tmpsY(:,j) = tmpalltimesy(:,s(2));
               j = j+1;
            end
         end
         Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
      end;
   end
end;

% handle other trial bootstrap types
% ----------------------------------
function [Boot, Rbootout] = bootcomppost(Boot, allRn, alltmpsX, alltmpsY);
trials    = size(alltmpsX, 1);
times     = size(alltmpsX, 2);
nb_points = size(alltmpsX, 3);
if ~isnan(Boot.alpha) & isnan(Boot.rboot)
   if strcmp(Boot.boottype, 'trials') % get g.naccu bootstrap estimates for each trial
      fprintf('\nProcessing trial bootstrap (of %d):',times(end));
      tmpsX = zeros(size(alltmpsX,3), Boot.naccu);
      tmpsY = zeros(size(alltmpsY,3), Boot.naccu );
      Boot.fullcoherboot = zeros(nb_points, Boot.naccu, times);

      for index=1:times
         if rem(index,10) == 0,  fprintf(' %d',index); end
         if rem(index,120) == 0, fprintf('\n'); end
         for allt=1:trials
            j=1;
            while j<=Boot.naccu
               t = ceil(rand([1 2])*trials); % random ints [1,g.timesout]
               if (allRn(t(1),index) == 1) & (allRn(t(2),index) == 1)
                  tmpsX(:,j) = squeeze(alltmpsX(t(1),index,:));
                  tmpsY(:,j) = squeeze(alltmpsY(t(2),index,:));
                  j = j+1;
               end
            end
            Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
         end;
         Boot.Coherboot = cohercomppost(Boot.Coherboot);  % CHECK IF NECSSARY FOR ALL BOOT TYPE
         Boot.fullcoherboot(:,:,index) = Boot.Coherboot.R;
         Boot.Coherboot = coherinit(nb_points, trials, Boot.naccu, Boot.Coherboot.type);
      end;
      Boot.Coherboot.R = Boot.fullcoherboot;
      Boot = rmfield(Boot, 'fullcoherboot');
   elseif strcmp(Boot.boottype, 'timestrials') % handle timestrials bootstrap
      fprintf('\nProcessing time and trial bootstrap (of %d):',trials);
      tmpsX = zeros(size(alltmpsX,3), Boot.naccu);
      tmpsY = zeros(size(alltmpsY,3), Boot.naccu );
      for allt=1:trials
         if rem(allt,10) == 0,  fprintf(' %d',allt); end
         if rem(allt,120) == 0, fprintf('\n'); end
         j=1;
         while j<=Boot.naccu
            t = ceil(rand([1 2])*trials); % random ints [1,g.timesout]
            goodbasewins = find((allRn(t(1),:) & allRn(t(2),:)) ==1);
            if Boot.baseboot % use baseline windows only
               goodbasewins = find(goodbasewins<=baselength);
            end
            ngdbasewins = length(goodbasewins);

            if ngdbasewins>1
               s = ceil(rand([1 2])*ngdbasewins); % random ints [1,g.timesout]
               s=goodbasewins(s);

               if all(allRn(t(1),s(1)) == 1) & all(allRn(t(2),s(2)) == 1)
                  tmpsX(:,j) = squeeze(alltmpsX(t(1),s(1),:));
                  tmpsY(:,j) = squeeze(alltmpsY(t(2),s(2),:));
                  j = j+1;
               end
            end
         end
         Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
      end
      Boot.Coherboot = cohercomppost(Boot.Coherboot);
   elseif strcmp(Boot.boottype, 'times') % boottype is 'times'
      Boot.Coherboot = cohercomppost(Boot.Coherboot);
   end;
end;

% test if precomputed
if ~isnan(Boot.alpha) & isnan(Boot.rboot) % if bootstrap analysis included . . .
   % 'boottype'='times' or 'timestrials', size(R)=nb_points*naccu
   % 'boottype'='trials',                 size(R)=nb_points*naccu*times
   Boot.Coherboot.R = abs (Boot.Coherboot.R);
   Boot.Coherboot.R = sort(Boot.Coherboot.R,2);

   % compute bootstrap significance level
   i = round(Boot.naccu*Boot.alpha);
   Boot.Rsignif = mean(Boot.Coherboot.R(:,Boot.naccu-i+1:Boot.naccu),2); % significance levels for Rraw
   Boot.Coherboot.R = squeeze(mean(Boot.Coherboot.R(:,Boot.naccu-i+1:Boot.naccu),2));
   if size(Boot.Coherboot.R, 2) == 1
	   Rbootout(:,2) = Boot.Coherboot.R;
   else
	   Rbootout(:,:,2) = Boot.Coherboot.R;
   end;
   % BEFORE
   %Rboot = [mean(Rboot(1:i,:)) ; mean(Rboot(g.naccu-i+1:g.naccu,:))];
elseif ~isnan(Boot.rboot)
	Boot.Coherboot.R = Boot.rboot;
	Boot.Rsignif     = Boot.rboot;
	Rbootout         = Boot.rboot;
else
	Boot.Coherboot.R = [];
	Boot.Rsignif     = [];
	Rbootout         = [];
end % NOTE: above, mean ?????

function w = hanning(n)
if ~rem(n,2)
   w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
   w = [w; w(end:-1:1)];
else
   w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
   w = [w; w(end-1:-1:1)];
end

author: PatternRecognition

repo: OpenBMI


text:
function J=JacobPolarCoordTurn2D(T,x,turnType,discPoint,tauTurn,tauLinAccel)
%%JACOBCOORDTURN2D Evaluate the Jacobian (matrix of partial derivatives
%              with respect to the state elements) of a discrete-time state
%              transition function for a two-dimensional coordinated turn
%              model where the velocity is specified in terms of a heading
%              and a speed. The turn rate can be specified in terms of a
%              turn rate in radians per second, or in terms of a
%              transversal acceleration. Additionally, a linear
%              acceleration can be given. The turn rate and linear
%              acceleration can optionally have time constants associated
%              with them, like in the Singer model, modelling a tendancy to
%              eventually want to return to non-accelerating, straight-line
%              motion. This Jacobian is associated with the continuous-time
%              model fTransPolarCoordTurn2D.
%
%INPUTS: T The time-duration of the propagation interval in seconds.
%        x The  target state for 2D motion where the velocity is given in
%          terms of heading and speed components. If there is no linear
%          acceleration (acceleration along the direction of motion),
%          then x can either be x=[x;y;h;v;omega], where h is the heading
%          in terms of radians counterclockwise from the x-axis, v is the
%          speed, and omega is the turn rate (the derivative of h with
%          respect to time) or  x=[x;y;h;v;at] where at is the
%          transversal acceleration, which is orthogonal to the velocity
%          and is defined such that positive values of at map to positive
%          values of omega. If there is a linear acceleration, then the
%          target state is either x=[x;y;h;v;omega;al] where omega is the
%          turn rate and al is the linear acceleration or the target
%          state is x=[x;y;h;v;at;al] if the turn is expressed in terms
%          of a transversal acceleration. The dimensionality of the state
%          is used to determine whether a linear acceleration component
%          is present. The linear acceleration component changes the
%          speed. That means that it is the derivative of the speed.
% turnType A string specifying whether the turn is given in terms of a
%          turn rate in radians per second or a transversal acceleration
%          in m/s^2. Possible values are
%          'TurnRate'   The turn is specified in terms of a turn rate
%                       (The default if this parameter is omitted).
%          'TransAccel' The turn is specified in terms of a transversal
%                       acceleration.
% discPoint This optional parameter specified what value of the turn rate
%           is used for the discretized state prediction. The thre
%           possible values were suggested in Li's paper cited
%           below. Possible values are
%           0 (The default if omitted) use omega=x(5), or the equivalent
%             value when specifying the turn rate using a transverse
%             acceleration, from the non-predicted target state for
%             building the state transition matrix. \omega_k
%           1 Use the average value of the predicted omega (or the average
%             value of the transverse acceleration) over the interval T
%             for building the state transition matrix. \bar{\omega}
%           2 Use the approximate average value of the predicted omega
%             over the interval T for building the state transition
%             matrix. \bar{\omega} This is the suggestion of using half
%             the prior and prediction, as was given in Li's paper. When
%             given a transverse acceleration instead of a turn rate, half
%             of the prior and predicted accelerations is used.
%           3 Use the forward-predicted omega (or transverse acceleration)
%             for building the state transition matrix. \omega_{k+1}
%   tauTurn The correlation time constant for the turn rate in seconds.
%           tau must be positive but does not have to be finite. If this
%           parameter is omitted, then tauTurn is set to infinity.
% tauLinAccel The correlation time constant for the linear acceleration (if
%            present) in seconds. This parameter is not needed if there is
%            no linear acceleration. If a linear acceleration is present
%            and this parameter is omitted, then tauLinAccel is set to
%            infinity.
%
%OUTPUTS: J The 5X5 or 6X6 Jacobian matrix of the function
%           fTransPolarCoordTurn2D, where the partial derivatives are
%           ordered [df/dx(1), df/dx(2),...,df/dx(xDim)]. That is, column i
%           consists of partial derivatives with respect to element i of
%           the x vector.
%
%More information on the discrete-time turning model is given in the
%comments to the function fTransPolarCoordTurn2D.
%
%April 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The time constant for the linear acceleration.
if(nargin<6||isempty(tauLinAccel))
    tauLinAccel=Inf;
end

%The time constant for the turn.
if(nargin<5||isempty(tauTurn))
   tauTurn=Inf;
end

if(nargin<4||isempty(discPoint))
    discPoint=0;
end

if(nargin<3||isempty(turnType))
    turnType='TurnRate';
end

xDim=size(x,1);

beta=exp(-T/tauTurn);
switch(discPoint)
    case 0%Use the prior value of at/omega for the linearization.
        turnVal=x(5);
    case 1
        %Use the average value of at.omega over the prediction
        %interval for the linearization.
        turnVal0=x(5);
        turnVal=(tauTurn-beta*tauTurn+T*turnVal0)/T;

        %Assume tau=Inf or T=0 was used, in which case the
        %asymptotic average value of at/omega should be used.
        turnVal(~isfinite(tauTurn)||~isfinite(turnVal))=turnVal0(~isfinite(tauTurn)||~isfinite(turnVal));
    case 2%Use the simple mean of at/omega for the linearization.
        turnVal=(beta+1)*x(5)/T;
    case 3%Use the predicted value of at/omega for the linearization.
        turnVal=beta*x(5);
    otherwise
        error('Invalid value entered for discPoint');
end

v=x(4);%The speed

switch(turnType)
    case 'TransAccel'%The turn is expressed in terms of a transverse
        %acceleration.
        omega=turnVal./v;
        omega(~isfinite(omega))=0;%Deal with zero velocity
    case 'TurnRate'%The turn is expressed in terms of a turn rate.
        omega=turnVal;
    otherwise
        error('Unknown turn type specified.');
end
%We now have the omega term for the turn.

theta=x(3);%The heading (counterclockwise from the x axis).

sinHead=sin(theta);
cosHead=cos(theta);

sinVal=sin(omega*T);
cosVal=cos(omega*T);

sinRat=sinVal./omega;
sinRat(~isfinite(sinRat))=T;%The limit as omega goes to zero.
cosRat=(1-cosVal)/omega;
cosRat(~isfinite(cosRat))=0;%The limit as omega goes to zero.

switch(turnType)
    case 'TransAccel'
        at=x(5);
        J=[1,0,-v*(sinHead*sinRat+cosHead*cosRat),2*(cosHead*sinRat-sinHead*cosRat)-T*(cosHead*cosVal-sinHead*sinVal),cosHead*(T*cosVal/omega-sinRat/omega)-sinHead*(T*sinRat-cosRat/omega);
           0,1,v*(cosHead*sinRat-sinHead*cosRat), 2*(sinHead*sinRat+cosHead*cosRat)-T*(sinHead*cosVal+cosHead*sinVal),sinHead*(T*cosVal/omega-sinRat/omega)+cosHead*(T*sinRat-cosRat/omega);
           0,0,1,                                 -T*at/v^2,                                                          T/v;
           0,0,0,                                 1,                                                                  0;
           0,0,0,                                 0,                                                                  beta];
    case 'TurnRate'
        J=[1,0,-v*(sinHead*sinRat+cosHead*cosRat),cosHead*sinRat-sinHead*cosRat,v*(cosHead*(T*cosVal/omega-sinRat/omega)-sinHead*(T*sinRat-cosRat/omega));
           0,1,v*(cosHead*sinRat-sinHead*cosRat), sinHead*sinRat+cosHead*cosRat,v*(sinHead*(T*cosVal/omega-sinRat/omega)+cosHead*(T*sinRat-cosRat/omega));
           0,0,1,                                 0,                            T;
           0,0,0,                                 1,                            0;
           0,0,0,                                 0,                            beta];
end
if xDim==6
    betaLinAccel=exp(-T/tauLinAccel);
    J(4,6)=T;
    J(6,6)=betaLinAccel;
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function poly2 = clipPolygon(polygon, w)
%CLIPPOLYGON Clip a polygon with a rectangular box
%
%   POLY2 = clipPolygon(POLY, BOX);
%   POLY is [Nx2] array of points
%   BOX has the form: [XMIN XMAX YMIN YMAX].
%   Returns the polygon created by the itnersection of the polygon POLY and
%   the bounding box BOX.
%
%   Note: Works only for convex polygons at the moment.
%
%   See also:
%   polygons2d, boxes2d, clipPolygonHP
%

% ---------
% author : David Legland
% created the 14/05/2005.
% Copyright 2010 INRA - Cepia Software Platform.

%   HISTORY
%   2007/09/14 fix doc

% check case of polygons stored in cell array
if iscell(polygon)
    poly2 = cell(1, length(polygon));
    for i=1:length(polygon)
        poly2{i} = clipPolygon(polygon{i}, w);
    end
    return;
end

% check case of empty polygon
N = size(polygon, 1);
if N==0
    poly2 = zeros(0, 2);
    return
end

% create edges array of polygon
edges = [polygon polygon([2:N 1], :)];

% clip edges
edges = clipEdge(edges, w);

% select non empty edges, and get their vertices
ind = sum(abs(edges), 2)>1e-14;
pts = unique([edges(ind, 1:2); edges(ind, 3:4)], 'rows');

% add vertices of window corner
corners = [w(1) w(3); w(1) w(4);w(2) w(3);w(2) w(4)];
ind = inpolygon(corners(:,1), corners(:,2), polygon(:,1), polygon(:,2));
pts = [pts; corners(ind, :)];

% polygon totally outside the window
if size(pts, 1)==0
    poly2 = pts;
    return;
end

% compute centroid of visible polygon
pc = centroid(pts);

% sort vertices around polygon
angle = edgeAngle([repmat(pc, [size(pts, 1) 1]) pts]);
[dummy, I] = sort(angle); %#ok<ASGLU>

% create resulting polygon
poly2 = pts(I, :);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function L = isinf(x)
%PREAL/ISINF    Overloaded ISINF function for class preal.

L=isinf(double(x));

author: Sable

repo: mcbench-benchmarks


text:
% Use is m.file to test the frequency calculater function for that you will
% need the Angle.mat file which contain a real measured signal with
% freuqncy changing between 50 Hz to 55% in three differnt ways:
% 1- A jump (step change) of the frequency
% 2- ramp change of the freuqnecy with 10 Hz/s
% 3- ramp change of the frequency with 20 Hz/s
% for using this file you will need the: "Frequency_Calculation_Function"
% Author: Aubai Al Khatib datum: 19.09.2013
clear
clc
if 1 == 0
    Sample = 10000;%Hz
    f_test = 50;%Hz
    t = 0:1/Sample:1;
    y = sin(2*pi*f_test*t);
    [f,f_time] = Frequency_Calculation_Function(y,t);
    figure(1);
    subplot(1,1,2);
    plot(f_time,f);ylim([49,56]);grid on;xlim([80,120]);xlabel('Time in S');ylabel('Frequency in Hz');
    subplot(1,2,2);
    plot(t,y);grid on;xlabel('Time in S');ylabel('Amplitude in V');
else
    load('Angle.mat')
    [f,f_time] = Frequency_Calculation_Function(U1,U1_time);
    figure(1);
    subplot(2,1,1);
    plot(f_time,f);ylim([49,56]);grid on;xlim([80,120]);xlabel('Time in S');ylabel('Frequency in Hz');
    subplot(2,1,2);
    plot(U1_time,U1);grid on;xlim([80,120]);xlabel('Time in S');ylabel('Amplitude in V');
    figure(2);
    subplot(2,1,1);
    plot(f_time,f);ylim([49,56]);grid on;xlim([104.5,108.5]);xlabel('Time in S');ylabel('Frequency in Hz');
    subplot(2,1,2);
    plot(U1_time,U1);grid on;xlim([104.5,108.5]);xlabel('Time in S');ylabel('Amplitude in V');
end


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the target point positions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function targets = update_Target_Point_Positions(dt,current_time,targets)

IDs = targets(:,1);                 % Stores Lag-Pt IDs in col vector
xPts= targets(:,2);                 % Original x-Values of x-Target Pts.
yPts= targets(:,3);                 % Original y-Values of y-Target Pts.
kStiffs = targets(:,4);             % Stores Target Stiffnesses

N_target = length(targets(:,1));    %Gives total number of target pts!

% Coefficients for Polynomial Phase-Interpolation
a = 2.739726027397260;  % y1(t) = at^2
b = 2.739726027397260;  % y3(t) = -b(t-1)^2+1
c = -2.029426686960933; % y2(t) = ct^3 + dt^2 + gt + h
d = 3.044140030441400;
g = -0.015220700152207;
h = 0.000253678335870;

% Period Info
tP1 = 0.025;                        % Contraction time
tP2 = 0.025;                        % Expansion time
period = tP1+tP2;                   % Period
t = rem(current_time,period);       % Current time in simulation ( 'modular arithmetic to get time in period')

% Read In Pts!
[xP1,yP1,xP2,yP2] = read_File_In('All_Positions.txt');

for i=1:N_target                    % Loops over all target points!

    if (t <= tP1)

			%PHASE 1 --> PHASE 2

			tprev = 0.0;
			t1 = 0.1*tP1;
			t2 = 0.9*tP1;
			if (t<t1) 							%For Polynomial Phase Interp.
				g1 = a*power((t/tP1),2);
            elseif ((t>=t1)&&(t<t2))
				g1 = c*power((t/tP1),3) + d*power((t/tP1),2) + g*(t/tP1) + h;
            elseif (t>=t2)
				g1 = -b*power(((t/tP1) - 1),2) + 1;
            end

			xPts(IDs(i)) = xP1(i) + g1*( xP2(i) - xP1(i) );
			yPts(IDs(i)) = yP1(i) + g1*( yP2(i) - yP1(i) );

    elseif ((t>tP1)&&(t<=(tP1+tP2)))

			%PHASE 2 --> PHASE 1

			tprev = tP1;
			t1 = 0.1*tP2 + tP1;
			t2 = 0.9*tP2 + tP1;
			if (t<t1) 							%//For Polynomial Phase Interp.
				g2 = a*power( ( (t-tprev)/tP2) ,2);
            elseif ((t>=t1)&&(t<t2))
				g2 = c*power( ( (t-tprev)/tP2) ,3) + d*power( ((t-tprev)/tP2) ,2) + g*( (t-tprev)/tP2) + h;
            elseif (t>=t2)
				g2 = -b*power( (( (t-tprev)/tP2) - 1) ,2) + 1;
            end

			xPts(IDs(i)) = xP2(i) + g2*( xP1(i) - xP2(i) );
			yPts(IDs(i)) = yP2(i) + g2*( yP1(i) - yP2(i) );

    end

    targets(IDs(i),2) = xPts(IDs(i)); % Store new xVals
    targets(IDs(i),3) = yPts(IDs(i)); % Store new yVals

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: reads in info from file that contains probability distributions
% (rows) for each game of Bingo w/ N players (columns)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x1,y1,x2,y2] = read_File_In(file_name)

filename = file_name;  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

mat_info = C{1};   %Stores all read in data

%Store all elements in matrix
mat = mat_info(1:end,1:end);

x1 =  mat(:,1); %store regular bingo expectation values
y1 =  mat(:,2); %store inner bingo expectation values
x2 =  mat(:,3); %store outer bingo expectation values
y2 =  mat(:,4); %store cover all bingo expectation values

author: nickabattista


text:
function A = area(grains,varargin)
% calculates the area of a list of grains
%
% Input
%  grains - @grain2d
%
% Output
%  A  - list of areas (in measurement units)
%

A = zeros(length(grains.poly),1);

poly = grains.poly;
V = grains.V;

for ig = 1:length(poly)
  A(ig) = polySgnArea(V(poly{ig},1),V(poly{ig},2));
end

author: mtex-toolbox

repo: mtex

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


text:
% DEMO4 - Playing with NetCDF-JAVA functions
echo('on')
% Starting DEMO4 ----------------------------------------------------------
%% ---- Example of data access using the NetCDF-Java API
ds = ncdataset('http://geoport.whoi.edu/thredds/dodsC/examples/OS_MBARI-M1_20091020_R_TS.nc');
nc = ds.netcdf;

%% ---- Query NetCDF using Java API
nc.getConventionUsed()
nc.getDetailInfo()
nc.hasUnlimitedDimension()
nc.getUseNaNs()

%% ---- Set properties using Java API
% The following line turns off the auto conversion of missing_values to NaN
% NOTE: This switches ALL netcdf files not just the one you are working on.
% nc.setUseNaNs(0)

%% ---- Dump NcML to file using Java API
out = java.io.BufferedWriter(java.io.FileWriter('trash_me.xml'));
nc.writeNcML(out, '')
out.close()
echo('off') % Ending DEMO4 ------------------------------------------------

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function v3d_write_volume(I,fname,scales)
% Function for writing V3D volume files version R6.1
%
% v3d_write_volume(volume, filename, voxelsize in mm)
%
% examples:
% I=uint16(rand(64,64,64)*65536);
%
% 1: v3d_write_volume(I);
% 2: v3d_write_volume(I,'random.v3d',[2 2 2];

% Filename
    if((exist('fname','var')==0)),
        [filename, pathname] = uiputfile('*.v3d', 'Write v3d-file');
        fname = [pathname filename];
    end

% Sizes
    sizes=size(I);
% Scales
    if(exist('scales','var')==0), scales=ones(1,3); end;
% Offset
    offset=204; % Header size R6.1
% Filesize
    fsize=numel(I)*2+offset;
% Format
    fform='3D-RA';
% Version
    vers='6.1';
% Voxelbits
    bits=16;
% par1;
    par1 = 1072693248; % The meaning of this number is unknown

disp(['filename : ' num2str(fname)]);
disp(['format : ' fform]);
disp(['version : ' vers]);
disp(['filesize : ' num2str(fsize)]);
fprintf('sizes : %i, %i, %i\n',sizes);
fprintf('scales : %2.6f, %2.6f, %2.6f\n',scales);
disp(['voxelbits : ' num2str(bits)]);
disp(['offset : ' num2str(offset)]);
disp(['par1 : ' num2str(par1)]);
fprintf('\n');

fout=fopen(fname,'wb');
fwrite(fout,[fform ' R' vers],'char');
fwrite(fout,zeros(1,30),'uint8'); %seek
fwrite(fout,sizes,'int');
fwrite(fout,scales,'double');

fwrite(fout,zeros(1,28),'uint8'); %seek
fwrite(fout,par1,'int');
fwrite(fout,zeros(1,28),'uint8'); %seek
fwrite(fout,par1,'int');
fwrite(fout,zeros(1,28),'uint8'); %seek
fwrite(fout,par1,'int');

fwrite(fout,bits,'int'); % Seems the new (R6.1) number of bits position

% Meaning of this part of the header is unknown
fwrite(fout,140,'int');
fwrite(fout,1,'int');
fwrite(fout,1075789855,'int');
fwrite(fout,1068449823,'int');
fwrite(fout,0,'int');
fwrite(fout,0,'uint8');
fwrite(fout,64,'uint8');
fwrite(fout,143,'uint8');
fwrite(fout,192,'uint8');
fwrite(fout,1,'int');

% Write uint16 volume
fwrite(fout,uint16(I),'uint16');

fclose('all');

author: tomdoel

repo: pulmonarytoolkit


text:
function SO3F = times(SO3F1,SO3F2)
% overloads |SO3F1 .* SO3F2|
%
% Syntax
%   sF = SO3F1 .* SO3F2
%   sF = a .* SO3F2
%   sF = SO3F1 .* a
%
% Input
%  SO3F1, SO3F2 - @SO3FunCBF
%  a - double
%
% Output
%  SO3F - @SO3Fun
%

if isnumeric(SO3F1)
  SO3F = SO3F2;
  SO3F.weights = SO3F.weights .* SO3F1;
  return
end

SO3F = times@SO3Fun(SO3F1,SO3F2);

end


text:
function HierarchicalClusteringBiox2(FiberGroupFile, numNodesToResample, NumFibersPerPartition, maxDistanceOfInterest, NClusters)

%Make all parameters numerical;
if ischar(numNodesToResample)
    numNodesToResample=str2num(numNodesToResample);
end

%Make all parameters numerical;
if ischar(NumFibersPerPartition)
    NumFibersPerPartition=str2num(NumFibersPerPartition);
end

%Make all parameters numerical;
if ischar(maxDistanceOfInterest)
    maxDistanceOfInterest=str2num(maxDistanceOfInterest);
end

%Make all parameters numerical;
if ischar(NClusters)
    NClusters=str2num(NClusters);
end

warning off;

[pathstr, FiberGroupFileName, ext, versn] = fileparts(FiberGroupFile);
if isempty(pathstr)
    pathstr=['.' filesep ];
end

%0 Reorder fibers by distance among seeds
load(FiberGroupFile);
%fibersReorderedByDistance=dtiReorderFibersByDistance(fg)
fg=dtiReorderFibersByDistance(fg);
ReorderedFiberGroupFileName=[FiberGroupFileName 'ReorderedByDistance'];
save(ReorderedFiberGroupFileName, 'fg', 'versionNum', 'coordinateSpace');
 display('Reordered by Distance');

%1. Resample originalFibersFile ->ResampledFibers
%Make all parameters numerical;
if ischar(numNodesToResample)
    numNodesToResample=str2num(numNodesToResample);
end

fg=resampleFiberGroup(fg, numNodesToResample);

ResampledFolder=['FibersResampledTo' num2str(numNodesToResample)];
mkdir(pathstr, ResampledFolder);
ResampledFiberGroupFile=fullfile(pathstr, ResampledFolder, [ReorderedFiberGroupFileName 'Resampled']);
save(ResampledFiberGroupFile, 'fg', 'coordinateSpace', 'versionNum');
clear fg coordinateSpace versionNum; %Saving memoryl
display('Resampled');

%2. Partition ResampledFiberGroupFile by chunks of
%size=NumFibersPerPartition
cd(ResampledFolder)
dtiFractionizeFiberDataSet([ReorderedFiberGroupFileName 'Resampled'], NumFibersPerPartition);

cd('fgFile_Parts')

%3. Go through FIBER partitions and cluster them, creating master clusterLabels
clusterlabelsAll=[];
allPartitions=dir([ReorderedFiberGroupFileName 'Resampled' '*.mat']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%THIS PART IS SPECIFIC TO BIOX2CLUSTER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for PartitionId=1:size(allPartitions, 1)

%[clusterlabels]=ClusterFibers(allPartitions(PartitionId).name, maxDistanceOfInterest, NClusters);
%the two lines below perform clustering by submitting a job ^^
unix(['/home/elenary/project_scripts/full_brain_clustering/cluster_fibers.batch.job ' allPartitions(PartitionId).name ' ' num2str(maxDistanceOfInterest) ' ' num2str(NClusters)]);
display(['Fibers clustered for partition ' allPartitions(PartitionId).name]);
end

%Check if the jobs are completed, then load computed cluster labels
a=[];
while size(a, 1)<size(allPartitions, 1)
    a=dir('*clusterlabels.mat');
    pause(2);
end
        unix(['ls *clusterlabels.mat > ' ReorderedFiberGroupFileName 'Resampled.txt']);

for PartitionId=1:size(allPartitions, 1)

[pathstrCL, NameCL, extCL, versnCL] = fileparts(allPartitions(PartitionId).name);
load([NameCL 'clusterlabels']);
clusterlabelsAll=[clusterlabelsAll; clusterlabels+NClusters*(PartitionId-1)];
end

%cleanup temp files
unix(['rm *clusterlabels*']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%END OF CLUSTER SPECIFIC PART%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%RECURSIVE PART START
%3. Create COMPUTE SUPERFIBERREPRESENTATION for the full original FiberSet
%DECIDE ON: levels; NumSuperFibersPerPartition, NClusters on every level
level=0;
    NumSuperFibersPerPartition=NumFibersPerPartition;  NSFClusters=NClusters; %maxDistanceOfInterest=50; %5cm
display('Start recursion');

while size(allPartitions, 1)>1
    level=level+1;

    cd('../../'); %To go back up from fg_FileParts -up, ResampledFolder -->up
    load(ResampledFiberGroupFile);
    display(['Computing superfiber representations at level ' num2str(level)]);
    tic; SuperFibersGroup = dtiComputeSuperFiberRepresentation(fg, clusterlabelsAll); toc;
    display('Superfibers Created');

    fg=SuperFibersGroup;
    ParentFiberGroupFile=[ResampledFiberGroupFile '.mat'];
    SuperFiberGroupFileName=['SuperFibersOf_' ReorderedFiberGroupFileName 'Resampled_Level' num2str(level)];
    SuperFiberGroup=fullfile(pathstr, ResampledFolder, SuperFiberGroupFileName);
    save(SuperFiberGroup, 'ParentFiberGroupFile', 'clusterlabelsAll', 'maxDistanceOfInterest', 'fg', 'coordinateSpace', 'versionNum');

    %Cluster the SuperFibers and update master clusterlabelsAll
    cd(ResampledFolder)
    FractionizeFiberDataSet(SuperFiberGroupFileName, NumSuperFibersPerPartition); %Partitions will be saved un fgFile_Parts dir
    cd('fgFile_Parts')
    %Go through FIBER partitions and cluster them, creating master clusterLabels
    %We have here clusterlabelsALL
    allPartitions=dir([SuperFiberGroupFileName '*mat']);
    clusterlabelsCurrAll=clusterlabelsAll; %Labels from the current level

    %%%%%%%%%%%THIS PART IS SPECIFIC TO BIOX2 CLUSTER
    %%%HERE WE ARE SUBMITTING BJOBS

    for SFPartitionId=1:size(allPartitions, 1)

       % [clusterlabels]=ClusterFibers(allPartitions(SFPartitionId).name, maxDistanceOfInterest, NSFClusters);
       unix(['/home/elenary/project_scripts/full_brain_clustering/cluster_fibers.batch.job ' allPartitions(SFPartitionId).name ' ' num2str(maxDistanceOfInterest) ' ' num2str(NSFClusters)]);
       display(['SuperFibers clustered for partition ' allPartitions(SFPartitionId).name]);
        %Udate the portion of clusterlabelsAll pertaining to current partition
        %of superfiberset only

    end

    %%%%%%%%%%%THIS PART IS SPECIFIC TO BIOX2 CLUSTER
    %%%HERE WE ARE READING RESULS OF SUBMITTED JOBS

        %Check if the jobs are completed, then load computed cluster labels
        a=[];
        while size(a, 1)<size(allPartitions, 1)
            a=dir('*clusterlabels.mat');
            pause(2);
        end
        unix(['ls *clusterlabels.mat > SuperFibersOf_' ReorderedFiberGroupFileName 'Resampled_Level' num2str(level) '.txt']);

    for SFPartitionId=1:size(allPartitions, 1)

        %Udate the portion of clusterlabelsAll pertaining to current partition
        %of superfiberset only
        [pathstrCLSF, NameCLSF, extCLSF, versnCLSF] = fileparts(allPartitions(SFPartitionId).name);
        load([NameCLSF 'clusterlabels']);

        for SuperFiberId=1:min(NumSuperFibersPerPartition, size(clusterlabels, 1))
            clusterlabelsAll(find(clusterlabelsCurrAll==SuperFiberId+(SFPartitionId-1)*NumSuperFibersPerPartition), 1)=clusterlabels(SuperFiberId)+NSFClusters*(SFPartitionId-1);
        end

    end

        %cleanup temp files
        unix(['rm *clusterlabels*']);

end %while size(allPartitions, 1)>1

%So we have a single partition of SuperFibers for which clustering has just
%been performed. Save the results.

    level=level+1;

    cd('../../'); %To go back up from fg_FileParts -up, ResampledFolder -->up
    load(ResampledFiberGroupFile);
    SuperFibersGroup = dtiComputeSuperFiberRepresentation(fg, clusterlabelsAll);
    display('Superfibers Created');

    fg=SuperFibersGroup;
    ParentFiberGroupFile=[ResampledFiberGroupFile '.mat'];
    SuperFiberGroupFileName=['SuperFibersOf_' ReorderedFiberGroupFileName 'Resampled_Level' num2str(level)];
    SuperFiberGroup=fullfile(pathstr, ResampledFolder, SuperFiberGroupFileName);
    save(SuperFiberGroup, 'ParentFiberGroupFile', 'clusterlabelsAll', 'maxDistanceOfInterest', 'fg', 'coordinateSpace', 'versionNum');
%Note: we are creating superfiber representation for
%ResampledFiberGroupFileName, not ReorderedFiberGroupFileName, because for
%computing SF representation the fibers need to be equisampled.

%However, to create a good visualization of clustering results, we are
%encouraged to use
%fg from ReorderedFiberGroupFileName
%and clusterlabelsAll from     SuperFiberGroupFileName

%fibergroupsvector=1:50;
%load(ReorderedFiberGroupFileName);
%dtiSaveFiberClusters(clusterlabelsAll, fibergroupsvector, fg, versionNum, coordinateSpace)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function w_lcmv = sphLCMV(sphCOV, constraint_dirs, constraints)
%SPHLCMV Get beamforming weights for a LCMV beamformer in the SHD
%
%   Computes the beamweights for a linearly constrained minimum variance
%   (LCMV) beamformer using SH signals.
%
%   Inputs:
%       sphCOV:     (order+1)^2x(order+1)^2 covariance/correlation matrix
%           of SH signals
%       constraint_dirs:    Kx2 [azi elev] directions that the beamformer
%           should achieve the user-specified constraints
%       constraints:        Kx1 constraint values
%
%   Outputs:
%       w_lcmv: (order+1)^2x1 vector of beamforming weights
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SPHLCMV.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nSH = size(sphCOV,1);
order = sqrt(nSH)-1;

% steering vectors for distortionless constraints
constraint_dirs2 = [constraint_dirs(:,1) pi/2-constraint_dirs(:,2)]; % from azi-elevation to azi-inclination
Y_lcmv = getSH(order, constraint_dirs2, 'real');

% LCMV weights vector
stVecMtx = Y_lcmv.';
invA_B = sphCOV\stVecMtx;
B_invA_B = stVecMtx' * invA_B;
w_lcmv = (invA_B / B_invA_B) * constraints;

end

author: polarch

repo: Spherical-Array-Processing


text:
function [ a, b ] = p03_ab ( m )

%*****************************************************************************80
%
%% P03_AB evaluates the limits of the optimization region for problem 03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Marcin Molga, Czeslaw Smutnicki,
%    Test functions for optimization needs.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real A(M), B(M), the lower and upper bounds.
%
  a(1:m) = - 65.536;
  b(1:m) = + 65.536;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_suite = test_daubcqf
initTestSuite;

function test_daubcqf_min
  [a, b] = daubcqf(4);
  ax = [0.482962913144534   0.836516303737808   0.224143868042013  -0.129409522551260];
  bx = [0.129409522551260   0.224143868042013  -0.836516303737808   0.482962913144534];
assertVectorsAlmostEqual(a, ax, 'relative', 0.001);
assertVectorsAlmostEqual(b, bx, 'relative', 0.001);

function test_daubcqf_max
  [a, b] = daubcqf(4, 'max');
  ax = [-0.129409522551260   0.224143868042013   0.836516303737808   0.482962913144534];
  bx = [-0.482962913144534   0.836516303737808  -0.224143868042013  -0.129409522551260];
assertVectorsAlmostEqual(a, ax, 'relative', 0.001);
assertVectorsAlmostEqual(b, bx, 'relative', 0.001);

function test_daubcqf_mid_even_k
  [a, b] = daubcqf(4, 'mid');
  ax = [0.482962913144534   0.836516303737808   0.224143868042013  -0.129409522551260];
  bx = [0.129409522551260   0.224143868042013  -0.836516303737808   0.482962913144534];
assertVectorsAlmostEqual(a, ax, 'relative', 0.001);
assertVectorsAlmostEqual(b, bx, 'relative', 0.001);

function test_daubcqf_mid_odd_k
  [a, b] = daubcqf(6, 'mid');
  ax = [0.332670552950083   0.806891509311093   0.459877502118491  -0.135011020010255  -0.085441273882027   0.035226291885710];
  bx = [-0.035226291885710  -0.085441273882027   0.135011020010255   0.459877502118491 -0.806891509311093   0.332670552950083];
assertVectorsAlmostEqual(a, ax, 'relative', 0.001);
assertVectorsAlmostEqual(b, bx, 'relative', 0.001);

function test_daubcqf_odd
  handle = @() daubcqf(9);
assertExceptionThrown(handle, '');

author: ricedsp


text:
function [omckk,Tckk,Rckk,H,x,ex,JJ] = compute_extrinsic(x_kk,X_kk,fc,cc,kc,alpha_c,MaxIter,thresh_cond),

%compute_extrinsic
%
%[omckk,Tckk,Rckk,H,x,ex] = compute_extrinsic(x_kk,X_kk,fc,cc,kc,alpha_c)
%
%Computes the extrinsic parameters attached to a 3D structure X_kk given its projection
%on the image plane x_kk and the intrinsic camera parameters fc, cc and kc.
%Works with planar and non-planar structures.
%
%INPUT: x_kk: Feature locations on the images
%       X_kk: Corresponding grid coordinates
%       fc: Camera focal length
%       cc: Principal point coordinates
%       kc: Distortion coefficients
%       alpha_c: Skew coefficient
%
%OUTPUT: omckk: 3D rotation vector attached to the grid positions in space
%        Tckk: 3D translation vector attached to the grid positions in space
%        Rckk: 3D rotation matrices corresponding to the omc vectors
%        H: Homography between points on the grid and points on the image plane (in pixel)
%           This makes sense only if the planar that is used in planar.
%        x: Reprojections of the points on the image plane
%        ex: Reprojection error: ex = x_kk - x;
%
%Method: Computes the normalized point coordinates, then computes the 3D pose
%
%Important functions called within that program:
%
%normalize_pixel: Computes the normalize image point coordinates.
%
%pose3D: Computes the 3D pose of the structure given the normalized image projection.
%
%project_points.m: Computes the 2D image projections of a set of 3D points

if nargin < 8,
   thresh_cond = inf;
end;

if nargin < 7,
   MaxIter = 20;
end;

if nargin < 6,
   alpha_c = 0;
	if nargin < 5,
   	kc = zeros(5,1);
   	if nargin < 4,
      	cc = zeros(2,1);
      	if nargin < 3,
         	fc = ones(2,1);
         	if nargin < 2,
            	error('Need 2D projections and 3D points (in compute_extrinsic.m)');
            	return;
         	end;
      	end;
   	end;
	end;
end;

% Initialization:

[omckk,Tckk,Rckk] = compute_extrinsic_init(x_kk,X_kk,fc,cc,kc,alpha_c);

% Refinement:
[omckk,Tckk,Rckk,JJ] = compute_extrinsic_refine(omckk,Tckk,x_kk,X_kk,fc,cc,kc,alpha_c,MaxIter,thresh_cond);

% computation of the homography (not useful in the end)

H = [Rckk(:,1:2) Tckk];

% Computes the reprojection error in pixels:

x = project_points2(X_kk,omckk,Tckk,fc,cc,kc,alpha_c);

ex = x_kk - x;

% Converts the homography in pixel units:

KK = [fc(1) alpha_c*fc(1) cc(1);0 fc(2) cc(2); 0 0 1];

H = KK*H;

return;

% Test of compte extrinsic:

Np = 4;
sx = 10;
sy = 10;
sz = 5;

om = randn(3,1);
T = [0;0;100];

noise = 2/1000;

XX = [sx*randn(1,Np);sy*randn(1,Np);sz*randn(1,Np)];
xx = project_points(XX,om,T);

xxn = xx + noise * randn(2,Np);

[omckk,Tckk] = compute_extrinsic(xxn,XX);

[om omckk om-omckk]
[T Tckk T-Tckk]

figure(3);
plot(xx(1,:),xx(2,:),'r+');
hold on;
plot(xxn(1,:),xxn(2,:),'g+');
hold off;

repo: ekfmonoslam


text:
function performBatchRFvasari_LGG(pathRF,nTrees,seed,matlabPATH)
% -------------------------------------------------------------------------
% function performBatchRFvasari_LGG(pathRF,nTrees,seed,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function performs prediction estimation using out-of-bag estimates
% for random forests constructed using (VASARI) features.
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathRF: Full path to the RANDOM_FORESTS results folder.
%            --> Ex: '/myProject/WORKSPACE/RANDOM_FORESTS'
% 2. nTrees: Number of decision-trees in the random forests.
%            --> Ex: 500
% 3. seed: Numerical number to use as seed for random-forest construction.
%          --> Ex: 54288
% 6. matlabPATH: Full path to the MATLAB executable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% -------------------------------------------------------------------------
% OUTPUTS: Final out-of-bag prediction estimation results saved in a table
% in 'pathRF'.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

cd(pathRF), load('dataStruct')
nameOutcomes = fieldnames(dataStruct); nOutcomes = numel(nameOutcomes);

time = 30; testCost = 0.25:0.25:5;
mkdir('batchLog_RFvasari'), cd('batchLog_RFvasari'), pathBatch = pwd;

% PRODUCE BATCH COMPUTATIONS
save('workspace','nTrees','dataStruct','nameOutcomes','seed','testCost'), batch = 0;
for o = 1:nOutcomes
    batch = batch + 1;
    nameScript = ['batch',num2str(batch),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    fprintf(fid,['nameOutcome = nameOutcomes{',num2str(o),'};\n']);
    fprintf(fid,['nVas = dataStruct.(nameOutcome).nVas;\n']);
    fprintf(fid,['nText = dataStruct.(nameOutcome).nText;\n']);
    fprintf(fid,['X = dataStruct.(nameOutcome).data(:,1:nVas);\n']);
    fprintf(fid,['Y = dataStruct.(nameOutcome).outcome;\n']);
    fprintf(fid,['cat = dataStruct.(nameOutcome).categories(1:nVas);\n']);
    fprintf(fid,['[auc,sensitivity,specificity,accuracy] = predictionEstimationRF_LGG(X,Y,logical(cat),nTrees,seed,testCost);\n']);
    fprintf(fid,['save(''result_batch',num2str(batch),''',''auc'',''sensitivity'',''specificity'',''accuracy'')\n']);
    fprintf(fid,['system(''touch batch',num2str(batch),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nOutcomes)
delete('workspace.mat')

% GROUP RESULTS
cd(pathBatch)
batch = 0; AUC = zeros(nOutcomes,1); Sensitivity = zeros(nOutcomes,1); Specificity = zeros(nOutcomes,1); Accuracy = zeros(nOutcomes,1);
for o = 1:nOutcomes
    batch = batch + 1;
    load(['result_batch',num2str(batch)]) % Variable 'auc', 'sensitivity', 'specificity', 'accuracy' gets out of here.
    AUC(o) = auc; Sensitivity(o) = sensitivity; Specificity(o) = specificity; Accuracy(o) = accuracy;
    delete(['result_batch',num2str(batch),'.mat'])
    clear auc sensitivity specificity accuracy
end
tab = table(AUC,Sensitivity,Specificity,Accuracy,'RowNames',nameOutcomes);
cd(pathRF), save('tableSummary_Vasari','tab')

cd(startpath)
end

author: mvallieres

repo: radiomics