text:
function [ anim info ] = bundleAdjustment( anim, varargin )
% Bundle Adjustment (using the excellent SBA)
%
% Can perform BA on full camera matrices (3 x 4) or K, R, t
%
% This can work for homographies too in case S0 is [ 2 x nPoint ]
%
% The calibration parameters are given in 3 or 5 by nFrame matrices.
% Each line corresponds to an element:
% Affine
%   K1  K2  0
%   0   K3  0
%   0   0   1
% Projective
%   K1  K2  K4
%   0   K3  K5
%   0   0   1
%
% If K is common to all cameras, it will not be optimized upon :( due to a
% limitation in SBA
%
% USAGE
%  [ anim ] = bundleAdjustment( anim, sfmCase, varargin )
%
% INPUTS
%  anim      - Animation object
%  varargin   - list of paramaters in quotes alternating with their values
%       - sfmCase   - 'motstr' for motion+struct, 'mot' motion only, 'str'
%              for structure only (this option is only for the rigid case
%       - 'KMask', [ 3 x 1 ] or [ 5 x 1 ] contains 1 when the calibration
%               parameter is fixed. By default, all of K is optimized upon
%       - 'nItr' number of BA iterations
%       - 'nFrameFixed' [1], number of frames (starting from the first)
%                       whose camera parameters are fixed
%
% OUTPUTS
%  anim      - optimized animation with optimized S, P, K, R, t
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

[ KMask nItr nFrameFixed sfmCase ] = ...
  getPrmDflt( varargin,{ 'KMask', [], 'nItr', 500, ...
  'nFrameFixed', 1, 'sfmCase', 'motstr'}, 1 );
KMaskOri=KMask;

% Figure out where the files are located
switch computer
  case {'PCWIN'},
    projPath=[ '@' fileparts(mfilename('fullpath')) ...
      '\private\sba\sbaProjection32.dll' ];
  case {'PCWIN64'},
    projPath=[ '@' fileparts(mfilename('fullpath')) ...
      '\private\sba\sbaProjection64.dll' ];
  otherwise,
    % Linux/Mac Matlab/Octave
    projPath=[ '@' fileparts(mfilename('fullpath')) ...
      '/private/sba/sbaProjection.so' ];
end

% set some parameters from the anim obejct
if isempty(anim.mask); WMask=ones(anim.nPoint, anim.nFrame);
else WMask=anim.mask; end

KAll=anim.K;
if size(KAll,2)==1 && ~isempty(KAll)
  KAll=repmat(anim.K,1,anim.nFrame);
end

% Initialize the camera parameters with what already exists
% if there is nothing, use the full projection matrices
if ~isempty(anim.R) && ~isempty(anim.t)
  % create P0 for rotation and stuff
  if isempty(KAll)
    K0=zeros(0,anim.nFrame); nK = 0; doK=0;
    if anim.isProj; KMask = zeros(5,1); else KMask = zeros(3,1); end
  else
    doK=1;
    if isempty(KMaskOri)
      if anim.isProj; KMask = ones(5,1); else KMask = ones(3,1); end
    end
    if size(anim.K,2)==1
      % do not optimize over K if it is a common K (we will optimize later)
      K0=zeros(0,anim.nFrame); nK = 0;
    else
      K0=KAll(~KMask,:); nK = length(find(~KMask));
    end
  end

  % Deal with NRSFM
  switch size(anim.l,1)
    case anim.nBasis, %Xiao
      isFirstCoeff1=false;
    case anim.nBasis-1, %Torresani
      isFirstCoeff1=true;
    otherwise
      if ~isempty(anim.l); error([ 'Problem with the dimension of l ' ...
          'and SBasis' ]);end
  end

  % Get the quaternion from the rotation matrices
  Q = quaternion(anim.R);

  %  sba(n, m, mcon, vmask, p0, cnp, pnp, x, covx, mnp, proj, projac, ...
  %itmax, verbose, opts, reftype, varargin)
  pnp=3;
  if anim.isProj % Projective camera
    P0 = [ Q; anim.t; K0 ]; cnp=7+nK;
  else % Affine camera
    P0 = [ Q; anim.t(1:2,:); K0 ]; cnp=6+nK;
  end

  if ~isempty(anim.l) % NRSFM
    P0 = [ P0; anim.l ];
    P0 = [ P0(:)' reshape( permute( anim.SBasis, [ 1 3 2 ] ), [], 1 )' ];
    if anim.isProj
      error('has to be an affine camera for NRSFM');
    end
    cnp=6+size(anim.l,1);
    pnp=3*anim.nBasis;
    proj='affineNRSFM';
    sfmCase='motstr';
  else
    P0 = [ P0(:)' anim.S(:)' ];
    if anim.isProj % Projective camera
      if doK
        proj='projectivekap1kap2pp1pp2Ignored';
      else
        proj='projectivek1k2k3k4k5kap1kap2pp1pp2Ignored';
      end
    else % Affine camera
      if doK
        proj='affinetr3k4k5kap1kap2pp1pp2Ignored';
      else
        proj='affinetr3k1k2k3k4k5kap1kap2pp1pp2Ignored';
      end
    end
  end
  fullP = false;
else
  fullP = true;
  if size(anim.P,2)==3
    %homography case
    proj = 'homography';
    P0 = reshape( anim.P, [ 9 anim.nFrame ] );
    P0 = bsxfun(@rdivide, P0(1:8,:), P0(9,:)); cnp=8;
    P0 = [ P0(:)' anim.S(:)' ];
    pnp = 2;
    nFrameFixed = 0;
  else
    if anim.isProj
      P0 = reshape( permute( anim.P, [ 2 1 3 ] ), [ 12 anim.nFrame ] ); cnp=12;
      proj = 'projectiveFull';
    else
      P0 = reshape( permute( anim.P(1:2,:,:), [ 2 1 3 ] ), [ 8 anim.nFrame ] ); cnp=8;
      proj = 'affineFull';
    end
    P0 = [ P0(:)' anim.S(:)' ];
    pnp = 3;
  end
end

% Initialize the point variables
WPermute=permute(anim.W, [1,3,2]);
x=WPermute(repmat(reshape(permute(WMask,[2,1]),1,...
  anim.nFrame,anim.nPoint), [2,1,1])>0);

% launch sba
projac=[ proj 'Jac' projPath ];
proj=[ proj projPath ];

if ~isempty(anim.l)
  % NRSfM
  isFirstCoeff1=double(isFirstCoeff1);
  nFrameFixed = 0;
  [ ret P info ] = sba( anim.nPoint, 0, anim.nFrame, nFrameFixed, WMask+0.0, ...
    P0, cnp, pnp, x, 2, proj, projac, nItr, 0, [], sfmCase, isFirstCoeff1,...
    anim.nBasis);
else
  % rigid SfM
  [ ret P info ] = sba( anim.nPoint, 0, anim.nFrame, nFrameFixed, WMask+0.0, ...
    P0, cnp, pnp, x, 2, proj, projac, nItr, 0, [], sfmCase,...
    KAll, KMask );
end

%  info(1:2)
%info(1:2)/anim.nFrame

% create the output information
if fullP
  if size(anim.P,2)==3
    % homography
    anim.S = reshape( P(8*anim.nFrame + 1 : end ), 2, anim.nPoint );
    P = reshape( P(1:8*anim.nFrame ), [ 8 anim.nFrame ] );
    P(end+1,:) = 1;
    anim.P = reshape( P, [ 3 3 anim.nFrame ] );
  else
    if anim.isProj
      P1 = permute( reshape( P(1:12*anim.nFrame ), ...
        [ 4 3 anim.nFrame ] ), [ 2 1 3 ] );
      anim.S = reshape( P(12*anim.nFrame + 1 : end ), 3, anim.nPoint );
    else
      P1 = permute( reshape( P(1:8*anim.nFrame ), ...
        [ 4 2 anim.nFrame ] ), [ 2 1 3 ] );
      P1(3,4,:) = 1;
      anim.S = reshape( P(8*anim.nFrame + 1 : end ), 3, anim.nPoint );
    end
    anim.P = P1;
  end
else
  if anim.isProj % Projective camera
    P1 = reshape( P(1:(7+nK)*anim.nFrame ), [], anim.nFrame );
    % the following line is just to prevent automatic update until
    % anim.t, anim.K and anim.R are all full
    anim.t = [];
    anim.R = quaternion( P1(1:4,:) );
    anim.t = P1(5:7,:);
    if size(P1,1)>=8; KAll( ~KMask, : ) = P1( 8:end, : ); end
    anim.S = reshape( P(( 7 + nK)*anim.nFrame + 1 : end ), 3, ...
      anim.nPoint );
  else % Affine camera
    if isempty(anim.l)
      P1 = reshape( P(1:(6+nK)*anim.nFrame ), [], anim.nFrame );
      % the following line is just to prevent automatic update until
      % both anim.t and anim.R are full
      anim.R = [];
      anim.t = P1(5:6,:); anim.t(3,:) = 0;
      anim.R = quaternion( P1(1:4,:) );
      if size(P1,1)>=7; KAll( ~KMask, : ) = P1( 7:end, : ); end
      anim.S = reshape( P(( 6 + nK)*anim.nFrame + 1 : end ), 3, ...
        anim.nPoint );
    else
      P1 = reshape( P(1:cnp*anim.nFrame ), cnp, anim.nFrame );
      % the following line is just to prevent automatic update until
      % both anim.t and anim.R (or anim.l and anim.SBasis) are full
      anim.t = []; anim.SBasis = [];
      anim.R = quaternion( P1(1:4,:) );
      anim.t = P1(5:6,:); anim.t(3,:) = 0;
      anim.l = P1( 7 : end, : );
      anim.SBasis = permute( reshape( P(cnp*anim.nFrame + 1 : end ), 3, ...
        anim.nBasis, anim.nPoint ), [ 1 3 2 ] );
    end
  end
end

if isempty(KMaskOri)
  if anim.isProj; KMaskOri = zeros(5,1); else KMaskOri = zeros(3,1); end
else
  KMaskOri = KMaskOri>0;
end

if ~isempty(anim.K) && any(~KMaskOri)
  if size(anim.K,2)==1
    % average all the K found as a first estimate
    anim.K=mean(KAll,2);
    err=anim.computeError(); errPrev=err(1);

    % optimize K and the rest through alternate gradient descent
    while 1
      S=reshape(anim.generateSAbsolute(),3,[]);
      % we want to solve K*S=W or kron(S',eye(2)) vec(K)=vec(W)
      % K is now only its first two rows
      if anim.isProj
        W=bsxfun(@times,reshape(anim.W,2,[]),S(3,:));
        A=kron(S',eye(2));
      else
        W=anim.W;
        A=kron(S(1:2,:)',eye(2));
      end

      % remove the column where K has a value of 0
      W=reshape(W,[],1);
      A(:,2)=[];

      % remove the columns for which we know the values
      if any(KMaskOri)
        W=W-A(:,KMaskOri)*anim.K(KMaskOri);
        A(:,KMaskOri)=[];
      end

      % solve for the best K
      anim.K(~KMaskOri)=A\W;

      % optimize the other parameters
      anim=bundleAdjustment(anim,'KMask',ones(size(anim.K,1),1),...
        'nItr',nItr,'nFrameFixed',nFrameFixed,'sfmCase',sfmCase);
      err=anim.computeError(); err=err(1);

      % exit when the error changes by less than 1%
      if (errPrev-err<0.01*errPrev); break; end
      errPrev=err;
    end
  else
    anim.K=KAll;
  end
end

author: vrabaud


text:
function I2= iir(filename,f,varargin)
%IIR Increases the resolution of an image by interpolation
% B= IIR(inputfile,f) returns the image stored in file 'inputfile' with
% resolution increased by factor f in both dimensions. 'filename' must be a
% valid graphic file (jpg, gif, tiff, etc.). It can be grayscale or color.
% Parameter 'f' is the size increase ratio, so to increase by 50%
% use f= 1.5, to double size (in each dimension) use f= 2.
%
% Additional parameters:
% B= IIR(A,f,'Display','off') eliminates display of both images, the original
% and the modified. Deafult 'on'
%
% B= IIR(A,f,'Method',method) Allows to choose between five methods of
% interpolation: linear, spline, pchip, cubic or v5cubic. 'method' must
% be a string character. Default 'linear'
%
% Example:
% B= iir('myimage.jpg',2,'Method','cubic');
%
% Last modified: Sep. 2010

% Defaults
method= 'linear';
displ= 'on';
f= max(f,1);
npass= 1;

% Extract optional arguments
for j= 1:2:length(varargin)
	switch varargin{j}
		case 'method'
			method= varargin{j+1};
		case 'display'
			displ= varargin{j+1};
		otherwise
			error('Unknown parameter name');
	end
end

% Read image file
I= imread(filename);

% -------------------------------------------
% Do the math
% Sizes and new image array
nrow= ceil(size(I,1)+ size(I,1)*(f-1)/npass);
ncol= ceil(size(I,2)+ size(I,2)*(f-1)/npass);
I2= uint8(zeros(nrow,ncol,size(I,3)));

% Loop through rows & cols
for j= 1:size(I,1)
	for c= 1:size(I,3)
		I2(j,:,c)= expand(double(I(j,:,c)),ncol,method);
	end
end
for j= 1:size(I2,2)
	for c= 1:size(I,3)
		I2(:,j,c)= expand(double(I2(1:size(I,1),j,c)),nrow,method);
	end
end

% Plot final images
if strcmp(displ,'on')
	if size(I,3) == 1
		figure(1), imagesc(I); colormap(gray); axis image, axis off
		figure(2), imagesc(I2);colormap(gray); axis image, axis off
	else
		figure(1), image(I); axis image, axis off
		figure(2), image(I2);axis image, axis off
	end
end

% ##################################################
% Support function
% ##################################################

function yy= expand(y,ndot,method)

x = 1:length(y);
xx = linspace(1,length(y),ndot);
switch method
	case 'linear'
		yy = uint8(interp1(y,xx,'linear'));
	case 'spline'
		yy = uint8(interp1(y,xx,'spline'));
	case 'pchip'
		yy = uint8(interp1(y,xx,'pchip'));
	case 'cubic'
		yy = uint8(interp1(y,xx,'cubic'));
	case 'v5cubic'
		yy = uint8(interp1(y,xx,'v5cubic'));
end

author: Sable

repo: mcbench-benchmarks


text:
% Description:
%
%     Enumerate N choose 2 combinations.
%
% Syntax:
%
%     c = mNC2(n)
%
% Inputs:
%
%     n - [ 1 x 1 ] (int)
%
% Outputs:
%
%     c - [ N x 2 ] (int)
%
% Details:
%
% Examples:
%
% Notes:
%
% Author(s):
%
%     William Gruner (williamgruner@gmail.com)
%
% References:
%
% Acknowledgements:
%
%     Many thanks to Dr. Erik Erhardt and Dr. Elena Allen of the Mind Research
%     Network (www.mrn.org) for their continued collaboration.
%
% Version:
%
%     $Author: williamgruner $
%     $Date: 2010-04-01 11:39:09 -0600 (Thu, 01 Apr 2010) $
%     $Revision: 482 $

function c = mNC2(n)

    c = [];

    for i = 1 : n
        for j = i + 1 : n
            c(end + 1, :) = [ i j ];
        end
    end


text:
function op = smooth_logdet(q,C)
% SMOOTH_LOGDET   The -log( det( X ) ) function.
%   (Note the minus sign)
%   FUNC = SMOOTH_LOGDET( q ) returns a function handle that
%   provides a TFOCS-compatible implementation of the funciton
%       -q*log( det( X ) )
%
%   FUNC = SMOOTH_LOGDET( q, C ) represents
%       -q*log( det( X ) ) + < C, X >, where C is symmetric/Hermitian
%
%   X must be symmetric/Hermitian and positive definite,
%   and q must be a positive real number (if not provided,
%   the default value is q = 1).
%
%   N.B. it is the user's responsibility to ensure
%   that X is Hermitian and pos. def., since
%   automatically checking is expensive.
%
%   This function is differentiable, but the gradient
%   is not Lipschitz on the domain X > 0
%
%   This function does support proximity operations, and so
%   it may be used as a nonsmooth function.
%   However, the input must be symmetric positive definite
%   (and if C is used, then it must also be > tC )

% SRB: have not yet tested this.
% SRB: I think we CAN compute the proximity operator to logdet.
%       Will implement this in prox_logdet
if nargin < 1, q = 1; end
if nargin < 2, C = []; end
if ~isreal(q) || q <= 0
    error('First argument must be real and positive');
end

%op = @smooth_logdet_impl;
if isempty(C)
    op = @(varargin)smooth_logdet_impl( q, varargin{:} );
else
    op = @(varargin)smooth_logdet_impl_C( q, C, varargin{:} );
end

function [ v, g ] = smooth_logdet_impl( q, x, t )
if size(x,1) ~= size(x,2)
    error('smooth_logdet: input must be a square matrix');
end
switch nargin
    case 2
        % the function is being used in a "smooth" fashion
        %v = -log(det(x));
        v = -2*q*sum(log(diag(chol(x))));  % chol() takes half the time as det()
                                 % and it is easier to avoid overflow errors
                                 % since we sum the logs.
                                 % Also, chol() will warn if not pos. def.
        if nargout > 1
            g = -q*inv(x);
            % it would be nice to make g a function handle
            % that calculates g(y) = -x\y
        end

    case 3
        % the function is being used in a "nonsmooth" fashion
        % i.e. return g = argmin_g  -q*log(det(g)) + 1/(2t)||g-x||^2
        x = full(x+x')/2;  % March 2015, project it to be symmetric
        [V,D]   = safe_eig(x);
        d       = diag(D);
        % This is OK: input need not be pos def
        %if any(d<=0),
            %v   = Inf;
            %g   = nan(size(x));
            %return;
%%             error('log_det requires a positive definite point');
        %end
        l       = ( d + sqrt( d.^2 + 4*t*q ) )/2;
        g       = V*diag(l)*V';
        v       = -q*sum(log(l));
    otherwise
        error('Wrong number of arguments');
end

function [ v, g ] = smooth_logdet_impl_C( q, C, x, t )
if size(x,1) ~= size(x,2)
    error('smooth_logdet: input must be a square matrix');
end
if size(C,1) ~= size(C,2)
    error('smooth_logdet: input must be a square matrix');
end
switch nargin
    case 3
        % the function is being used in a "smooth" fashion
        %v = -log(det(x));
        v = -2*q*sum(log(diag(chol(x))));  % chol() takes half the time as det()
                                 % and it is easier to avoid overflow errors
                                 % since we sum the logs.
                                 % Also, chol() will warn if not pos. def.
        v = v + tfocs_dot( C, x );
        if nargout > 1
            g = -q*inv(x) + C;
            % it would be nice to make g a function handle
            % that calculates g(y) = -x\y
        end

    case 4
        % the function is being used in a "nonsmooth" fashion
        % i.e. return g = argmin_g  -q*log(det(g)) + 1/(2t)||g-x||^2
        x       = x - t*C;
        x = full(x+x')/2;  % March 2015, project it to be symmetric
        [V,D]   = safe_eig(x);
        d       = diag(D);
        % This is OK: input need not be pos def
        l       = ( d + sqrt( d.^2 + 4*t*q ) )/2;
        g       = V*diag(l)*V';
        v       = -q*sum(log(l));
        v       = v + tfocs_dot( C, g );
    otherwise
        error('Wrong number of arguments');
end

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

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [y]=lars(p,ind,mat,rhs,n,d,smin,h,i0,j0)
%The dimension is p^2 (thus w
k0=p*p;
x=zeros(k0,1);
pp=2.^(0:d-1);
for i=1:k0
   ix=ind(d*(i-1)+1:d*i);
   ix=dot((ix-1),pp);
   x(i)=smin+h*ix;
end
mm=x(1)*mat{1};
for i=2:numel(mat)
  mm=mm+x(i)*mat{i};
end
   sol=mm \ rhs;
   sol=reshape(sol,n,n);
y=sol(i0,j0);
return
end

author: oseledets

repo: TT-Toolbox


text:
function varargout = expression( nm, varargin )

%EXPRESSION Declares a single CVX object for storing subexpressions.
%   EXPRESSION x
%   where x is a valid MATLAB variable nm, declares a scalar expression
%   holder for the current cvx model. Like a variable, an expression holder
%   can be used in constraints and objectives, according to the DCP ruleset.
%   However, unlike a variable, an expression holder is initialized to
%   zero, because the intent is for it to hold intermediate computations.
%
%   EXPRESSION x(n1,n2,...,nk)
%   declares a vector, matrix, or array expression holder with dimensions
%   n1, n2, ..., nk, each of which must be nonnegative integers. The value
%   of the expression holder is initialized to zero.
%
%   Examples:
%      variable x y
%      expression z
%      z = 2 * x - y;
%
%   See also EXPRESSIONS.

global cvx___
prob = evalin( 'caller', 'cvx_problem', '[]' );
if ~isa( prob, 'cvxprob' ),
    error( 'No CVX model exists in this scope.' );
elseif isempty( cvx___.problems ) || cvx___.problems( end ).self ~= prob,
    error( 'Internal CVX data corruption. Please CLEAR ALL and rebuild your model.' );
elseif nargin > 1,
    error( 'Too many input arguments.\nTrying to declare multiple expression holders? Use the EXPRESSIONS keyword instead.', 1 ); %#ok
end

%
% Step 1: separate the name from the parenthetical, verify the name
%

toks = regexp( nm, '^\s*([a-zA-Z]\w*)\s*(\(.*\))?\s*$', 'tokens' );
if isempty( toks ),
    error( 'Invalid variable specification: %s', nm );
end
toks = toks{1};
x.name = toks{1};
x.size = toks{2};
if x.name(end) == '_',
    error( 'Invalid expression specification: %s\n   Names ending in underscores are reserved for internal use.', nm );
end

%
% Step 2: Parse the size. In effect, all we do here is surround what is
% replace the surrounding parentheses with square braces and evaluate. All
% that matters is the result is a valid size vector. In particular, it
% need to be a simple comma-delimited list.
%

if isempty( x.size ),
	x.size = [1,1];
else
    try
        x.size = evalin( 'caller', [ '[', x.size(2:end-1), '];' ] );
    catch exc
        error( exc.identifier, exc.message );
    end
    [ temp, x.size ] = cvx_check_dimlist( x.size, true );
    if ~temp,
        error( 'Invalid expression specification: %s\n   Dimension list must be a vector of finite nonnegative integers.', nm );
    end
end

%
% Step 3. Initialize
%

v = cvx( x.size, [] );
if nargout > 0,
    varargout{1} = v;
else
    assignin( 'caller', x.name, v );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
  function out = padn(mat, newdim)
%|function out = padn(mat, newdim)
%| pad input matrix to newdim, preserving 'center point'
%| using this avoids matlab's padarray.m which is in image proc toolbox

if nargin < 1, ir_usage, end
if streq(mat, 'test'), padn_test, return, end

% default: round to next up power of 2
if nargin < 2
	newdim = 2 .^ ceil(log2(size(mat)))
end

olddim = size(mat);
if any(newdim < olddim), error('must be bigger'), end

idx = cell(1, length(newdim));
for ii=1:length(newdim)
	pad = newdim(ii) - olddim(ii);
	if ~rem(pad,2) % even
		offset = pad/2;
	else
		if rem(olddim(ii),2) % odd
			offset = ceil(pad/2);
		else
			offset = floor(pad/2);
		end
	end
	idx{ii} = [1:olddim(ii)] + offset;
end
out = zeros1(newdim);
out(idx{:}) = mat;

 function z = zeros1(dim)
%function z = zeros1(dim)
% make array of zeros(dims(1), dims(2), ...)
% works logically even if the input is a scalar
% jeff fessler

if nargin < 1, ir_usage, end

if length(dim) == 1
	dim = [dim 1]
end
z = zeros(dim);

function padn_test
x = ones(3,4);
jf_equal(x, unpadn(padn(x, [4 6]), size(x)))
jf_equal(x, unpadn(padn(x, [5 7]), size(x)))

jf_equal(padn([1 2 1], [1 5]), [0 1 2 1 0])
jf_equal(padn([0 1 2 1], [1 5]), [0 1 2 1 0])
jf_equal(padn([0 1 2 1], [1 6]), [0 0 1 2 1 0])
jf_equal(padn([1 2 1], [1 6]), [0 0 1 2 1 0])

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Ex0201: Perform SPHARM-PDM topology fix')

SpharmMatDir = '../';
if ~ispc
    curP = pwd;
    cd(SpharmMatDir);
    SpharmMatDir = [pwd '/'];
    cd(curP);
end
CodeDir = [SpharmMatDir 'code'];
if exist([CodeDir '/sysConfig.mat'])
    load([CodeDir '/sysConfig.mat']);
else
    disp('PDM_path variable (sysConfig.mat) does not exist in SPHARM-MAT directory.')
    disp('Please set the path to SHARM-PDM for SPHARM-MAT and use this script.');
    return;
end
addpath(CodeDir);

DataFolder = 'data';

%% (1) Input data directory
InDataDir = [SpharmMatDir DataFolder '/Ex0202/hip01_bim'];

%% (2) List input file names
inFiles = dir([InDataDir '/*.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [InDataDir '/' inFiles(i).name];
end

%% (3) Display input objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'orig';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'mesh';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'none';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

%% (4) Format Conversion bim2gipl
confs.OutDirectory = [SpharmMatDir DataFolder '/Ex0202/hip02_gipl'];

if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end

SpharmMatUtilFormatConvert(confs, inNames, 'bim2gipl');

%% (5) Perform SPHRM-PDM topology fix
inFiles = dir([confs.OutDirectory '/*.gipl']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [confs.OutDirectory '/' inFiles(i).name];
end

PDM_Fix.command = 'SegPostProcess';
if ispc
    PDM_path = PDM_path_pc;
elseif isunix
    PDM_path = PDM_path_unix;
end

com = [PDM_path '/' PDM_Fix.command];
PDM_Fix.exist = 0;
if ispc
    if exist([com '.exe'],'file')
        PDM_Fix.exist = 1;
    end
else
    if exist(com, 'file')
        PDM_Fix.exist = 1;
    end
end

confs.vars = {'space','OutDirectory', 'others'};
confs.args = [3 200 20];
confs.space = [0.75,0.75,0.75];
confs.OutDirectory = [SpharmMatDir DataFolder '/Ex0202/hip03_gipl_fix'];
confs.others = '';
confs.command = PDM_Fix.command;
confs.path = PDM_path;

if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end

if (PDM_Fix.exist == 0)
    errordlg(sprintf('%s function does not exist in the PDM directory, %s',PDM_Fix.command,PDM_path));
    rmpath(Codedir);
    return;
else
    outNames = SpharmMatUtilTopologyFix(confs, inNames, 'PDM_Fix');
end

%% (6) Format Conversion gipl2bim

confs.OutDirectory = [SpharmMatDir DataFolder '/Ex0202/hip04_fix'];

if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end

SpharmMatUtilFormatConvert(confs, outNames, 'gipl2bim');

%% (7) Display output objects (optional)
inFiles = dir([confs.OutDirectory '/*.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [confs.OutDirectory '/' inFiles(i).name];
end

    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'orig';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'mesh';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'none';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

rmpath(CodeDir);

clear;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% (Internal) Filtered calculation of RR time series from QRS detections
%
% This function calculates a gap-removed version of the RR interval
% sequence.
%
%       RR = RR_calculation(QRS_detections, sampling_rate)
%
% Arguments:
%
%	   QRS_detections: the QRS detections
%
%	   sampling_rate: the sampling frequency of the detections.
%
% Example
%
% See also resample_sequences, MedianFiltSequence
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 25/4/2017
% Last update: 25/4/2017
% Copyright 2008-2017
%
function [ RR, RR_filt] = RR_calculation(QRS_detections, sampling_rate, RR_filt)

    RR = diff(QRS_detections);
    RR = [ RR(1,:); RR ];

    if( nargin < 3 || length(RR_filt) ~= length(QRS_detections) )
        % resample the RR sequence. Heavy in computation
        RR_filt = MedianFiltSequence(QRS_detections, RR, 2*sampling_rate);
    end

    % remove gaps using next heartbeats
    gap_relative_time = 3; % times the median RR interval
    aux_val = RR ./ RR_filt;
    aux_idx = find(aux_val >= gap_relative_time);

    for ii = rowvec(aux_idx)

        aux_idx = find( QRS_detections > QRS_detections(ii) & aux_val < gap_relative_time, 1, 'first');

        if ( isempty(aux_idx) )
            aux_idx = find( QRS_detections < QRS_detections(ii) & aux_val < gap_relative_time, 1, 'last');
        end

        RR(ii) = RR( aux_idx );

    end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [C, grps, G] = grp2cell(X, G, dim)
%GRP2CELL Splits an array into a cell using a grouping variable.
% Usage:
%   C = grp2cell(X, G)
%   C = grp2cell(X, G, dim)
%   [C, grps] = grp2cell(_)
%
% Args:
%   X: array with at least one dimension of the same size as G
%   G: grouping variable vector
%   dim: dimension along which to slice X (default: [])
%        if empty, first dimension of the same size as G is used
%
% Returns:
%   C: cell array with as many cells as unique elements in G
%   grps: values in G that correspond to the grouping in C
%
% Example:
%   >> C = grp2cell(X,G,dim);
%   >> isequal(cellcat(C, dim), X)
%   ans =
%     logical
%      1
%
% See also: arr2cell, cellcat, celljoin

if nargin < 3; dim = []; end
if isscalar(G)
    % Default to first non-singleton dimension
    dim = find(size(X) > 1, 1);
    if isempty(dim); dim = 1; end % or 1

    % Create grouping vector
    G = ceil((1:size(X,dim)) / G);
end

if isempty(dim)
    dim = find(size(X) == numel(G),1);
end

subs = af(@(x) 1:x, size(X));

grps = unique(G);
C = cell1(numel(grps),dim);
for i = 1:numel(grps)
    subs_i = subs;
    subs_i{dim} = find(G == grps(i));

    C{i} = X(subs_i{:});
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
% pop_icathresh() - main menu for choosing threshold for component
%                      rejection in EEGLAB.
%
% Usage:
%   >> [OUTEEG rej] = pop_icathresh(INEEG, threshval, rejmethod,
%                                      rejvalue, interact);
%
% Inputs:
%   INEEG      - input dataset
%   threshval  - values of thresholds for each of the 3 statistical
%                measures. Default is [] and the program uses the value
%                in the dataset.
%   rejmethod  - either 'percent', 'dataset' or 'current'. 'percent'
%                will reject a given percentage of components with
%                the highest value in one or several statistical
%                measure. 'dataset' will use an other dataset for
%                calibration. 'current' will use the current dataset
%                for calibration. Default is 'current'.
%   rejvalue   - percentage if rejmethod is 'percent', dataset number
%                if rejmethod is 'dataset' (no input if 'current'). If it
%                is a percentage, '25' for instance means that the 25%
%                independent components with the highest values of one
%                statistical measure are  rejected. Note that, one can
%                also enter one percentage per statistical value (such as
%                25 20 30).
%   interact   - interactive windows or just rejection
%
% Inputs:
%   OUTEEG     - output dataset with updated thresholds
%   rej        - new rejection array
%
% Graphic interface:
%   The graphic interface is divided into 3 parts. On the top, the
% experimenter chooses the basis for the component rejection (see
% rejmethod input). On the middle panel, the user can tune thresholds
% manually and plot the distribution of statistical values. Each time
% that setting of one of these two top panels is modified, the curve
% on the third panel are redrawn.
%   The third panel is composed of 3 graphs, one for each statistical
% measure. The blue curve on each graph indicates the accuracy of the
% measure to detect artifactual components and non-artifactual
% components of a dataset. Note that 'artifactual components are
% defined with the rejection methods fields on the top (if you use
% the percentage methods, these artifactual components may not actually
% be artifactual). For accurate rejection, one should first reject
% artifactual component manually and then set up the threshold to
% approximate this manual rejection. The green curve indicate the
% rejection when all measure are considered. Points on the blue and
% on the green curves indicate the position of the current threshold on
% the parametric curve. Not that for the green cumulative curve, this
% point is at the same location for all graphs.
%   How to tune the threshold? To tune the threshold, one should try to
% maximize the detection of non-artifactual components (because it is
% preferable to miss some artifactual components than to classify as
% artifactual non-artifact components).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eeglab()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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

%PROBLEM: when the % is set and we change manually yhe threshold,
% the software comes back to the percentage rejection

% 01-25-02 reformated help & license -ad

function [EEG, rej, com] = pop_icathresh( EEG, threshval, rejmethod, rejvalue, interact);

com = [];
rej = EEG.reject.gcompreject;
if nargin < 1
	help pop_icathresh;
	return;
end;

if nargin < 2
	threshval = [];
end;
if nargin < 3
	rejmethod = 'current';
end;
if nargin < 4
	rejvalue = 25;
end;
if nargin < 5
	interact = 1;
end;
if ~isempty(threshval)
	EEG.reject.threshentropy = threshval(1);
	EEG.reject.threshkurtact = threshval(2);
	EEG.reject.threshkurtdist = threshval(3);
end;

tagmenu = 'pop_icathresh';

if ~isempty( findobj('tag', tagmenu))
	error('cannot open two identical windows, close the first one first');
end;

% the current rejection will be stored in userdata of the figure
% --------------------------------------------------------------
gcf = figure('visible', 'off', 'numbertitle', 'off', 'name', 'Choose thresholds', 'tag', tagmenu, 'userdata', rej);
pos = get(gca, 'position');
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)]./100;
axis off;

% definition of callbacks
% -----------------------

cb_cancel  = [  'userdat = get(gcbo, ''userdata'');' ... % restore thresholds
				'EEG.reject.threshentropy  = userdat{1};' ...
				'EEG.reject.threshkurtact  = userdat{2};' ...
				'EEG.reject.threshkurtdist = userdat{3};' ...
				'clear userdat; close(gcbf);' ];

drawgraphs  = [ 'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'if get( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'') == 1, ' ... % test if percentage
				'   perc = str2num(get( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''String''));' ...
				'   if length(perc < 2), perc = [perc perc perc]; end;' ...
				'   perc = round((100-perc) * length( EEG.stats.compenta) / 100);' ... % convert the percentage
				'   method = zeros( size( EEG.stats.compenta ) );' ...
				'   [tmprej tmpindex] = sort(EEG.stats.compenta(:));' ...
				'   method( tmpindex(perc(1)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''entstring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold
				'   [tmprej tmpindex] = sort(EEG.stats.compkurta(:));' ...
				'   method( tmpindex(perc(2)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''kurtstring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold
				'   [tmprej tmpindex] = sort(EEG.stats.compkurtdist(:));' ...
				'   method( tmpindex(perc(3)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''kurtdiststring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold
				'   allvalues = [EEG.stats.compenta(:) EEG.stats.compkurta(:) EEG.stats.compkurtdist(:) ];' ...
				'   clear perc tmprej tmpindex;' ...
                'end;' ...
				'if get( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'') == 1, ' ... % test if other dataset
				'   di = str2num(get( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''String''));' ...
				'   if isempty( di ), clear fig; return; end;' ...
				'   method    = ALLEEG(di).reject.gcompreject'';' ...
				'   allvalues = [ALLEEG(di).stats.compenta(:) ALLEEG(di).stats.compkurta(:) ALLEEG(di).stats.compkurtdist(:) ];' ...
				'   clear di;' ...
                'end;' ...
				'if get( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'') == 1, ' ... % test if current dataset
				'   method = EEG.reject.gcompreject'';' ...
				'   allvalues = [EEG.stats.compenta(:) EEG.stats.compkurta(:) EEG.stats.compkurtdist(:) ];' ...
                'end;' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphent'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [0 EEG.reject.threshkurtact EEG.reject.threshkurtdist ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshentropy, ''Entropy of activity'',''Artifact detection (%)'', ''Non-artifact detection (%)'');' ...
				'set(gca, ''tag'', ''graphent'');' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphkurt'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [EEG.reject.threshentropy 0 EEG.reject.threshkurtdist ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshkurtact, ''Kurtosis of activity'', ''Artifact detection (%)'', '''');' ...
				'set(gca, ''tag'', ''graphkurt'');' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphkurtdist'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [EEG.reject.threshentropy EEG.reject.threshkurtact 0 ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshkurtdist, ''Kurtosis of topography'', ''Artifact detection (%)'', '''');' ...
				'set(gca, ''tag'', ''graphkurtdist'');' ...
				'clear method allvalues fig;' ...
				];

cb_percent = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''on'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''off''); clear fig;' drawgraphs ];

cb_other   = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''off'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''on''); clear fig;' drawgraphs ];

cb_current = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''off'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''off''); clear fig;' drawgraphs ];

cb_entropy  = [ 'EEG.reject.threshentropy = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];
cb_kurtact  = [ 'EEG.reject.threshkurtact = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];
cb_kurtdist  = [ 'EEG.reject.threshkurtdist = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];

if interact
	cb_calrej = [ 'ButtonName=questdlg2( ''This will erase previous projections'', ''Confirmation'', ''CANCEL'', ''OK'', ''OK'');' ]
else
	cb_calrej = [ 'ButtonName= ''OK''' ];
end;
cb_calrej = [ cb_calrej ...
    		  'switch ButtonName,' ...
       	 	  '   case ''OK'',' ...
			  '    rej1 = find(EEG.stats.compenta > EEG.reject.threshentropy);' ...
			  '    rej2 = find(EEG.stats.compkurta > EEG.reject.threshkurtact);' ...
			  '    rej3 = find(EEG.stats.compkurtdist > EEG.reject.threshkurtdist);' ...
			  '    EEG.reject.gcompreject = (rej1 & rej2) | rej3;' ...
			  '    clear rej1 rej2 rej3;' ...
			  'end; clear ButtonName;' ];

% default value for rejection methods
% -----------------------------------
rejvalother = '';
rejvalpercent = '25';
switch rejmethod
	case 'percent', rejvalpercent = num2str( rejvalue);
	case 'dataset', rejvalother = num2str( rejvalue);
end;

% -----------------------------------------------------
allh = supergui(gcf, { 	[1] ...
					[2 1] [2 1] [2 1] ...
					[1] ...
					[1] ...
					[1.5 0.5 1] [1.5 0.5 1] [1.5 0.5 1] ...
					[1] [1] [1] [1] [1] [1] [1] [1] [1] [1] ...
					[1 1 1 1 1] ...
					}, [], ...
	{ 'Style', 'text', 'string', 'Calibration method', 'FontSize', 13, 'fontweight', 'bold' }, ...
	...
	{ 'style', 'checkbox', 'String', '% of artifactual components (can also put one % per rejection)', 'tag', 'Ipercent', 'value', 1, 'callback', cb_percent}, ...
	{ 'style', 'edit', 'String', rejvalpercent, 'tag', 'Ipercenttext', 'callback', drawgraphs }, ...
	...
	{ 'style', 'checkbox', 'String', 'Specific dataset (enter dataset number)', 'tag', 'Iother', 'value', 0, 'callback', cb_other}, ...
	{ 'style', 'edit', 'String', rejvalother, 'tag', 'Iothertext', 'enable', 'off', 'callback', drawgraphs }, ...
	...
	{ 'style', 'checkbox', 'String', 'Current dataset', 'tag', 'Icurrent', 'value', 0, 'callback', cb_current}, ...
	{ }, ...
	...
	{ }, ...
	...
	{ 'Style', 'text', 'string', 'Threshold values', 'FontSize', 13, 'fontweight', 'bold' }, ...
	...
	{ 'style', 'text', 'String', 'Entropy of activity threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshentropy), 'tag', 'entstring', 'callback',  cb_entropy }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compenta,20); title(''Entropy of activity'');' }, ...
	...
	{ 'style', 'text', 'String', 'Kurtosis of activity threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshkurtact), 'tag', 'kurtstring', 'callback', cb_kurtact }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compkurta,200); title(''Kurtosis of activity'');' },  ...
	...
	{ 'style', 'text', 'String', 'Kurtosis of topography threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshkurtdist), 'tag', 'kurtdiststring', 'callback', cb_kurtdist }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compkurtdist,20); title(''Kurtosis of topography'');' }, ...
	...
	{ }, { }, { }, { }, { }, { }, { }, { }, { }, { }, ...
	...
	{ 'style', 'pushbutton', 'String', 'Cancel', 'callback', cb_cancel, 'userdata', { EEG.reject.threshentropy EEG.reject.threshkurtact EEG.reject.threshkurtdist }}, ...
	{ 'style', 'pushbutton', 'String', 'Auto thresh'  , 'callback', '', 'enable', 'off' }, ...
	{ 'style', 'pushbutton', 'String', 'Help'  , 'callback', 'pophelp(''pop_icathresh'');' }, ...
	{ 'style', 'pushbutton', 'String', 'Accept thresh.' , 'callback', 'close(gcbf);' }, ...
	{ 'style', 'pushbutton', 'String', 'Calc. rejection' , 'callback', 'close(gcbf);' } ...
	...
	);

h = axes('position', [0  15 30 30].*s+q, 'tag', 'graphent');
h = axes('position', [35 15 30 30].*s+q, 'tag', 'graphkurt');
h = axes('position', [70 15 30 30].*s+q, 'tag', 'graphkurtdist');

rejmethod
switch rejmethod
	case 'dataset', eval([ 'gcbf = [];' cb_other]);
	case 'current', eval([ 'gcbf = [];' cb_current]);
	otherwise, eval([ 'gcbf = [];' cb_percent]);
end;

return;

author: PatternRecognition

repo: OpenBMI


text:
classdef DualUpdater_NullSpace < handle

    properties (Access = private)
       constraint
       cost
       designVariable
       dualVariable
       problem
       options
       constraintCase
       nConstr
       dualOld
       isInequality
       position
       parameter
    end

    methods (Access = public)

        function obj = DualUpdater_NullSpace(cParams)
            obj.init(cParams);
            obj.defineConstraintCases();
            obj.defineRangeStepParameterValue(cParams);
        end

        function defineConstraintCases(obj)
            obj.isInequality = false;
            k                = 1;
            for i = 1:obj.nConstr
                switch obj.constraintCase{i}
                    case {'EQUALITY'}

                    case {'INEQUALITY'}
                        obj.isInequality = true;
                        obj.position(k)  = i;
                        k                = k + 1;
                end
            end
        end

        function defineRangeStepParameterValue(obj,cParams)
            switch cParams.optimizerNames.primal
                case {'PROJECTED GRADIENT','HAMILTON JACOBI'}
                    obj.parameter = inf;
                case {'SLERP'}
                    obj.parameter = 5;
                otherwise

            end
        end

        function update(obj)
            if obj.isInequality
                obj.computeQuadraticProblem();
            else
                obj.computeDirectDual();
            end
        end

        function reset(obj)
            obj.dualVariable.value = obj.dualOld;
        end

        function updateOld(obj)
            obj.dualOld = obj.dualVariable.value;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.cost           = cParams.cost;
            obj.constraint     = cParams.constraint;
            obj.designVariable = cParams.designVar;
            obj.constraintCase = cParams.constraintCase;
            obj.dualVariable   = cParams.dualVariable;
            obj.dualOld        = obj.dualVariable.value;
            obj.nConstr        = cParams.constraint.nSF;
        end

        function computeDirectDual(obj)
            obj.constraint.computeFunctionAndGradient();
            obj.cost.computeFunctionAndGradient();
            DJ = obj.cost.gradient;
            Dh = obj.constraint.gradient;
            h  = obj.constraint.value;
            S  = (Dh'*Dh)^-1;
            aJ = 1;
            aC = 1;
            f  = obj.parameter;
            f2 = 1;
            AC = min(f,f2*aC/aJ*S*h);
            AJ = -aC/aJ*S*Dh'*DJ;
            l  = AC + AJ;
            obj.dualVariable.value = l;
        end

        function computeQuadraticProblem(obj)
            obj.constraint.computeFunctionAndGradient();
            obj.cost.computeFunctionAndGradient();
            obj.computeDualProblemParameters();
            obj.computeDualProblemOptions();
            PROBLEM         = obj.problem;
            PROBLEM.options = obj.options;
            l = quadprog(PROBLEM);
            obj.dualVariable.value = l;
        end

         function computeDualProblemParameters(obj)
            Dg          = obj.constraint.gradient;
            DJ          = obj.cost.gradient;
            g           = obj.constraint.value;
            f           = obj.parameter;
            c           = min(f*(Dg'*Dg),g);
            i           = length(g);
            p           = obj.position;
            prob.H      = Dg'*Dg;
            prob.f      = -c + Dg'*DJ;
            prob.A      = [];
            prob.b      = [];
            prob.Aeq    = [];
            prob.beq    = [];
            prob.lb     = -inf*ones(i,1);
            prob.lb(p)  = 0;
            prob.ub     = inf*ones(length(g),1);
            prob.x0     = zeros(length(prob.H),1);
            prob.solver = 'quadprog';
            obj.problem = prob;
        end

        function computeDualProblemOptions(obj)
            opts = optimoptions("quadprog");
            opts = struct( ...
                'Algorithm','interior-point-convex', ...
                'Diagnostics','off', ...
                'Display','none', ...
                'HessMult',[], ...
                'MaxIter',1e3, ...
                'MaxPCGIter','max(1,floor(numberOfVariables/2))', ...
                'PrecondBandWidth',0, ...
                'ProblemdefOptions', struct, ...
                'TolCon',1e-5, ...
                'TolFun',[], ...
                'TolFunValue', [], ...
                'TolPCG',0.1, ...
                'TolX',100*eps, ...
                'TypicalX','ones(numberOfVariables,1)', ...
                'LinearSolver', 'auto', ...
                'ObjectiveLimit', -1e20 ...
                );
            obj.options = opts;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [node,elem,face]=volmap2mesh(img,ix,iy,iz,elemnum,maxvol,thickness,Amat,Bvec)
%
% [node,elem,face]=volmap2mesh(img,ix,iy,iz,thickness,elemnum,maxvol,A,B)
%
% convert a binary volume to tetrahedral mesh followed by an Affine transform
%
% author: Qianqian Fang (q.fang at neu.edu)
% date:   2008/01/12
%
% input:
%        img, ix,iy,iz, elemnum and  maxvol: see vol2mesh.m
%        thickness: scale z-dimension of the mesh to specified thickness,
%                   if thickness==0, scaling is bypassed
%        Amat: a 3x3 transformation matrix
%        Bvec: a 3x1 vector
%        Amat and Bvec maps the image index space to real world coordnate system by
%                   [x,y,z]_new=Amat*[x,y,z]_old+Bvec
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

[node,elem,face]=vol2mesh(img,ix,iy,iz,elemnum,maxvol);

node(:,1:3)=(Amat*node(:,1:3)'+repmat(Bvec(:),1,size(node,1)))';

if(thickness)
	zmin=min(node(:,3));
	zmax=max(node(:,3));
	node(:,3)=(node(:,3)-zmin)/(zmax-zmin)*thickness;
end

author: fangq


text:
function varargout = process_canoltymap( varargin )
% This function generates Canolty like maps (Science 2006, figure 1) for the input signal.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Canolty maps';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Frequency';
    sProcess.Index       = 661;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw',      'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Resting#Canolty_maps';
    % ==== INPUT ====
    sProcess.options.label_in.Comment = 'Input options:';
    sProcess.options.label_in.Type    = 'label';
    % === TIME WINDOW
    sProcess.options.timewindow.Comment = 'Time window: ';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
    % === SENSOR SELECTION
    sProcess.options.target_data.Comment    = 'Sensor types or names (empty=all): ';
    sProcess.options.target_data.Type       = 'text';
    sProcess.options.target_data.Value      = 'MEG, EEG';
    sProcess.options.target_data.InputTypes = {'data', 'raw'};
    % === SCOUTS SELECTION
    sProcess.options.scouts.Comment    = 'Use scouts';
    sProcess.options.scouts.Type       = 'scout_confirm';
    sProcess.options.scouts.Value      = {};
    sProcess.options.scouts.InputTypes = {'results'};
    % === SCOUT FUNCTION ===
    sProcess.options.scoutfunc.Comment    = {'Mean', 'Max', 'PCA', 'Std', 'All', 'Scout function:'};
    sProcess.options.scoutfunc.Type       = 'radio_line';
    sProcess.options.scoutfunc.Value      = 1;
    sProcess.options.scoutfunc.InputTypes = {'results'};
    % === SCOUT TIME ===
    sProcess.options.scouttime.Comment    = {'Before', 'After', 'When to apply the scout function:'};
    sProcess.options.scouttime.Type       = 'radio_line';
    sProcess.options.scouttime.Value      = 1;
    sProcess.options.scouttime.InputTypes = {'results'};
    % === ROW NAMES
    sProcess.options.target_tf.Comment    = 'Row names or indices (empty=all): ';
    sProcess.options.target_tf.Type       = 'text';
    sProcess.options.target_tf.Value      = '';
    sProcess.options.target_tf.InputTypes = {'timefreq', 'matrix'};

    % ==== ESTIMATOR ====
    sProcess.options.label_method.Comment = '
Estimator options:';
    sProcess.options.label_method.Type    = 'label';
    % === EPOCH TIME
    sProcess.options.epochtime.Comment = 'Epoch time: ';
    sProcess.options.epochtime.Type    = 'range';
    sProcess.options.epochtime.Value   = {[-1, 1], 'ms', []};
    % === NESTING FREQ
    sProcess.options.lowfreq.Comment = 'Nesting frequency (low): ';
    sProcess.options.lowfreq.Type    = 'value';
    sProcess.options.lowfreq.Value   = {4, 'Hz', 2};
    % === MAX_BLOCK_SIZE
    sProcess.options.max_block_size.Comment = 'Number of signals to process at once: ';
    sProcess.options.max_block_size.Type    = 'value';
    sProcess.options.max_block_size.Value   = {100, ' ', 0};

    % ==== OUTPUT ====
    sProcess.options.label_out.Comment = '
Output configuration:';
    sProcess.options.label_out.Type    = 'label';
    % === SAVE AVERAGED LOW-FREQ SIGNALS
    sProcess.options.save_erp.Comment = 'Save averaged low frequency signals';
    sProcess.options.save_erp.Type    = 'checkbox';
    sProcess.options.save_erp.Value   = 1;
end

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

%% ===== RUN =====
function OutputFile = Run(sProcess, sInput) %#ok<DEFNU>
    % Get options
    if ~isempty(sProcess.options.timewindow.Value) && iscell(sProcess.options.timewindow.Value)
        OPTIONS.TimeWindow = sProcess.options.timewindow.Value{1};
    else
        OPTIONS.TimeWindow = [];
    end
    OPTIONS.EpochTime  = sProcess.options.epochtime.Value{1};
    OPTIONS.LowFreq    = sProcess.options.lowfreq.Value{1};
    OPTIONS.MaxSignals = sProcess.options.max_block_size.Value{1};
    OPTIONS.SaveErp    = sProcess.options.save_erp.Value;
    % Get target
    if ismember(sInput.FileType, {'data','raw'}) && isfield(sProcess.options, 'target_data') && ~isempty(sProcess.options.target_data.Value)
        OPTIONS.Target = sProcess.options.target_data.Value;
    elseif ismember(sInput.FileType, {'timefreq', 'matrix'}) && isfield(sProcess.options, 'target_tf') && ~isempty(sProcess.options.target_tf.Value)
        OPTIONS.Target = sProcess.options.target_tf.Value;
    else
        OPTIONS.Target = [];
    end
    OutputFile = {};

    % ===== GET SCOUTS OPTIONS =====
    if strcmpi(sInput.FileType, 'results') && isfield(sProcess.options, 'scouts') && isfield(sProcess.options.scouts, 'Value')
        % Override scouts function
        switch (sProcess.options.scoutfunc.Value)
            case 1, OPTIONS.ScoutFunc = 'mean';
            case 2, OPTIONS.ScoutFunc = 'max';
            case 3, OPTIONS.ScoutFunc = 'pca';
            case 4, OPTIONS.ScoutFunc = 'std';
            case 5, OPTIONS.ScoutFunc = 'all';
        end
        % Scout function order
        switch (sProcess.options.scouttime.Value)
            case 1, OPTIONS.ScoutTime = 'before';
            case 2, OPTIONS.ScoutTime = 'after';
        end
        % Perform some checks
        if strcmpi(OPTIONS.ScoutTime, 'before') && ismember(OPTIONS.ScoutFunc, {'max', 'std'})
            bst_report('Error', sProcess, [], 'Scout functions MAX and STD should not be applied before estimating the PAC.');
            return;
        end
        if strcmpi(OPTIONS.ScoutTime, 'after') && strcmpi(OPTIONS.ScoutFunc, 'pca')
            bst_report('Error', sProcess, [], 'Scout function PCA cannot be applied after estimating the PAC.');
            return;
        end
        % Selected scouts
        AtlasList = sProcess.options.scouts.Value;
        % Set input/output scouts functions
        if ~isempty(AtlasList)
            OPTIONS.Target = AtlasList;
            % Apply function before: get all the scouts time series in advance
            if strcmpi(OPTIONS.ScoutTime, 'before')
                LoadOptions.TargetFunc = OPTIONS.ScoutFunc;
            % Apply function after: Get all the time series of all the scouts
            elseif strcmpi(OPTIONS.ScoutTime, 'after')
                LoadOptions.TargetFunc = 'all';
            end
        end
    end

    % === READ INPUT DATA ===
    % Options for LoadInputFile()
    LoadOptions.LoadFull    = 0;  % Load kernel-based results as kernel+data
    LoadOptions.IgnoreBad   = 0;  % Ignore the bad channels
    LoadOptions.ProcessName = func2str(sProcess.Function);
    % Load input signals
    [sMat, nSignals, iRows] = bst_process('LoadInputFile', sInput.FileName, OPTIONS.Target, OPTIONS.TimeWindow, LoadOptions);
    if isempty(sMat.Data)
        return;
    end
    % Default time window
    if isempty(OPTIONS.TimeWindow)
        OPTIONS.TimeWindow = [sMat.Time(1), sMat.Time(end)];
    end
    % Get sampling frequency
    sRate = 1 / (sMat.Time(2) - sMat.Time(1));
    % Low frequency comment
    if (length(OPTIONS.LowFreq) == 1)
        strFreq = [num2str(round(OPTIONS.LowFreq * 100) / 100), 'Hz'];
    else
        strFreq = '';
    end
    % Replicate low frequency if only one was provided for all the signals
    if (length(OPTIONS.LowFreq) == 1) && (nSignals > 1)
        OPTIONS.LowFreq = repmat(OPTIONS.LowFreq, nSignals, 1);
    % Select only a subset of the low frequencies
    elseif ~isempty(OPTIONS.Target) && ~isempty(iRows) && (length(OPTIONS.LowFreq) ~= length(iRows)) && (max(iRows) < length(OPTIONS.LowFreq))
        OPTIONS.LowFreq = OPTIONS.LowFreq(iRows);
    % Check for compatible size of LowFreq array
    elseif (length(OPTIONS.LowFreq) ~= nSignals)
        bst_report('Error', sProcess, sInput, sprintf('The size of the low-frequency array (%d) does not match the number of signals to process (%d).', length(OPTIONS.LowFreq), nSignals));
        return;
    end

    % ===== CALCULATE CANOLTY MAPS =====
    % Number of blocks of signals
    MAX_BLOCK_SIZE = OPTIONS.MaxSignals;
    nBlocks = ceil(nSignals / MAX_BLOCK_SIZE);
    TF  = [];
    ERP = [];
    % Process each block of signals
    for iBlock = 1:nBlocks
        bst_progress('inc', round(iBlock/nBlocks*100));
        % Indices of the signals
        iSignals = (iBlock-1)*MAX_BLOCK_SIZE+1 : min(iBlock*MAX_BLOCK_SIZE, nSignals);
        % Get signals to process
        if ~isempty(sMat.ImagingKernel)
            Fblock = sMat.ImagingKernel(iSignals,:) * sMat.Data;
        else
            Fblock = sMat.Data(iSignals,:);
        end
        % Calculate canolty map signal
        [TFblock, ERPblock, TimeVectorOut, Freqs, errMsg] = Compute(Fblock, sRate, OPTIONS.LowFreq(iSignals), OPTIONS.EpochTime);
        if ~isempty(errMsg)
            bst_report('Error', sProcess, sInput, errMsg);
            return;
        end
        % Initialize output variable
        if isempty(TF)
            TF  = zeros(nSignals, size(TFblock,2), size(TFblock,3));
            ERP = zeros(nSignals, size(TFblock,2));
        end
        % Copy block results to output variable
        TF(iSignals,:,:) = TFblock;
        ERP(iSignals,:)  = ERPblock;
    end

    % ===== UNCONSTRAINED SOURCES =====
    % Unconstrained sources => SUM for each point
    if ismember(sMat.DataType, {'results','scout','matrix'}) && ~isempty(sMat.nComponents) && (sMat.nComponents ~= 1)
        [TF,  tmp, tmp]                      = bst_source_orient([], sMat.nComponents, sMat.GridAtlas, TF,  'sum',  sMat.DataType, sMat.RowNames);
        [ERP, sMat.GridAtlas, sMat.RowNames] = bst_source_orient([], sMat.nComponents, sMat.GridAtlas, ERP, 'mean', sMat.DataType, sMat.RowNames);
    end

    % ===== PROCESS SCOUTS =====
    % Get scouts
    isScout = ~isempty(OPTIONS.Target) && (isstruct(OPTIONS.Target) || iscell(OPTIONS.Target)) && isfield(sMat, 'Atlas') && isfield(sMat.Atlas, 'Scouts') && ~isempty(sMat.Atlas.Scouts);
    if isScout
        sScouts = sMat.Atlas.Scouts;
    end
    % If the scout function has to be applied AFTER the PAC computation
    if isScout && strcmpi(OPTIONS.ScoutTime, 'after') && ~strcmpi(OPTIONS.ScoutFunc, 'all')
        nScouts      = length(sScouts);
        TF_scouts    = zeros(nScouts, size(TF,2), size(TF,3));
        ERP_scouts   = zeros(nScouts, size(ERP,2));
        % Vertices
        iVerticesAll = [1, cumsum(cellfun(@length, {sScouts.Vertices})) + 1];
        % For each unique row name: compute a measure over the clusters values
        for iScout = 1:nScouts
            iScoutVert = iVerticesAll(iScout):iVerticesAll(iScout+1)-1;
            % Process scouts for the TF maps
            F = reshape(TF(iScoutVert,:,:), length(iScoutVert), []);
            F = bst_scout_value(F, OPTIONS.ScoutFunc);
            TF_scouts(iScout,:,:) = reshape(F, [1, size(TF,2), size(TF,3)]);
            % Process scouts for the ERP map
            ERP_scouts(iScout,:) = bst_scout_value(ERP(iScoutVert,:), OPTIONS.ScoutFunc);
        end
        % Save only the requested rows
        sMat.RowNames = {sScouts.Label}';
        TF  = TF_scouts;
        ERP = ERP_scouts;
    end

    % ===== COMMENT =====
    Comment = 'Canolty maps';
    % Time window (RAW only)
    if ~isempty(strfind(sInput.Condition, '@raw'))
        if ~isempty(strFreq)
            Comment = [Comment, sprintf('(%ds-%ds,%s)', round(OPTIONS.TimeWindow), strFreq)];
        else
            Comment = [Comment, sprintf('(%ds-%ds)', round(OPTIONS.TimeWindow))];
        end
    elseif ~isempty(strFreq)
        Comment = [Comment '(' strFreq ')'];
    end
    % Scouts
    if isScout && (length(sScouts) < 6)
        Comment = [Comment, ':'];
        for is = 1:length(sScouts)
            Comment = [Comment, ' ', sScouts(is).Label];
        end
        Comment = [Comment, ', ', OPTIONS.ScoutFunc];
        if ~strcmpi(OPTIONS.ScoutFunc, 'All')
             Comment = [Comment, ' ' OPTIONS.ScoutTime];
        end
    % Single input
    elseif (length(sMat.RowNames) == 1)
        if iscell(sMat.RowNames)
            Comment = [Comment, ': ' sMat.RowNames{1}];
        else
            Comment = [Comment, ': #', num2str(sMat.RowNames(1))];
        end
    end

    % ===== SAVE TF FILE =====
    % Get the study filename
    sStudy = bst_get('Study', sInput.iStudy);
    % Convert row indices to row names
    if iscell(sMat.RowNames)
        Description = sMat.RowNames;
    else
        Description = cell(length(sMat.RowNames), 1);
        for iRow = 1:length(sMat.RowNames)
            Description{iRow} = num2str(sMat.RowNames(iRow));
        end
    end
    % Create new output structure
    sOutput = db_template('timefreqmat');
    sOutput.TF          = TF;
    sOutput.CanoltyERP  = ERP;
    sOutput.Comment     = Comment;
    sOutput.DataType    = sMat.DataType;
    sOutput.RowNames    = sMat.RowNames;
    sOutput.Time        = TimeVectorOut;
    sOutput.Freqs       = round(Freqs .* 100) / 100;
    sOutput.Measure     = 'other';
    sOutput.Method      = 'canolty';
    sOutput.DataFile    = sInput.FileName;
    sOutput.SurfaceFile = sMat.SurfaceFile;
    sOutput.Atlas       = [];
    sOutput.Options     = OPTIONS;
    sOutput.CanoltyERP  = ERP;
    if ~isempty(sMat.GridLoc)
        sOutput.GridLoc = sMat.GridLoc;
    end
    if ~isempty(sMat.GridAtlas)
        sOutput.GridAtlas = sMat.GridAtlas;
    end
    % Output filename
    OutputFile = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), 'timefreq_canolty');
    % Save file
    bst_save(OutputFile, sOutput, 'v6');
    % Add file to database structure
    db_add_data(sInput.iStudy, OutputFile, sOutput);

    % ===== SAVE MATRIX FILE =====
    if OPTIONS.SaveErp
        % Comment
        ErpComment = sMat.Comment;
        iBar = find(ErpComment == '|', 1, 'last');
        if ~isempty(iBar)
            ErpComment = ErpComment(1:iBar-1);
        end
        ErpComment = [ErpComment ' | ' strrep(Comment, 'Canolty maps', 'Canolty ERP')];
        % Create new output structure
        sOutput = db_template('matrixmat');
        sOutput.Value       = ERP;
        sOutput.Comment     = ErpComment;
        sOutput.Time        = TimeVectorOut;
        sOutput.Description = Description;
        sOutput.Options     = OPTIONS;
        % Output filename
        OutputFileERP = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), 'matrix_canolty');
        % Save file
        bst_save(OutputFileERP, sOutput, 'v6');
        % Add file to database structure
        db_add_data(sInput.iStudy, OutputFileERP, sOutput);
    end
end

%% ===== COMPUTE CANOLTY MAPS =====
function [TF, ERP, TimeOut, chirp_center_high, errMsg] = Compute( F, sRate, lowfreq, EpochTime)
    nTime = size(F,2);
    errMsg = '';

    % ===== CREATE CHIRPLETS =====
    % Definitions
    fmin = 1;
    fmax = 250;
    numfreqs = 70;
    fstep = 0.75;
    % Calculate center frequencies
    temp1 = (0:numfreqs-1) * fstep;
    temp2 = logspace(log10(fmin), log10(fmax), numfreqs);
    temp2 = (temp2-temp2(1)) * ((temp2(end)-temp1(end)) / temp2(end)) + temp2(1);
    chirp_center_high = temp1 + temp2;
    % Delete the frequencies that are too high for the sampling frequency
    chirp_center_high(chirp_center_high >= sRate/3) = [];
    % Calculate chirplets
    [chirpF_high, Freqs] = bst_chirplet(sRate, nTime, chirp_center_high);

    % ===== INITIALIZE RETURNED VARIABLES =====
    % Generate epoch indices
    iEpochTime = round(EpochTime(1)*sRate):round(EpochTime(2)*sRate);
    if isempty(iEpochTime)
        errMsg = 'Invalid epoch time';
    end
    % Output time vector
    TimeOut = iEpochTime / sRate;
    % Initialize returned variable
    TF  = zeros(size(F,1), length(TimeOut), length(chirp_center_high));
    ERP = zeros(size(F,1), length(TimeOut));

    % ===== FFT OF SIGNALS =====
    % Transform sensor time series into analytic signals
    F_fft = fft(F, length(Freqs), 2);
    % This step scales analytic signal such that: real(analytic_signal) = raw_signal
    % but note that analytic signal energy is double that of raw signal energy
    F_fft(:,Freqs<0) = 0;
    F_fft(:,Freqs>0) = 2 * F_fft(:,Freqs>0);

    % ===== LOOP ON SIGNALS =====
    for iSource = 1:size(F,1)
        % === DETECT MIN/MAX LOW FREQ ===
        % Calculate one chirplet for the low frequency
        [chirpF_low, Freqs] = bst_chirplet(sRate, nTime, lowfreq(iSource));
        % Filter again: Positive version of the signal
        fs_low = bst_freqfilter(F(iSource,:), chirpF_low, Freqs);
        % Detection of phase maxima of theta filtered signal (POSITIVE)
        [tmp, iMaxTheta] = find_maxima(angle(fs_low));

        % === FILTER GAMMA ===
        % Filter source signal using all the low-frequency chirplet
        fs_high = bst_freqfilter(F(iSource,:), chirpF_high, Freqs, F_fft(iSource,:));
        % Magnitude
        fs_high = abs(fs_high);
        % Zscore normalization
        fs_high = process_baseline_norm('Compute', fs_high, fs_high, 'zscore');

        % === EPOCH ===
        % Makes sure all triggers allow full epoch
        iMaxTheta(iMaxTheta <= abs(iEpochTime(1))) = [];
        iMaxTheta(iMaxTheta >= nTime - iEpochTime(end)) = [];
        % Loop on every peak of theta
        for i = 1:length(iMaxTheta)
            % Find epoch indices
            iTime = iMaxTheta(i) + iEpochTime;
            % Phase-triggered ERP of raw signal
            ERP(iSource,:) = ERP(iSource,:) + F(iSource, iTime) ./ length(iMaxTheta);
            % Phase-triggered time-frequency amplitude values (normalized)
            TF(iSource,:,:) = TF(iSource,:,:) + fs_high(1,iTime,:) ./ length(iMaxTheta);
        end
    end
end

author: brainstorm-tools


text:
function linplus_test05 ( )

%*****************************************************************************80
%
%% TEST05 tests R83_NP_DET, R83_NP_FA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  For a tridiagonal matrix that can be factored\n' );
  fprintf ( 1, '    with no pivoting,\n' );
  fprintf ( 1, '  R83_NP_FA factors,\n' );
  fprintf ( 1, '  R83_NP_DET computes the determinant.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r83_random ( n, seed );
%
%  Copy the matrix into general storage.
%
  b = r83_to_r8ge ( n, a );
%
%  Factor the matrix.
%
  [ a_lu, info ] = r83_np_fa ( n, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST05 - Warning!\n' );
    fprintf ( 1, '  R83_NP_FA returns INFO = %d\n', info );
  end

  r83_print ( n, a_lu, '  The factored R83 matrix:' );
%
%  Compute the determinant.
%
  det = r83_np_det ( n, a_lu );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R83_NP_DET computes determinant = %f\n', det );
%
%  Factor the matrix in R8GE storage.
%
  [ b_lu, info ] = r8ge_np_fa ( n, b );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST05 - Warning!\n' );
    fprintf ( 1, '  R8GE_NP_FA returns INFO = %d\n', info );
  end
%
%  Compute the determinant of the R8GE matrix.
%
  det = r8ge_np_det ( n, b_lu );

  fprintf ( 1, '  R8GE_NP_DET computes determinant = %f\n', det );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [points] = UniformPointsOnASphere(p,r)
% This function generates a point cloud equally distributed across the
% surface of a sphere.

% Input sanity check
if nargin < 2
    r = 1;
end

% Constants
viewDir = 3;
goldenRatio     = (1 + sqrt(5))/2;
angleIncrement  = pi*2*goldenRatio;

points = [];
for i = 1:numel(viewDir)
    t = i/viewDir;
    inclination = acos(1-2*t);
    azimuth = angleIncrement*i;

    % Convert spherical to cartesian
    x = p(1) + r*sin(inclination)*cos(azimuth);
    y = p(2) + r*sin(inclination)*sin(azimuth);
    z = p(3) + r*cos(inclination);

    % Assign point to matrix
    points(i,:) = [x,y,z];
end

end

%     static BoidHelper () {
%         directions = new Vector3[BoidHelper.numViewDirections];
%
%         float goldenRatio = (1 + Mathf.Sqrt (5)) / 2;
%         float angleIncrement = Mathf.PI * 2 * goldenRatio;
%
%         for (int i = 0; i < numViewDirections; i++) {
%             float t = (float) i / numViewDirections;
%             float inclination = Mathf.Acos (1 - 2 * t);
%             float azimuth = angleIncrement * i;
%
%             float x = Mathf.Sin (inclination) * Mathf.Cos (azimuth);
%             float y = Mathf.Sin (inclination) * Mathf.Sin (azimuth);
%             float z = Mathf.Cos (inclination);
%             directions[i] = new Vector3 (x, y, z);
%         }
%     }

repo: OpenMAS


text:
%kckhuemap 'creates a colormap for displaying Hue images with color'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ckhuemap.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% Double: saturation 'Saturation', default: 100: 'saturation for hsv->rgb conversion'
% Double: value 'Value ', default: 100: 'value for hsv->rgb conversion'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kckhuemap(i, {'i','';'saturation',100;'value',100;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ckhuemap - creates a colormap for displaying Hue images with color
%
%  DESCRIPTION
% Creates a colormap (map segment) for an hue image in order to display correctly the image with colors. User can enter value and saturation values.
% The user can also choose a hue model so the generated values will be displayed closer to the expected.
% For other hue models than HSV (so far only IHS) the colors will be approximated, (not correct and in most cases, plain wrong), because the other models cannot give correct hues for fixed saturation and values. If you're looking for a correct way to display only the hues of an image, HSV model will be better.
% If you want to use IHS, better results can be obtained by setting saturation < 60 and intensity < 90. Even so the results will not be correct.
%
%
%
%  EXAMPLES
% $COLOR/workspaces/HSV.wk
%
%  "SEE ALSO"
% iRGBtoHSV, iRGBtoHSI, iRGBtoHLS, iRGBtoSOWstar
%
%  RESTRICTIONS
% Value segment in output object is forced to KUINT for correct assignment of map values.
%
%  REFERENCES
% All references for the toolbox COLOR are on the $COLOR/Readme.html file.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

function varargout = kckhuemap(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kckhuemap(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'saturation', 100;'value', 100;'o', '__output'};
maxval={0,100,100,0};
minval={0,0,0,0};
istoggle=[0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','Double','Double','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'ckhuemap"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef MOMBIII < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Many objective metaheuristic based on the R2 indicator II
% alpha   ---   0.5 --- Threshold of variances
% epsilon --- 0.001 --- Tolerance threshold
% record  ---     5 --- The record size of nadir vectors

%------------------------------- Reference --------------------------------
% R. Hernandez Gomez and C. A. Coello Coello, Improved metaheuristic based
% on the R2 indicator for many-objective optimization, Proceedings of the
% Annual Conference on Genetic and Evolutionary Computation, 2015, 679-686.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [alpha,epsilon,recordSize] = Algorithm.ParameterSet(0.5,0.001,5);

            %% Generate random population
            [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
            Population    = Problem.Initialization();
            % Ideal and nadir points
            zmin = min(Population.objs,[],1);
            zmax = max(Population.objs,[],1);
            % For storing the nadir vectors of a few generations
            Record = repmat(zmax,recordSize,1);
            % For storing whether each objective has been marked for a few
            % generations
            Mark = false(recordSize,Problem.M);
            % R2 ranking procedure
            [Rank,Norm] = R2Ranking(Population.objs,W,zmin,zmax);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool  = TournamentSelection(2,Problem.N,Rank,Norm);
                Offspring   = OperatorGA(Problem,Population(MatingPool));
                Population  = [Population,Offspring];
                [Rank,Norm] = R2Ranking(Population.objs,W,zmin,zmax);
                [~,rank]    = sortrows([Rank,Norm]);
                Population  = Population(rank(1:Problem.N));
                Rank        = Rank(rank(1:Problem.N));
                Norm        = Norm(rank(1:Problem.N));
                [zmin,zmax,Record,Mark] = UpdateReferencePoints(Population.objs,zmin,zmax,Record,Mark,alpha,epsilon);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [plane, inlier_count, inliers] = ransac_depth_plane(x, depth, alpha,...
    iter, p)
%RANSAC_DEPTH_PLANE  Fit a 3D depth plane to a set of input points using RANSAC.
%The distance of a point from the plane is defined as the vertical offset.
%
%   INPUTS:
%
%   -|x|: |P|-by-2 matrix with image-plane coordinates of input points.
%
%   -|depth|: |P|-by-1 vector with depth values for input points.
%
%   -|alpha|: positive constant regulating the inlier threshold used in RANSAC.
%
%   -|iter|: maximum number of RANSAC iterations.
%
%   -|p|: target probability for having obtained a pure inlier set, used to cut
%   down on number of RANSAC iterations.
%
%   OUTPUTS:
%
%   -|plane|: 3-by-1 vector holding the parameters of the fitted plane, so that
%   for each point [x, y, d] it holds d = [x, y, 1] * plane.
%
%   -|inlier_count|: cardinality of the final inlier set.
%
%   -|inliers|: |inlier_count|-by-1 vector with indices of the members of the
%   final inlier set.

% Determine inlier threshold |theta| based on the depth values. Goal: invariance
% to range of depth values.
theta = alpha * median(depth);

% Number of plane parameters.
N = 3;

% Total number of points.
P = length(depth);

% Represent pixel coordinates homogeneously.
x_h = [x, ones(P, 1)];

% Initialize variables to save the best fit across iterations.
inliers = [];
inlier_count = 0;
inlier_ratio = 0;

% Main RANSAC loop.
for i = 1:iter
    colinear = true;
    while colinear
        % Select a sample at random.
        permut = randperm(P);
        sample = permut(1:N);
        x_h_sample = x_h(sample, :);

        % Proceed to plane fitting only if the points of the sample are not
        % colinear, i.e. square matrix with homogeneous coordinates of
        % points in the sample is of full rank.
        if rank(x_h_sample) == N
            colinear = false;
        end
    end

    % Fit a plane to the current sample.
    plane_tmp = fit_depth_plane(x_h_sample, depth(sample));

    % Identify inliers for current fit.

    % Vertical distance.
    inliers_tmp = abs(x_h * plane_tmp - depth) <= theta;
    % Orthogonal distance.
    % plane_tmp_normal = [plane_tmp(1:N-1); 1];
    % inliers_tmp = abs([points(known_pixels_finite, :), depth_known_finite] *...
    %     plane_tmp_normal - plane_tmp(N)) / norm(plane_tmp_normal) <= theta;
    inlier_tmp_count = nnz(inliers_tmp);

    % Update best fit.
    if inlier_tmp_count > inlier_count
        inlier_count = inlier_tmp_count;
        inlier_ratio = inlier_count / P;
        inliers = find(inliers_tmp);
    end

    % An adaptive version of RANSAC is implemented. If a pure inlier sample
    % has been already examined with probability greater that |p|, avoid
    % performing more iterations.
    p_bound = 1 - (1 - inlier_ratio ^ N) ^ i;
    if p <= p_bound
        break;
    end
end

% Perform a final least squares fit to compute the best depth plane with respect
% to all points in the maximum inlier set that has been found with RANSAC.
plane = fit_depth_plane(x_h(inliers, :), depth(inliers));

end

author: sakaridis


text:
function ret = spm_ov_goto_max(varargin)
% Goto maximum intensity tool - plugin for spm_orthviews
%
% This tool provides capabilities similar to the "Goto ... maximum"
% functionality in spm_mip_ui.m. When the tool is called for the first
% time, it has to read the whole image data file. This might result in a
% slow response depending on the image dimensions.
%
% This routine is a plugin to spm_orthviews. For general help about
% spm_orthviews and plugins type
%             help spm_orthviews
% at the MATLAB prompt.
%__________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Volkmar Glauche
% $Id: spm_ov_goto_max.m 4996 2012-10-11 18:28:37Z guillaume $

global st;
if isempty(st)
    error('goto_max: This routine can only be called as a plugin for spm_orthviews!');
end

if nargin < 2
    error('goto_max: Wrong number of arguments. Usage: spm_orthviews(''goto_max'', cmd, volhandle, varargin)');
end

cmd = lower(varargin{1});
volhandle = varargin{2};

switch cmd

    % Context menu and callbacks
    case 'context_menu'
        item0 = uimenu(varargin{3}, 'Label', 'Goto maximum');
        item1 = uimenu(item0, 'Label', 'Goto global maximum',...
            'Callback', ['spm_orthviews(''goto_max'', ''global'',' num2str(volhandle) ');']);
        item2 = uimenu(item0, 'Label', 'Goto nearest local maximum',...
            'Callback', ['spm_orthviews(''goto_max'', ''local'',' num2str(volhandle) ');']);
        ret = item0;

    case 'global'
        if ~isfield(st.vols{volhandle}, 'goto_max')
            [dat, xyz] = spm_read_vols(st.vols{volhandle});
            [unused, mxind] = max(dat(:));
            st.vols{volhandle}.goto_max.globalmm = xyz(:, mxind);
        end
        posmm = st.vols{volhandle}.premul*[st.vols{volhandle}.goto_max.globalmm;1];
        spm_orthviews('reposition', posmm(1:3));

    case 'local'
        % Read volume, goto local maxima higher or equal current voxel
        % intensity
        dat       = spm_read_vols(st.vols{volhandle});
        [x, y, z] = ndgrid(1:st.vols{volhandle}.dim(1), 1:st.vols{volhandle}.dim(2), 1:st.vols{volhandle}.dim(3));
        xyz       = [x(:) y(:) z(:)]';
        posvx     = round((st.vols{volhandle}.premul*st.vols{volhandle}.mat)\[spm_orthviews('pos');1]);
        try
            if isfinite(dat(posvx(1), posvx(2), posvx(3)))
                sel = isfinite(dat(:)) & dat(:) >= dat(posvx(1), posvx(2), posvx(3)) - eps;
            else
                sel = isfinite(dat(:));
            end
        catch
            sel = isfinite(dat(:));
        end
        [unused, unused, XYZ]  = spm_max(dat(sel), xyz(:,sel));
        XYZdist      = bsxfun(@minus,XYZ,posvx(1:3));
        [unused, nmaxind] = min(sum(XYZdist.^2));
        posmm        = st.vols{volhandle}.premul*st.vols{volhandle}.mat*[XYZ(:, nmaxind); 1];
        spm_orthviews('reposition', posmm(1:3));
end

author: spm


text:
%% Clean the slate
% The demo here in the main.html document can be run directly from main.m
% In this section the memory is cleared, some parameters are set, and the
% original images are loaded and displayed.
close all
clear
clc

% Setting up parameters and filenames
params.VERBOSE          = 1;                    %Get more plots out (True or False)
params.PLANE            = 3;                    %Which plane of the RGB is of interest to us for opperations on a single plane
params.BLACK_BACKGROUND = 1;                    %Images have a black background (True or False)
params.map                      = copper(256);  %Set colormap to something that is similar to wood
params.number_angular_divisions = 2^7;          %Powers of two are faster FFT
params.number_radial_divisions  = 40;           %Arbitrary constant
base_filename = 'dowel01.jpg';                  %Filename of non-moving image
move_filename = 'dowel02.jpg';                  %Filename of image that will move to the base image

base_image = imread(base_filename);    %reads in the image
move_image = imread(move_filename);    %reads in the image

% Data visualization.
subplot(1,2,1)
subimage(base_image);
title('Base image')
subplot(1,2,2)
subimage(move_image);
title('Image to move')
set (gcf, 'color', 'w')

%% Correct the X and Y displacements
% The two images are displaced from each other in the X and Y direcions.
% First the X and Y direction will be corrected by lining up the centroids
% of the two images.
%
% Method: First mask the image from background then find
% centroid of the image. Finally center the image based
% on its centroid.

base_plane_of_interest              = base_image(:,:,params.PLANE);
%Must use a single layer grayscale for most operations.  Blue plane is best
%contrast for these images
base_bw_plane_of_interest           = im2bw(base_plane_of_interest, graythresh(base_plane_of_interest));
%Turn to binary based on threshold gotten from graythresh
base_plane_of_interest_segmented    = bwmorph(base_bw_plane_of_interest, 'open');
%morphology
base_binary_mask                    = imfill (base_plane_of_interest_segmented, 'holes');
%fill it in

% Get some data about the new region
base_properties                     = regionprops(real(base_binary_mask),'all');
base_centroid_row                   = round(base_properties.Centroid(2));
base_centroid_col                   = round(base_properties.Centroid(1));

% Place a dot at the centroid, one for each layer
base_image(base_centroid_row, base_centroid_col, 1) = 255;
base_image(base_centroid_row, base_centroid_col, 2) = 255;
base_image(base_centroid_row, base_centroid_col, 3) = 255;

% Nice to have image be 0dd number of pixels in each diension
base_image       = make_odd_by_odd(base_image);
base_binary_mask = make_odd_by_odd(base_binary_mask);
% Grab the new size
[base_num_rows, base_num_cols, base_num_layers] = size(base_image);

% Where is the center of the image?
base_goal_row   = (base_num_rows - 1) / 2 + 1;
base_goal_col   = (base_num_cols - 1) / 2 + 1;

% how much do I need to move to center this?
base_delta_rows = base_goal_row - base_centroid_row;
base_delta_cols = base_goal_col - base_centroid_col;

%shift the images to be centered
base_image       = circshift(base_image      , [base_delta_rows, base_delta_cols]);
base_binary_mask = circshift(base_binary_mask, [base_delta_rows, base_delta_cols]);

% Same thing for the second image.  In production code this would be in
% a function, but this script was made for seminar presentation where all
% the code should be visible.
% Begin repeated code -------------------------------------
move_plane_of_interest              = move_image(:,:,params.PLANE);
move_bw_plane_of_interest           = im2bw(move_plane_of_interest, graythresh(move_plane_of_interest));
move_plane_of_interest_segmented    = bwmorph(move_bw_plane_of_interest, 'open');
move_binary_mask                    = imfill (move_plane_of_interest_segmented, 'holes');
move_properties                     = regionprops(real(move_binary_mask),'all');
move_centroid_row                   = round(move_properties.Centroid(2));
move_centroid_col                   = round(move_properties.Centroid(1));
move_image(move_centroid_row, move_centroid_col, 1) = 255;
move_image(move_centroid_row, move_centroid_col, 2) = 255;
move_image(move_centroid_row, move_centroid_col, 3) = 255;
move_image                          = make_odd_by_odd(move_image);
move_binary_mask                    = make_odd_by_odd(move_binary_mask);
[move_num_rows, move_num_cols, move_num_layers]     = size(move_image);
move_goal_row                       = (move_num_rows - 1) / 2 + 1;
move_goal_col                       = (move_num_cols - 1) / 2 + 1;
move_delta_rows                     = move_goal_row - move_centroid_row;
move_delta_cols                     = move_goal_col - move_centroid_col;
move_image                          = circshift(move_image      , [move_delta_rows, move_delta_cols]);
move_binary_mask                    = circshift(move_binary_mask, [move_delta_rows, move_delta_cols]);
% End of repeated code -------------------------------------

% Data visualization
subplot(1,2,1)
subimage(base_image);
title('Centered base image')
subplot(1,2,2)
subimage(move_image);
title('Centered moveable image')
set (gcf, 'color', 'w')

%%

% Data visualization
subplot(1,2,1)
subimage(base_binary_mask);
title('Centered base binary mask')
subplot(1,2,2)
subimage(move_binary_mask);
title('Centered movable binary mask')
set (gcf, 'color', 'w')

%% Convert the data from X-Y coordinates to Theta-Radius coordinates
% This section transforms the image such that each set of pixels that radiate
% out at a given angle from the centroid become a vertical line in a new image.
% These images can then be 'slid' left and right across each other to
% find the highest correlation between the images.

[base_radii, base_theta, base_intensity] = image2angularintesity(base_binary_mask, rgb2gray(base_image), params);
set (gcf, 'color', 'w')

%%
[move_radii, move_theta, move_intensity] = image2angularintesity(move_binary_mask, rgb2gray(move_image), params);
set (gcf, 'color', 'w')

%% Find the angular displacement by using frequency domain mathematics
% The Fast Fourier Tranform has some nice properties that allow the
% angle of highest correlation be found quickly and cleanly.
% The three dimensional plot of correlation as a function of
% angle and radius can be used to figure out which angular shift
% yields the highest correlation between the two images.
%
% Looking at the graph, it can be seen that most of the radius
% have reached a consenus as to which angular shift should be used, however
% some of the angles of smaller radius tend to disagree.  Because the
% smaller radii have less pixels, their votes are less informed, so must be
% given less weight in the final decision as to the proper angular shift.

BI = fft(base_intensity);  %Base Intensity FFT
MI = fft(move_intensity);  %Move Intensity FFT
CI = conj(BI) .* (MI);     %Correlation Intensity FFT
ci = real(ifft(CI));       %Correlation out of FFT space
nci = ci;                  %normalized correlation being built
nci = nci - repmat(mean(nci), size(nci,1), 1); % normalized correlation removes the mean value

%%
h = surf(move_radii(:,2:end), move_theta(:,2:end), nci(:,2:end));
set (h, 'linestyle', 'none')
title('Correlation of the two images')
xlabel('radius')
ylabel('angle')
zlabel('correlation')
colormap(jet)
shading interp
set (gcf, 'color', 'w')

%%
% Find the peak correlation
[value_of_correlation, indice_of_peak_correlation_at_given_radius] = max(ci);
% Do voting to find out how much to shift by
index_of_peak_correlation = best_correlation_index(indice_of_peak_correlation_at_given_radius, params);

%% Correct the theta displacement
% The angle decided upon in the previous section is now used as the input
% argument to an IMROTATE command which will rotate the image as needed.

% mathmatical housekeeping
degrees_per_step = 360 / params.number_angular_divisions;
angle_to_shift = index_of_peak_correlation * degrees_per_step;

% Finally rotate the images (mask and the image itself)
rotated_image       = imrotate(move_image      , angle_to_shift,'bicubic','crop');
rotated_binary_mask = imrotate(move_binary_mask, angle_to_shift,'bicubic','crop');

% Data visualization
subplot(1,3,1)
subimage(base_image)
title('Base Image')

subplot(1,3,2)
subimage(rotated_image)
title(['Image rotated by ' num2str(angle_to_shift) ' degrees'])

subplot(1,3,3)
subimage(move_image)
title('Image to rotate')
set (gcf, 'color', 'w')

%%
%    Doug Hull <hull@mathworks.com>     5/20/2002
%    Copyright 1984-2001 The MathWorks, Inc.
%    This function is not supported by The MathWorks, Inc.
%    It is provided 'as is' without any guarantee of
%    accuracy or functionality.
% Copyright 2002 - 2009 The MathWorks, Inc.


text:
function out = myresize(img, dim, type)
% resize to fixed dimension according to short side or long side
% --------------------------------------------------------------

if nargin < 3
    error('input not enough!');
end
if ~strcmp(type, 'short') && ~strcmp(type, 'long')
    error('wrong type!');
end
H = size(img, 1);
W = size(img, 2);

if strcmp(type, 'short')
    if H <= W
        out = imresize(img, [dim, NaN], 'bilinear');
    else
        out = imresize(img, [NaN, dim], 'bilinear');
    end
else
    if H <= W
        out = imresize(img, [NaN, dim], 'bilinear');
    else
        out = imresize(img, [dim, NaN], 'bilinear');
    end
end

author: cihangxie

repo: DAG

save_path: github-repos/MATLAB/cihangxie-DAG


text:
% TEST_STOKES_3D_SYMDRIVCAV_MP_H_DRCHLT: data function for Dirichlet boundary condition.

function h  = test_stokes_3d_symdrivcav_mp_h_drchlt (x, y, z, iside)
switch (iside)
 case {1, 2, 5, 6}
  h = zeros ([3, size(x)]);
 case 3
  h1 = -ones ([1, size(x)]);
  h2 = zeros ([1, size(x)]);
  h3 = zeros ([1, size(x)]);
  h  = cat (1, h1, h2, h3);
 case 4
  h1 = ones ([1, size(x)]);
  h2 = zeros ([1, size(x)]);
  h3 = zeros ([1, size(x)]);
  h  = cat (1, h1, h2, h3);
end
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes