text:
% EDICON Multi-edit and condense a training set
%
% 	J = EDICON(D,NSETS,NITERS,NTRIES)
%
% INPUT
%   D       Distance matrix dataset
%   NSETS   Number of subsets for editing, or [] for no editing (default: 3)
%   NITERS  Number of iterations for editing (default: 5)
%   NTRIES  Number of tries for condensing, or [] for no condensing (dflt: 10)
%
% OUTPUT
%   J       Indices of retained samples
%
% DESCRIPTION
% Returns the set of objects J such that the nearest neigbour gives zero error
% on the remaining objects. If MODE = 0, multi-edit the dataset represented by
% distance matrix D first. D can be computed from a dataset A by A*proxm(A).
%
% REFERENCES
% Devijver, P. and Kittler, J. "Pattern recognition", Prentice-Hall, 1982.
%
% SEE ALSO (PRTools Guide)
% PRDATASET, KNNC, PROXM

% Copyright: E. Pekalska and D. de Ridder, {E.Pekalska,D.deRidder}@ewi.tudelft.nl
% Faculty of Electrical Engineering, Mathematics and Computer Science
% Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands

function J = edicon (D,no_subsets,max_iter,no_tries)

    if (nargin < 4 | isempty(no_tries))
        prwarning(4,'number of tries for condensing not specified, assuming 10');
        no_tries = 10;
    end;
    if (nargin < 3 | isempty(max_iter))
        prwarning(4,'number of iterations for editing not specified, assuming 5');
        max_iter = 5;
    end;
    if (nargin < 2 | isempty(no_subsets))
        prwarning(4,'number of subsets for editing not specified, assuming 3');
        no_subsets = 5;
    end;

    % Extract dataset information.

	nlab    = getnlab(D); lablist = getlablist(D);
	[m,k,c] = getsize(D); p       = getprior(D);

    if (~isdataset(D) | (m ~= k))
        error('require a square distance matrix dataset');
    end

    J = 1:m;            % Initialise index array.

    % If requested, apply the multi-edit algorithm.

    if (~isempty(no_subsets))
        iter = 1;
        while (iter <= max_iter)

      	    % 1. Create NO_SUBSETS distinct subsets out of the sample indices
            %    remaining in J. Store the subset indices in subset_ind{.}.
            %    Note: this is slightly more complicated than strictly
            %    necessary, to enforce that each class is represented equally
            %    in each subset.

    		m = length(J); J = J(randperm(m));
		    subset_size = floor(m/no_subsets);

            % Find the MC(J) indices of samples belonging to class J in CLASS_IND{J}.
            % Also create a random permutation for class J.

            for j = 1:c
			    class_ind{j} = find(nlab(J)==j);
                mc(j)        = length(class_ind{j});
			    perm{j}      = randperm(mc(j));
    		end;

            % Distribute the samples over the subsets, per class.

    		for i = 1:no_subsets
		        subset_ind{i} = [];
			    for j = 1:c
                    num           = min(floor(mc(j)/no_subsets),length(class_ind{j}));
        		    subset_ind{i} = [subset_ind{i}; class_ind{j}(1:num)];
				    class_ind{j}  = class_ind{j}(num+1:end);
    			end;
			    perm_subset_ind{i} = J(subset_ind{i});
    		end;

          	% 2. In turn, classify each subset I using 1-NN on subset I+1.

          	drop = [];
  	        for i = 1:no_subsets
    			L1 = perm_subset_ind{i};
                L2 = perm_subset_ind{mod(i,no_subsets)+1};
    			if (length(L2)>1)
	  	    	    [dummy,nearest] = min(D(L2,L1));
    			else
				    nearest = 1;
    			end;
                drop = [drop; subset_ind{i}(find(nlab(L1)~=nlab(L2(nearest))))];
            end;

         	% 3. Discard incorrectly classified samples from J.

  	        if (isempty(drop))
  		        iter = iter + 1;
         	else
   		        J(drop) = []; iter = 0;
          	end;
        end;

        D = D(J,J);
    end;

    % Condense.

    if (~isempty(no_tries))
        % Extract dataset information.

    	nlab    = getnlab(D); lablist = getlablist(D);
	    [m,k,c] = getsize(D); p       = getprior(D);

        D = +D;

        % The whole procedure starts from a random object and continues to add
        % objects, so it is not optimal. Therefore it is repeated NO_TRIES times
        % and the smallest set found is returned.

        K = zeros(m,no_tries); storesz = zeros(1,no_tries);
        for o = 1:no_tries

            % Start with 1 sample index in STORE, the others in GRABBAG.
            p = randperm(m); store = p(1); grabbag = p(2:m);

            % While there are changes...
            transfer = 0;
            while (~transfer) & (~isempty(grabbag))
    	        storelab    = nlab(store);
	            grabbaglab  = nlab(grabbag);
	            new_grabbag = [];

                % For all samples in GRABBAG...
            	transfer = 0;
	            for k = 1:length(grabbag)
                    % ... find the nearest sample in STORE...
   	                [dummy,z] = min(D(grabbag(k),store));

                    % ... if it has a different label, move it to STORE.
   	                if (storelab(z) ~= grabbaglab(k))
                        store    = [store; grabbag(k)];
                        storelab = [storelab; grabbaglab(k)];
                        transfer = 1;
               	    else
                        new_grabbag = [new_grabbag; grabbag(k)];
                   	end;
    	        end;

            	grabbag = new_grabbag;

              end

              % Remember the STORE and its size for this repetition.
              storesz(o)        = length(store);
              K(1:storesz(o),o) = store;

        end

        % Take the STORE with minimal size.
        [minsz,minszind] = min(storesz);
        J = J(K(1:minsz,minszind));
    end;

return

author: marianux

repo: ecg-kit

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


text:
function a = r8cc_indicator ( m, n, nz_num, colptr, rowind )

%*****************************************************************************80
%
%% R8CC_INDICATOR sets up a R8CC indicator matrix.
%
%  Discussion:
%
%    The R8CC format is the double precision sparse compressed column
%    format.  Associated with this format, we have an M by N matrix
%    with NZ_NUM nonzero entries.  We construct the column pointer
%    vector COL of length N+1, such that entries of column J will be
%    stored in positions COL(J) through COL(J+1)-1.  This indexing
%    refers to both the ROW and A vectors, which store the row indices
%    and the values of the nonzero entries.  The entries of the
%    ROW vector corresponding to each column are assumed to be
%    ascending sorted.
%
%    The R8CC format is equivalent to the MATLAB "sparse" format,
%    and the Harwell Boeing "real unsymmetric assembled" (RUA) format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%
%    Input, integer N, the number of columns of the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero elements in A.
%
%    Input, integer COLPTR(N+1), points to the first element of each column.
%
%    Input, integer ROWIND(NZ_NUM), contains the row indices of the elements.
%
%    Output, real A(NZ_NUM), the matrix.
%
  fac = 10^( i4_log_10 ( n ) + 1 );

  for j = 1 : n
    for k = colptr(j) : colptr(j+1) - 1
      i = rowind(k);
      a(k) = fac * i + j;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p04_story ( )

%*****************************************************************************80
%
%% P04_STORY prints the "story" for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Theoretically, the data is a step, 0 to the left of 5, and 1\n' );
  fprintf ( 1, '  to the right.  To keep things simple, the data is defined\n' );
  fprintf ( 1, '  to be 0 up to 5 - RADIUS, 1/2 at 5, 1 at 5 + RADIUS and beyond,\n' );
  fprintf ( 1, '  with RADIUS set to a "small" value, currently 0.01.\n' );
  fprintf ( 1, '  Some interpolation methods will violently overreact to this jump.\n' );

  return
end


text:
f = [143 60];
A = [120 210; 110 30; 1 1];
b = [15000; 4000; 75];
lp = lp_maker(f, A, b, [-1; -1; -1], [], [], [], 1, 0);
solvestat = mxlpsolve('solve', lp)
format bank
obj = mxlpsolve('get_objective', lp)
format short
x = mxlpsolve('get_variables', lp)
mxlpsolve('delete_lp', lp);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function varargout = scl_slope(varargin)
% Format
% For getting the value
% dat = scl_slope(obj)
%
% For setting the value
% obj = scl_slope(obj,dat)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: scl_slope.m 1143 2008-02-07 19:33:33Z spm $

if nargin==2,
    varargout{1} = asgn(varargin{:});
elseif nargin==1,
    varargout{1} = ref(varargin{:});
else
    error('Wrong number of arguments.');
end;
return;

function dat = ref(obj)
dat = obj.scl_slope;
return;

function obj = asgn(obj,dat)
if isnumeric(dat), % && numel(dat)<=1,
    obj.scl_slope = double(dat);
else
    error('"scl_slope" must be numeric.');
end;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% demo for creating an SPM M/EEG dataset from arbitrary data using
% conversion of simple Fieldtrip raw data struct.
% SPM8 internal format is quite complex but is transparent to the user
% as meeg object's methods take care of maintaining its consistency.
% The most straightforward way to convert arbitrary data that is available
% as an ASCII file or *.mat file with some variables to SPM8 is to create
% a quite simple Fieldtrip raw data struct and then use SPM's
% spm_eeg_ft2spm to convert this struct to SPM8 file. Missing information
% can then be supplemented using meeg methods and SPM functions.
% Fieldtrip raw struct must contain the following fields:
% .fsample - sampling rate (Hz)
% .trial - cell array of matrices with identical dimensions channels x time
% .time - cell array of time vectors (in sec), the same length as the
%         second dimension of the data. For SPM8 they must be identical.
% .label - cell array of strings, list of channel labels. Same length as
%         the first dimension of the data.
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_convert_arbitrary_data.m 5404 2013-04-12 15:08:57Z vladimir $

% Initialize SPM
%--------------------------------------------------------------------------
spm('defaults','EEG');

% Some details about the data
%--------------------------------------------------------------------------
Nchannels = 5;
Nsamples  = 500;
Ntrials   = 50;
TimeOnset = -0.1; % in sec
Fsample = 1000;

chlabels = {
            'LFP1'
            'LFP2'
            'LFP3'
            'LFP4'
            'LFP5'
            };

% define the output file name
%--------------------------------------------------------------------------
fname = 'arbitrary_data_example';

% create data array
%--------------------------------------------------------------------------
data = randn([Nchannels, Nsamples, Ntrials]);

% create the time axis (should be the same for all trials)
%--------------------------------------------------------------------------
timeaxis = [0:(Nsamples-1)]./Fsample + TimeOnset;

% Create the Fieldtrip raw struct
%--------------------------------------------------------------------------

ftdata = [];

for i = 1:Ntrials
   ftdata.trial{i} = squeeze(data(:, :, i));
   ftdata.time{i} = timeaxis;
end

ftdata.fsample = Fsample;
ftdata.label = chlabels;
ftdata.label = ftdata.label(:);

% Convert the ftdata struct to SPM M\EEG dataset
%--------------------------------------------------------------------------
D = spm_eeg_ft2spm(ftdata, fname);

% Examples of providing additional information in a script
% [] comes instead of an index vector and means that the command
% applies to all channels/all trials.
%--------------------------------------------------------------------------
D = type(D, 'single');                        % Sets the dataset type
D = chantype(D, ':', 'LFP');                   % Sets the channel type
D = conditions(D, 1:Ntrials, 'Condition 1');  % Sets the condition label

% save
%--------------------------------------------------------------------------
save(D);

author: spm


text:
function test_bug3007

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3007'))

load cfg
load eye_movs

[cfgout, movement] = ft_detect_movement(cfg, eye_movs);

assert(isstruct(cfgout));
assert(~isstruct(movement)); % should be Nx3 array


text:
function jacobi_test01 ( )

%*****************************************************************************80
%
%% JACOBI_TEST01 tests JACOBI1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'JACOBI_TEST01:\n' );

  it_num = 400;
  n = 20;

  x_exact = ( 1 : n )';
  a = dif2 ( n );
  b = a * x_exact;

  x = zeros ( n, 1 );
  x_plot(1:n,1) = x;

  step = 1 : it_num + 1;
  r = nan ( it_num+1, 1 );
  xm = nan ( it_num+1, 1 );

  r(1,1) = ( norm ( a * x - b ) ).^2;

  for it = 1 : it_num

    x_new = jacobi1 ( n, a, b, x );
%   x_new = jacobi2 ( n, a, b, x );

    r(it+1,1) = ( norm ( a * x_new - b ) ).^2;
    x_plot(1:n,it+1) = x_new;
%
%  Display the residual.
%
    figure ( 1 )
    plot ( step, log ( r ), 'm-*' )
    title ( 'Log (Residual^2)' )
    xlabel ( 'Step' )
    ylabel ( 'Residual' )
    grid
%
%  Display the motion.
%
    xm(it,1) = sum ( ( x_new(:) - x(:) ).^2 ) / n;

    figure ( 2 )
    plot ( step, log ( xm ), 'm-*' )
    title ( 'Log (Average generator motion)' )
    xlabel ( 'Step' )
    ylabel ( 'Motion' )
    grid
%
%  Update the solution
%
    x = x_new;

  end
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 3 )

  y = ( 0 : it_num );
  for k = 1 : n
    plot ( x_plot(k,1:it_num+1), y )
    hold on;
  end
  grid on
  hold off;

  title ( 'Generator evolution.' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );

  return
end


text:
function planC = gamma_transPlanC(planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

oldScanSet = length(planC{indexS.scan});

for scanSet = 1:oldScanSet

    % the Scan Set will always be 1 as the 1st scan is deleted.
    planC = gamma_transScan(planC, 1);

end

for i = 1:length(planC{indexS.structures})
    planC{indexS.structures}(i).rasterized = 0;
    planC{indexS.structures}.rasterSegments = [];
end

planC{indexS.structureArray} = [];

planC = setUniformizedData(planC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% PLOTFRAME - plots a coordinate frame specified by a homogeneous transform
%
% Usage: function plotframe(T, len, label)
%
% Arguments:
%    T     - 4x4 homogeneous transform
%    len   - length of axis arms to plot (defaults to 1)
%    label - text string to append to x,y,z labels on axes
%
%  len and label are optional and default to 1 and '' respectively
%
% See also: ROTX, ROTY, ROTZ, TRANS, INVHT

% Copyright (c) 2001 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk at csse uwa edu au
% http://www.csse.uwa.edu.au/

function plotframe(T, len, label, colr)

    if ~all(size(T) == [4,4])
        error('plotframe: matrix is not 4x4')
    end

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

    if ~exist('label','var')
        label = '';
    end

    if ~exist('colr','var')
        colr = [0 0 1];
    end

    % Assume scale specified by T(4,4) == 1

    origin = T(1:3, 4);             % 1st three elements of 4th column
    X = origin + len*T(1:3, 1);     % point 'len' units out along x axis
    Y = origin + len*T(1:3, 2);     % point 'len' units out along y axis
    Z = origin + len*T(1:3, 3);     % point 'len' units out along z axis

    line([origin(1),X(1)], [origin(2), X(2)], [origin(3), X(3)], 'color', colr);
    line([origin(1),Y(1)], [origin(2), Y(2)], [origin(3), Y(3)], 'color', colr);
    line([origin(1),Z(1)], [origin(2), Z(2)], [origin(3), Z(3)], 'color', colr);

    text(X(1), X(2), X(3), ['x' label], 'color', colr);
    text(Y(1), Y(2), Y(3), ['y' label], 'color', colr);
    text(Z(1), Z(2), Z(3), ['z' label], 'color', colr);

author: christianwengert


text:
function msg = apm_lti(sys,name)

if nargin < 1,
    error('Input State Space or Transfer Function Model')
end
if nargin == 1,
    name = 'sys';
end
if nargin == 2,
    % trim whitespace
    filename = strtrim(name);
    % parse file name and extension
    [path,name,ext] = fileparts(name);
end
filename = [name '.apm'];

% Convert to state space form
sys = ss(sys);
% Check if discrete
if (sys.Ts>0),
    discrete = true;
else
    discrete = false;
end

% extract A, B, C, D matrices in sparse form
[n,m] = size(sys.B);
[p,m] = size(sys.D);

[ai,aj,av] = find(sparse(sys.A));
a = [ai,aj,av]';
[bi,bj,bv] = find(sparse(sys.B));
b = [bi,bj,bv]';
[ci,cj,cv] = find(sparse(sys.C));
if (size(sys.C,1)==1),
    c = [ci;cj;cv];
else
    c = [ci,cj,cv]';
end
[di,dj,dv] = find(sparse(sys.D));
d = [di,dj,dv]';
if(size(d,2)==0),
    d = [1,1,0]';
end

fid = fopen(filename,'w');
fprintf( fid,'\n');
fprintf( fid,'Objects \n');
fprintf( fid,['  ' name ' = lti\n']);
fprintf( fid,'End Objects \n');
fprintf( fid,'\n');
fprintf( fid,'Connections\n');
fprintf( fid,['  u[1:%d] = ' name '.u[1:%d]\n'],m,m);
fprintf( fid,['  x[1:%d] = ' name '.x[1:%d]\n'],n,n);
fprintf( fid,['  y[1:%d] = ' name '.y[1:%d]\n'],p,p);
fprintf( fid,'End Connections\n');
fprintf( fid,'\n');
fprintf( fid,'Model \n');
fprintf( fid,'  Parameters \n');
fprintf( fid,'    u[1:%d] = 0\n',m);
fprintf( fid,'  End Parameters \n');
fprintf( fid,'\n');
fprintf( fid,'  Variables \n');
fprintf( fid,'    x[1:%d] = 0\n',n);
fprintf( fid,'    y[1:%d] = 0\n',p);
fprintf( fid,'  End Variables \n');
fprintf( fid,'\n');
fprintf( fid,'  Equations \n');
fprintf( fid,'    ! add any additional equations here \n');
fprintf( fid,'  End Equations \n');
fprintf( fid,'End Model \n');
fprintf( fid,'\n');
fprintf( fid,'! dimensions\n');
fprintf( fid,'! (nx1) = (nxn)*(nx1) + (nxm)*(mx1)\n');
fprintf( fid,'! (px1) = (pxn)*(nx1) + (pxm)*(mx1)\n');
fprintf( fid,'!\n');
if discrete,
    fprintf( fid,'! discrete form with sampling time = %f\n',sys.Ts);
    fprintf( fid,'! x[k+1] = A * x[k] + B * u[k]\n');
    fprintf( fid,'!   y[k] = C * x[k] + D * u[k]\n');
else
    fprintf( fid,'! continuous form\n');
    fprintf( fid,'! dx/dt = A * x + B * u\n');
    fprintf( fid,'!   y   = C * x + D * u\n');
end
fprintf( fid,['File ' name '.txt\n']);
if discrete,
    fprintf( fid,'  sparse, discrete  ! dense/sparse, continuous/discrete\n');
else
    fprintf( fid,'  sparse, continuous  ! dense/sparse, continuous/discrete\n');
end
fprintf( fid,'  %d      ! m=number of inputs\n',m);
fprintf( fid,'  %d      ! n=number of states\n',n);
fprintf( fid,'  %d      ! p=number of outputs\n',p);
fprintf( fid,'End File\n');
fprintf( fid,'\n');
fprintf( fid,['File ' name '.a.txt \n']);
fprintf( fid,'   %d %d %f\n', a );
fprintf( fid,'End File \n');
fprintf( fid,'\n');
fprintf( fid,['File ' name '.b.txt \n']);
fprintf( fid,'   %d %d %f\n', b );
fprintf( fid,'End File \n');
fprintf( fid,'\n');
fprintf( fid,['File ' name '.c.txt \n']);
fprintf( fid,'   %d %d %f\n', c );
fprintf( fid,'End File \n');
fprintf( fid,'\n');
fprintf( fid,['File ' name '.d.txt \n']);
fprintf( fid,'   %d %d %f\n', d );
fprintf( fid,'End File \n');
fclose(fid);

msg = ['Created Linear Time Invariant (LTI) Model: ' filename];

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function varargout = MU_Seg_Edit(varargin)
% MU_SEG_EDIT MATLAB code for MU_Seg_Edit.fig
%      MU_SEG_EDIT, by itself, creates a new MU_SEG_EDIT or raises the existing
%      singleton*.
%
%      H = MU_SEG_EDIT returns the handle to a new MU_SEG_EDIT or the handle to
%      the existing singleton*.
%
%      MU_SEG_EDIT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MU_SEG_EDIT.M with the given input arguments.
%
%      MU_SEG_EDIT('Property','Value',...) creates a new MU_SEG_EDIT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before MU_Seg_Edit_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to MU_Seg_Edit_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help MU_Seg_Edit

% Last Modified by GUIDE v2.5 22-Sep-2013 01:30:53

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @MU_Seg_Edit_OpeningFcn, ...
                   'gui_OutputFcn',  @MU_Seg_Edit_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before MU_Seg_Edit is made visible.
function MU_Seg_Edit_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to MU_Seg_Edit (see VARARGIN)

handles.main_h=varargin{1};
set(handles.Seg_uitable,'Data',handles.main_h.V.Segs(:,1:3));

set(handles.Edit_pushbutton,'Enable','off');
set(handles.Update_pushbutton,'Enable','off');

% Choose default command line output for MU_Seg_Edit
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes MU_Seg_Edit wait for user response (see UIRESUME)
% uiwait(handles.MU_seg_edit_figure);

% --- Outputs from this function are returned to the command line.
function varargout = MU_Seg_Edit_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in Edit_pushbutton.
function Edit_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Edit_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.main_h=guidata(handles.main_h.MU_matrix_display);
if handles.main_h.V.Segs{handles.SegInd,2}~=handles.main_h.V.Slice
    return;
end

set(handles.Update_pushbutton,'Enable','on');
set(handles.Edit_pushbutton,'Enable','off');
MU_enable('off',[],handles.main_h);
handles.main_h=guidata(handles.main_h.MU_matrix_display);

eval(['handles.ROI_h=' handles.main_h.V.Segs{handles.SegInd,1} '(handles.main_h.Matrix_display_axes,handles.main_h.V.Segs{handles.SegInd,4});']);
Temp=handles.main_h.Mask(:,:,handles.main_h.V.Slice);
BW=createMask(handles.ROI_h);
Temp(BW~=0)=0;
handles.main_h.Mask(:,:,handles.main_h.V.Slice)=Temp;

guidata(hObject, handles);
guidata(handles.main_h.MU_matrix_display,handles.main_h);

% --- Executes on button press in Update_pushbutton.
function Update_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Update_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.main_h=guidata(handles.main_h.MU_matrix_display);
if handles.main_h.V.Segs{handles.SegInd,2}~=handles.main_h.V.Slice
    return;
end

set(handles.Update_pushbutton,'Enable','off');
MU_enable('on',[],handles.main_h);
handles.main_h=guidata(handles.main_h.MU_matrix_display);

Data=get(handles.Seg_uitable,'Data');

%deselect uitable ?? dirty code
set(handles.Seg_uitable, 'Data', []);
set(handles.Seg_uitable, 'Data', Data);

Temp=handles.main_h.Mask(:,:,handles.main_h.V.Slice);
BW=createMask(handles.ROI_h);
handles.main_h.V.Segs{handles.SegInd,3}=max(0,Data{handles.SegInd,3});
handles.main_h.V.Segs{handles.SegInd,4}=getPosition(handles.ROI_h);
Temp(BW~=0)=max(0,Data{handles.SegInd,3});

handles.main_h.Mask(:,:,handles.main_h.V.Slice)=Temp;
handles.main_h=MU_update_image(handles.main_h.Matrix_display_axes,{handles.main_h.TMatrix,handles.main_h.Mask},handles.main_h,0);
guidata(hObject, handles);
guidata(handles.main_h.MU_matrix_display,handles.main_h);

% --- Executes when entered data in editable cell(s) in Seg_uitable.
function Seg_uitable_CellEditCallback(hObject, eventdata, handles)
% hObject    handle to Seg_uitable (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) edited
%	PreviousData: previous data for the cell(s) edited
%	EditData: string(s) entered by the user
%	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
%	Error: error string when failed to convert EditData to appropriate value for Data
% handles    structure with handles and user data (see GUIDATA)

% --- Executes when selected cell(s) is changed in Seg_uitable.
function Seg_uitable_CellSelectionCallback(hObject, eventdata, handles)
% hObject    handle to Seg_uitable (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
% Indices: row and column indices of the cell(s) currently selecteds
% handles    structure with handles and user data (see GUIDATA)
if isempty(eventdata.Indices)
    return;
end

handles.main_h=guidata(handles.main_h.MU_matrix_display);
handles.SegInd=eventdata.Indices(1);
if handles.main_h.V.Segs{handles.SegInd,2}~=handles.main_h.V.Slice
    guidata(hObject, handles);
    return;
end

if eventdata.Indices(2)==1
    Pos=handles.main_h.V.Segs{handles.SegInd,4};
    axes(handles.main_h.Matrix_display_axes);
    text(Pos(1,1),Pos(1,2),num2str(handles.SegInd),'Color','r','FontSize',18);
end

if strcmp(get(handles.Update_pushbutton,'Enable'),'off')
    set(handles.Edit_pushbutton,'Enable','on');
end

guidata(hObject, handles);

% --- Executes on mouse motion over figure - except title and menu.
function MU_seg_edit_figure_WindowButtonMotionFcn(hObject, eventdata, handles)
% hObject    handle to MU_seg_edit_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if ~ishandle(handles.main_h.MU_matrix_display);
    close(hObject);
    return;
end

handles.main_h=guidata(handles.main_h.MU_matrix_display);
Data=get(handles.Seg_uitable,'Data');
if length(handles.main_h.V.Segs(:,3)) > length(Data(:,3))
    Data(end+1,3)={1};
end
set(handles.Seg_uitable,'Data',[handles.main_h.V.Segs(:,1:2) Data(:,3)]);
guidata(hObject, handles);

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function A = spm_clusters(L,n)
% Return the cluster index for a point list
% FORMAT A = spm_clusters(L,n)
% L     - locations [x y x]' {in voxels} ([3 x m] matrix)
% n     - connectivity criterion (see spm_bwlabel) [Default: 18]
%
% A     - cluster index or region number ([1 x m] vector)
%__________________________________________________________________________
%
% spm_clusters characterises a point list of voxel values defined with
% their locations (L) in terms of edge, face and vertex connected
% subsets, returning a list of indices in A, such that the ith location
% belongs to cluster A(i) (using an 18 connectivity scheme).
%__________________________________________________________________________
% Copyright (C) 1994-2012 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: spm_clusters.m 4929 2012-09-17 14:21:01Z guillaume $

if isempty(L), A = []; return; end
if nargin < 2, n = 18; end

% Turn location list to binary 3D volume
%--------------------------------------------------------------------------
dim       = [max(L(1,:)) max(L(2,:)) max(L(3,:))];
vol       = zeros(dim(1),dim(2),dim(3));
indx      = sub2ind(dim,L(1,:)',L(2,:)',L(3,:)');
vol(indx) = 1;

% Label each cluster in 3D volume with its own label using an 18
% connectivity criterion
%--------------------------------------------------------------------------
cci = spm_bwlabel(vol,n);

% Map back to list
%--------------------------------------------------------------------------
A = cci(indx);
A = A(:)';


text:
function Cs = subgraph_centrality(CIJ)
% SUBGRAPH_CENTRALITY     subgraph centrality of a network
%
%   Cs = subgraph_centrality(CIJ)
%
%   The subgraph centrality of a node is a weighted sum of closed walks of
%   different lengths in the network starting and ending at the node. This
%   function returns a vector of subgraph centralities for each node of the
%   network.
%
%   Inputs:     CIJ,        adjacency matrix (binary)
%
%   Outputs:     Cs,        subgraph centrality
%
%   Reference: Estrada and Rodriguez-Velasquez (2005) Phys Rev E 71, 056103
%              Estrada and Higham (2010) SIAM Rev 52, 696.
%
%   Xi-Nian Zuo, Chinese Academy of Sciences, 2010
%   Rick Betzel, Indiana University, 2012

[V,lambda] = eig(CIJ);                 % Compute the eigenvectors and
lambda     = diag(lambda);             % eigenvalues.
V2         = V.^2;                     % Matrix of squares of the eigenvectors elements.
Cs         = real(V2 * exp(lambda));   % Compute eigenvector centrality. Lop off imaginary part remaining due to precision error.


text:
close all;
num = 100;
bias = 0.3;
% margin = 0.5;
class1 = randn(num,3)/4;
class1_bias = bsxfun(@plus, class1, [-bias, 0, bias]);
class2 = randn(num,3)/4;
class2_bias = bsxfun(@plus, class2, [bias, 0, bias]);

net = caffe.Net('normalize_layer.prototxt','train');
f = net.forward({reshape([class1_bias; class2_bias]',[1,1,3,num*2]),reshape([zeros(1,num) ones(1,num)],[1,1,1,num*2])});

all_norm = squeeze(net.blobs('norm1').get_data())';
c1_norm = all_norm(1:num,:);
c2_norm = all_norm(num+1:end,:);

c1_center = mean(c1_norm);
c1_center = c1_center ./ norm(c1_center);
c2_center = mean(c2_norm);
c2_center = c2_center ./ norm(c2_center);

figure(1);
hold on;
scatter3(c1_norm(:,1),c1_norm(:,2),c1_norm(:,3),200,'r.');
scatter3(c2_norm(:,1),c2_norm(:,2),c2_norm(:,3),200,'b.');
% plot3(c1_center(1),c1_center(2),c1_center(3),'r.','MarkerSize',40);
% plot3(c2_center(1),c2_center(2),c2_center(3),'b.','MarkerSize',40);
plot3(c1_center(1),c1_center(2),c1_center(3),'ro','MarkerSize',10, 'LineWidth',5);
plot3(c2_center(1),c2_center(2),c2_center(3),'bo','MarkerSize',10, 'LineWidth',5);

w = [c1_center;c2_center];
net.layers('id_weight').params(1).set_data(w');

for i=1:1000
    f = net.forward({reshape([class1_bias; class2_bias]',[1,1,3,num*2]),reshape([zeros(1,num) ones(1,num)],[1,1,1,num*2])});
    w2 = net.layers('id_weight').params(1).get_data()';
    assert(sum(abs(w(:)-w2(:))) < 0.001);
    w = w2;
    g = net.backward({[1], [1;1]});
    gw1 = net.layers('id_weight').params(1).get_diff()';

    w = w - 0.00001*gw1;
    net.layers('id_weight').params(1).set_data(w');
end;
w = squeeze(net.blobs('id_weight_normalize').get_data())';
w_norm = squeeze(net.blobs('id_weight_normalize').get_data())';
plot3(w_norm(1,1),w_norm(1,2),w_norm(1,3),'rx','MarkerSize',20, 'LineWidth',5);
plot3(w_norm(2,1),w_norm(2,2),w_norm(2,3),'bx','MarkerSize',20, 'LineWidth',5);

direct_point = reshape(g{1},[3, num*2])';
update_point = [c1_norm; c2_norm] - 10 * direct_point;
update_point = bsxfun(@rdivide, update_point, sqrt(sum(update_point.^2,2)));
arrow3([c1_norm; c2_norm],update_point, 'k-0.5', 0.5,1,0.1);
xlabel('dimension 1');
ylabel('dimension 2');
zlabel('dimension 3');
hold off;
box on;

author: happynear

repo: NormFace

save_path: github-repos/MATLAB/happynear-NormFace


text:
function [filterfunc,winbw] = helper_filtergeneratorfunc(wintype,winCell,fs,bwmul,min_win,trunc_at,audscale,do_subprec,do_symmetric,do_warped)
firwinflags=getfield(arg_firwin,'flags','wintype');
freqwinflags=getfield(arg_freqwin,'flags','wintype');
probelen = 10000;

subprecflag = 'pedantic';
if ~do_subprec, subprecflag = 'nopedantic'; end

switch wintype
    case firwinflags
        winbw=norm(firwin(wintype,probelen)).^2/probelen;
        % This is the ERB-type bandwidth of the prototype

        if do_symmetric
            filterfunc = @(fsupp,fc,scal)...
                         blfilter(winCell,fsupp,fc,'fs',fs,'scal',scal,...
                                  'inf','min_win',min_win,subprecflag);
        else
            fsupp_scale=1/winbw*bwmul;
            filterfunc = @(fsupp,fc,scal)...
                         warpedblfilter(winCell,fsupp_scale,fc,fs,...
                                        @(freq) freqtoaud(freq,audscale),...
                                        @(aud)  audtofreq(aud,audscale),...
                                        'scal',scal,'inf');
        end
        bwtruncmul = 1;
    case freqwinflags
        if do_warped
            error('%s: TODO: Warping is not supported for windows from freqwin.',...
                upper(mfilename));
        end

        probebw = 0.01;

        % Determine where to truncate the window
        H = freqwin(winCell,probelen,probebw);
        winbw = norm(H).^2/(probebw*probelen/2);
        bwrelheight = 10^(-3/10);

        if trunc_at <= eps
            bwtruncmul = inf;
        else
            try
                bwtruncmul = winwidthatheight(abs(H),trunc_at)/winwidthatheight(abs(H),bwrelheight);
            catch
                bwtruncmul = inf;
            end
        end

        filterfunc = @(fsupp,fc,scal)...
                     freqfilter(winCell, fsupp, fc,'fs',fs,'scal',scal,...
                                'inf','min_win',min_win,...
                                'bwtruncmul',bwtruncmul,subprecflag);
end

function width = winwidthatheight(gnum,atheight)

width = zeros(size(atheight));
for ii=1:numel(atheight)
    gl = numel(gnum);
    gmax = max(gnum);
    frac=  1/atheight(ii);
    fracofmax = gmax/frac;

    ind =find(gnum(1:floor(gl/2)+1)==fracofmax,1,'first');
    if isempty(ind)
        %There is no sample exactly half of the height
        ind1 = find(gnum(1:floor(gl/2)+1)>fracofmax,1,'last');
        ind2 = find(gnum(1:floor(gl/2)+1)<fracofmax,1,'first');
        rest = 1-(fracofmax-gnum(ind2))/(gnum(ind1)-gnum(ind2));
        width(ii) = 2*(ind1+rest-1);
    else
        width(ii) = 2*(ind-1);
    end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function sVF = minus(sVF1, sVF2)
%
% Syntax
%   sVF = sVF1-sVF2
%   sVF = a-sVF1
%   sVF = sVF1-a
%

sVF = plus(sVF1, -sVF2);

end

author: mtex-toolbox

repo: mtex

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


text:
T=1/f;

sl=3*T; % signal length

xl=get(handles.axes1,'Xlim');
yl=get(handles.axes1,'Ylim');
dxl=xl(2)-xl(1);
dyl=yl(2)-yl(1);

xc1=0:dxl/res:sl;
xc=mod(xc1,T); % turn to one period
if iscnt
    yc=interp1([xys(1,:)-dxl xys(1,:) xys(1,:)+dxl],[xys(2,:) xys(2,:) xys(2,:)],xc,mth,'extrap');
else
    yc=interp1(xys(1,:),xys(2,:),xc,mth,'extrap');
end

nf=max(abs(yc)); % normalization factor
yc=yc/nf; % mormalize

figure;
plot(xc1,yc,'b-',xys(1,:),xys(2,:)/nf,'rx',[T T],[-1 1],'k--',[2*T 2*T],[-1 1],'k--');
xlabel('time, s');
ylabel('signal, notmalized');
title('3 periods of the signal');
ylim([-1 1]);

author: Sable

repo: mcbench-benchmarks


text:
function option_num = p03_option_num ( )

%*****************************************************************************80
%
%% P03_OPTION_NUM returns the number of options for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer OPTION_NUM, the number of options.
%
  option_num = 4;

  return
end


text:
function [view,area3d] = measureFlatROIAreaMesh(view, grayThickness)
%
% [view,area3d] = measureFlatROIAreaMesh(view, [grayThickness])
%
% Measures the area on the 3D cortical surface of the current
% flat-map ROI.
%
%    1. Select ROI in the Flat window.
%
%    2. Use the Area Measure pull down in the Analysis pulldown of
% the Flat window.
%
% NOTE: this code must find the flat.mat file associated with the
% unfold, load it, and find the 'unfoldMeshSummary' field. This extra
% information is saved from within newer versions of mrFlatMesh when
% you check 'save extra info'.
%
% The method here is much preferred to that in measureFlatROIArea, which
% tries to create it's own triangulation based only on layer 1 nodes. This
% usually works OK for ROIs that are convex, but can severely over-estimate
% the area when the ROI has concavities. Here we just use the triangles
% from the gray/white interface surface mesh.
%
% See Also: measureFlatROIArea, sumTriangularArea
%
% HISTORY:
% 2003.01.07 RFD (bob@white.stanford.edu) wrote it.
% 2003.01.30 RFD: minor cleaning, tested some alternative vertex inclusion
% criteria.
% 2006.04.10 RFD: addes comments about our preference for this function
% over the older measureFlatROIArea.

global mrSESSION;

if(~exist('grayThickness','var') | isempty(grayThickness))
    grayThickness = 3;
end

% If someone want the area returned, then we assume that we are scripted
% and thus don't use gui stuff.
if(nargout<2),  gui = 1;
else,  gui = 0; end

% Get selpts from current ROI
if view.selectedROI,   ROIcoords = getCurROIcoords(view);
else  error('No current ROI'); end
ROIname = view.ROIs(view.selectedROI).name;

if isempty(ROIcoords), error('ROI is empty!'); end

hemi = ROIcoords(3,1);
if(~all(ROIcoords(3,:) == hemi))
    error('ROI spans both hemispheres!');
end

% load the flat.mat file
if(hemi==1)
    hemiName = 'left';
    unfoldFile = view.leftPath;
else
    hemiName  = 'right';
    unfoldFile = view.rightPath;
end
if(~exist(unfoldFile,'file'))
    % Try to fix broken paths
    warning([unfoldFile ' not found- trying to find it...']);
    anatPath = getAnatomyPath(mrSESSION.subject);
    indx = findstr(mrSESSION.subject,unfoldFile);
    if(~isempty(indx))
        unfoldFile = fullfile(fileparts(anatPath), unfoldFile(indx(1):end));
        unfoldFile = strrep(unfoldFile, '\', filesep);
        unfoldFile = strrep(unfoldFile, '/', filesep);
    end
    if(exist(unfoldFile,'file'))
        warning(['Using a similarly-named file (' unfoldFile ')...']);
    else
        disp(unfoldFile);
        [f,p] = myUiGetFile(anatPath, {'*.mat';'*.*'}, ['Select ' hemiName ' unfold file...']);
        unfoldFile = fullfile(p,f);
    end
end

if(~isfield(view,'mesh') | size(view.mesh) < hemi | isempty(view.mesh{hemi}))
    unfold = load(unfoldFile);
    disp(['loaded unfold from ', unfoldFile,'...']);
    if(~isfield(unfold,'unfoldMeshSummary'))
        disp('*************************************************************');
        disp('This flat.mat file does not have an unfoldMeshSummary field!');
        disp('Redo the unfold and tell mrFlatMesh to save the extra info.');
        disp('For now, you can get an area estimate from measureFlatROIArea.');
        error('Missing the unfoldMeshSummary field.');
    end
    view.mesh{hemi} = unfold.unfoldMeshSummary;
end

locs2d = view.mesh{hemi}.locs2d;

% To convert raw glocs2d values to those used in mrLoadRet, we do:
minLocs2d = min(locs2d);
locs2d(:,1) = locs2d(:,1) - minLocs2d(1) + 1;
locs2d(:,2) = locs2d(:,2) - minLocs2d(2) + 1;
locs2d = locs2d';

scaledVertices(:,1) = view.mesh{hemi}.uniqueVertices(:,1) .* view.mesh{hemi}.scaleFactor(1);
scaledVertices(:,2) = view.mesh{hemi}.uniqueVertices(:,2) .* view.mesh{hemi}.scaleFactor(2);
scaledVertices(:,3) = view.mesh{hemi}.uniqueVertices(:,3) .* view.mesh{hemi}.scaleFactor(3);
areaList3d = findFaceArea(view.mesh{hemi}.connectionMatrix, ...
    scaledVertices, ...
    view.mesh{hemi}.uniqueFaceIndexList);
% Lets also measure the 2d area, just for fun.
areaList2d = findFaceArea(view.mesh{hemi}.connectionMatrix, ...
    [locs2d',zeros(size(locs2d,2),1)], ...
    view.mesh{hemi}.uniqueFaceIndexList);

%
% Find the triangles that overlap with the ROI.
%

% First we find all the vertices that fall within the ROI.
%
% Intersect returns only the unique matches, but we want them all.
% So, we use ismember, which tells us for each row of the first matrix if
% it matches any row in the second matrix.
roiVertexIndices = find(ismember(round(locs2d)', ROIcoords(1:2,:)', 'rows'));

% Now we find the triangles (faces) within the ROI by a simple heuristic-
% any one of the triangle's vertices is in the ROI. What we should really
% do is select only those triangles whose AREA falls mostly within the ROI.
% But that's too hard (ie. would be slow in matlab).
roiFaceIndices = find(ismember(view.mesh{hemi}.uniqueFaceIndexList(:,1), roiVertexIndices) ...
    | ismember(view.mesh{hemi}.uniqueFaceIndexList(:,2), roiVertexIndices) ...
    | ismember(view.mesh{hemi}.uniqueFaceIndexList(:,3), roiVertexIndices));
% 2003.01.30 RFD: NOTE: I tried making the above inclusion criteria more
% conservative- so that two of it's vertices needed to be within the ROI.
% However, this wasn't much better, especially for thin ROIs, since it
% erred in the other direction and eliminated too many triangles.
% But that more conservative criteria would have the very desireable
% property that any triangle will only get included in one of
% two abutting ROIs, never in both like the old method sometimes did.
% vertOne = ismember(view.mesh{hemi}.uniqueFaceIndexList(:,1), RoiVertexIndices);
% vertTwo = ismember(view.mesh{hemi}.uniqueFaceIndexList(:,2), RoiVertexIndices);
% vertThree = ismember(view.mesh{hemi}.uniqueFaceIndexList(:,3), RoiVertexIndices);
% roiFaceIndices = find((vertOne & vertTwo) | (vertOne & vertThree) | (vertTwo & vertThree));

area3d = sum(areaList3d(roiFaceIndices));
area2d = sum(areaList2d(roiFaceIndices));

areaStr = addText([],sprintf('%s area = %.0f mm^2 (%d triangles)\n',ROIname,area3d,length(roiFaceIndices)));
areaStr = addText(areaStr,sprintf('2d = %.0f~mm^2 and %.0f voxels',area2d,size(ROIcoords,2)));

if(gui)
    %msgbox([areaStr],[ROIname ' area']);
    % Make the ROI coords point to the upper left of each ROI pixel
    ROIcoords = ROIcoords-0.5;
    % now define four the lines for each ROI coordinate.
    roiLinesX = [[ ROIcoords(2,:);   ROIcoords(2,:)+1 ],...
            [ ROIcoords(2,:);   ROIcoords(2,:)   ],...
            [ ROIcoords(2,:)+1; ROIcoords(2,:)+1 ],...
            [ ROIcoords(2,:)+1; ROIcoords(2,:)   ]];
    roiLinesY =-[[ ROIcoords(1,:);   ROIcoords(1,:)   ],...
            [ ROIcoords(1,:);   ROIcoords(1,:)+1 ],...
            [ ROIcoords(1,:)+1; ROIcoords(1,:)   ],...
            [ ROIcoords(1,:)+1; ROIcoords(1,:)+1 ]];

    roiFaces = view.mesh{hemi}.uniqueFaceIndexList(roiFaceIndices,:);
    % This will color each triangle with the z-depth of the surface
    % (something like curvature).
    %colors = scaledVertices(roiFaces(:,1),3) + scaledVertices(roiFaces(:,2),3) + scaledVertices(roiFaces(:,3),3);
    %colors = colors-mean(colors);
    % Use abs to cope with the fact that we sometimes get negative numbers here...
    colors = abs(log(areaList3d(roiFaceIndices)./areaList2d(roiFaceIndices)));
    %areaStr = [areaStr ' (Color shows log 3D/2D area)'];

    figure;  hold on; axis equal; colormap(hot);
    patch([locs2d(2,roiFaces(:,1));...
            locs2d(2,roiFaces(:,2));...
            locs2d(2,roiFaces(:,3))], ...
        -[locs2d(1,roiFaces(:,1));...
            locs2d(1,roiFaces(:,2));...
            locs2d(1,roiFaces(:,3))], colors');
    lineH = line(roiLinesX, roiLinesY,'Color','blue');
    colorbar;
    %plot(locs2d(2,roiVertexIndices), -locs2d(1,roiVertexIndices), '.r');
    %plot(round(locs2d(2,roiVertexIndices)), -round(locs2dInt(1,roiVertexIndices)), '.g');
    hold off;
    title(areaStr);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% CONVERTLOCS - Convert electrode locations between coordinate systems
%                 using the EEG.chanlocs structure.
%
% Usage: >> newchans = convertlocs( chanlocs, 'command');
%
% Input:
%   chanlocs  - An EEGLAB EEG dataset OR a EEG.chanlocs channel locations structure
%   'command' - ['cart2topo'|'sph2topo'|'sphbesa2topo'| 'sph2cart'|'topo2cart'|'sphbesa2cart'|
%               'cart2sph'|'sphbesa2sph'|'topo2sph'| 'cart2sphbesa'|'sph2sphbesa'|'topo2sphbesa'|
%               'cart2all'|'sph2all'|'sphbesa2all'|'topo2all']
%                These command modes convert between four coordinate frames: 3-D Cartesian
%                (cart), Matlab spherical (sph), Besa spherical (sphbesa), and 2-D polar (topo)
%               'auto' -- Here, the function finds the most complex coordinate frame
%                 and constrains all the others to this one. It searches first for Cartesian
%                 coordinates, then for spherical and finally for polar. Default is 'auto'.
%
% Optional input
%   'verbose' - ['on'|'off'] default is 'off'.
%
% Outputs:
%   newchans - new EEGLAB channel locations structure
%
% Ex:  CHANSTRUCT = convertlocs( CHANSTRUCT, 'cart2topo');
%      % Convert Cartesian coordinates to 2-D polar (topographic).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 22 Dec 2002
%
% See also: READLOCS

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 22 Dec 2002, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function chans = convertlocs(chans, command, varargin)

if nargin < 1
   help convertlocs;
   return;
end

if ~isfield(chans, 'theta') && ~isfield(chans, 'X') && ~isfield(chans, 'radius') && ~isfield(chans, 'sph_theta_besa')
    return
end

if nargin < 2
   command = 'auto';
end
if nargin == 4 && strcmpi(varargin{2}, 'on')
    verbose = 1;
else
    verbose = 0; % off
end

% test if value exists for default
% --------------------------------
if strcmp(command, 'auto')
    if isfield(chans, 'X') && any(~cellfun(@isempty, { chans.X }))
        command = 'cart2all';
        if verbose
            disp('Make all coordinate frames uniform using Cartesian coords');
        end
    else
        if isfield(chans, 'sph_theta') && ~isempty(chans(1).sph_theta)
            command = 'sph2all';
            if verbose
                disp('Make all coordinate frames uniform using spherical coords');
            end
        else
            if isfield(chans, 'sph_theta_besa') && ~isempty(chans(1).sph_theta_besa)
                command = 'sphbesa2all';
                if verbose
                    disp('Make all coordinate frames uniform using BESA spherical coords');
                end
            else
                command = 'topo2all';
                if verbose
                    disp('Make all coordinate frames uniform using polar coords');
                end
            end
        end
    end
end

% convert
% -------
switch command
 case 'topo2sph'
   theta  = {chans.theta};
   radius = {chans.radius};
   indices = find(~cellfun('isempty', theta));
   [sph_phi, sph_theta] = topo2sph( [ [ theta{indices} ]' [ radius{indices}]' ] );
   if verbose
       disp('Warning: electrodes forced to lie on a sphere for polar to 3-D conversion');
   end
   for index = 1:length(indices)
      chans(indices(index)).sph_theta  = sph_theta(index);
      chans(indices(index)).sph_phi    = sph_phi  (index);
   end
   if isfield(chans, 'sph_radius')
       meanrad = mean([ chans(indices).sph_radius ]);
       if isempty(meanrad)
           [chans(indices).sph_radius] = deal(85);
           meanrad = 85;
       end
   else
        [chans(indices).sph_radius] = deal(85);
        meanrad = 85;
   end
   sph_radius(1:length(indices)) = {meanrad};
case 'topo2sphbesa'
   chans = convertlocs(chans, 'topo2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2sphbesa', varargin{:}); % search for spherical coords
case 'topo2cart'
   chans = convertlocs(chans, 'topo2sph', varargin{:}); % search for spherical coords
   if verbose
       disp('Warning: spherical coordinates automatically updated');
   end
   chans = convertlocs(chans, 'sph2cart', varargin{:}); % search for spherical coords
case 'topo2all'
   chans = convertlocs(chans, 'topo2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2sphbesa', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2cart', varargin{:}); % search for spherical coords
case 'sph2cart'
   sph_theta  = {chans.sph_theta};
   sph_phi    = {chans.sph_phi};
   indices = find(~cellfun('isempty', sph_theta));
   if ~isfield(chans, 'sph_radius')
        [chans(indices).sph_radius] = deal(85);
       sph_radius(1:length(indices)) = {85};
   else
       sph_radius = {chans.sph_radius};
   end
   inde = find(cellfun('isempty', sph_radius));
   if ~isempty(inde)
       meanrad = mean( [ sph_radius{:} ]);
       sph_radius(inde) = { meanrad };
   end
   [x, y, z] = sph2cart([ sph_theta{indices} ]'/180*pi, [ sph_phi{indices} ]'/180*pi, [ sph_radius{indices} ]');
   for index = 1:length(indices)
      chans(indices(index)).X = x(index);
      chans(indices(index)).Y = y(index);
      chans(indices(index)).Z = z(index);
   end
case 'sph2topo'
 if verbose
     % disp('Warning: all radii constrained to one for spherical to topo transformation');
 end
 sph_theta  = {chans.sph_theta};
 sph_phi    = {chans.sph_phi};
 indices = find(~cellfun('isempty', sph_theta));
 [chan_num,angle,radius] = sph2topo([ ones(length(indices),1)  [ sph_phi{indices} ]' [ sph_theta{indices} ]' ], 1, 2); % using method 2
 for index = 1:length(indices)
     chans(indices(index)).theta  = angle(index);
     chans(indices(index)).radius = radius(index);
     if ~isfield(chans, 'sph_radius') || isempty(chans(indices(index)).sph_radius)
         chans(indices(index)).sph_radius = 85;
     end
 end
case 'sph2sphbesa'
   % using polar coordinates
   sph_theta  = {chans.sph_theta};
   sph_phi    = {chans.sph_phi};
   indices = find(~cellfun('isempty', sph_theta));
   [chan_num,angle,radius] = sph2topo([ones(length(indices),1)  [ sph_phi{indices} ]' [ sph_theta{indices} ]' ], 1, 2);
   [sph_theta_besa, sph_phi_besa] = topo2sph([angle radius], 1, 1);
   for index = 1:length(indices)
      chans(indices(index)).sph_theta_besa  = sph_theta_besa(index);
      chans(indices(index)).sph_phi_besa    = sph_phi_besa(index);
   end
case 'sph2all'
   chans = convertlocs(chans, 'sph2topo', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2sphbesa', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2cart', varargin{:}); % search for spherical coords
case 'sphbesa2sph'
   % using polar coordinates
   sph_theta_besa  = {chans.sph_theta_besa};
   sph_phi_besa    = {chans.sph_phi_besa};
   indices = find(~cellfun('isempty', sph_theta_besa));
   [chan_num,angle,radius] = sph2topo([ones(length(indices),1)  [ sph_theta_besa{indices} ]' [ sph_phi_besa{indices} ]' ], 1, 1);
   %for index = 1:length(chans)
   %   chans(indices(index)).theta  = angle(index);
   %   chans(indices(index)).radius = radius(index);
   %   chans(indices(index)).labels = int2str(index);
   %end;
   %figure; topoplot([],chans, 'style', 'blank', 'electrodes', 'labelpoint');

   [sph_phi, sph_theta] = topo2sph([angle radius], 2);
   for index = 1:length(indices)
      chans(indices(index)).sph_theta  = sph_theta(index);
      chans(indices(index)).sph_phi    = sph_phi  (index);
   end
case 'sphbesa2topo'
   chans = convertlocs(chans, 'sphbesa2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2topo', varargin{:}); % search for spherical coords
case 'sphbesa2cart'
   chans = convertlocs(chans, 'sphbesa2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2cart', varargin{:}); % search for spherical coords
case 'sphbesa2all'
   chans = convertlocs(chans, 'sphbesa2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2all', varargin{:}); % search for spherical coords
case 'cart2topo'
   chans = convertlocs(chans, 'cart2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2topo', varargin{:}); % search for spherical coords
case 'cart2sphbesa'
   chans = convertlocs(chans, 'cart2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2sphbesa', varargin{:}); % search for spherical coords
case 'cart2sph'
    if verbose
        disp('WARNING: If XYZ center has not been optimized, optimize it using Edit > Channel Locations');
    end
    X  = {chans.X};
    Y  = {chans.Y};
    Z  = {chans.Z};
    indices = find(~cellfun('isempty', X));
    [th, phi, radius] = cart2sph( [ X{indices} ], [ Y{indices} ], [ Z{indices} ]);
	for index = 1:length(indices)
		 chans(indices(index)).sph_theta     = th(index)/pi*180;
		 chans(indices(index)).sph_phi       = phi(index)/pi*180;
		 chans(indices(index)).sph_radius    = radius(index);
	end
case 'cart2all'
   chans = convertlocs(chans, 'cart2sph', varargin{:}); % search for spherical coords
   chans = convertlocs(chans, 'sph2all', varargin{:}); % search for spherical coords
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [model] = USRtrain(fea, options)
% USRtrain: Training Unsupervised Spectral Regression Model
%
%       [model] = USRtrain(fea, options)
%
%             Input:
%
%               fea     - data matrix. Each row is a data point.
%           options     - Struct value in Matlab. The fields in options
%                         that can be set:
%
%                      W       -  Affinity matrix. You can either call
%                                 "constructW" to construct the W, or
%                                 construct it by yourself.
%                                 If W is not provided, USRtrain will
%                                 build a k-NN graph with Heat kernel
%                                 weight, where k is a prameter.
%
%                      k       -  The parameter for k-NN graph (Default is 5)
%                                 If W is provided, this k will be ignored.
%
%                  ReducedDim  -  The number of reduced dimensions.
%                                 Default ReducedDim = 30.
%
%                  ReguType    -  'Ridge': L2-norm regularizer (default)
%                                 'Lasso': L1-norm regularizer
%                  ReguAlpha   -  regularization paramter for L2-norm regularizer
%                                 Default 0.1
%                  ReguGamma   -  regularization paramter for L1-norm regularizer
%                                 Default 0.05
%                  LASSOway    -  'LARs': use LARs to solve the LASSO
%                                 problem. You need to specify the
%                                 cardinality requirement in LassoCardi.
%                                 'SLEP': use SLEP to solve the LASSO
%                                 problem. Please see http://www.public.asu.edu/~jye02/Software/SLEP/
%                                  for details on SLEP. (The Default)
%
%                 bCenter = 0 | 1  whether to center the data. (In some
%                                  cases, e.g., text categorization, The
%                                  data is very spase, centering the data
%                                  will destroy the sparsity and consume
%                                  too much memory. In this case, bCenter
%                                  should be set to 0)
%                                   Default: 1
%
%
%             Output:
%               model   -  used for USRtest.m
%
%
%    Examples:
%
%
%
% See also SR, SR_caller
%
%Reference:
%
%   [1] Deng Cai, Xiaofei He, Wei Vivian Zhang, Jiawei Han, "Regularized
%   Locality Preserving Indexing via Spectral Regression", Proc. 2007 ACM
%   Int. Conf. on Information and Knowledge Management (CIKM'07), Lisboa,
%   Portugal, Nov. 2007.
%
%   [2] Deng Cai, "Spectral Regression: A Regression Framework for
%   Efficient Regularized Subspace Learning", PhD Thesis, Department of
%   Computer Science, UIUC, 2009.
%
%   version 3.0 --Jan/2012
%   version 2.0 --December/2011
%   version 1.0 --May/2006
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

if ~exist('options','var')
    options = [];
end

bCenter = 1;
if isfield(options,'bCenter')
    bCenter = options.bCenter;
end

k = 5;
if isfield(options,'k')
    k = options.k;
end

ReducedDim = 30;
if isfield(options,'ReducedDim')
    ReducedDim = options.ReducedDim;
end

if ~isfield(options,'ReguType')
    options.ReguType = 'Ridge';
end

LARs = false;
switch lower(options.ReguType)
    case {lower('Ridge')}
        if ~isfield(options,'ReguAlpha')
            options.ReguAlpha = 0.1;
        end
    case {lower('Lasso')}
        if isfield(options,'ReguAlpha') && options.ReguAlpha > 0
            options.RidgeAlpha = options.ReguAlpha;
            options.ReguType = 'RidgeLasso';
        end
        if isfield(options,'ReguGamma')
            options.ReguAlpha = options.ReguGamma;
        else
            options.ReguAlpha = 0.05;
        end
        if ~isfield(options,'LASSOway')
            options.LASSOway = 'SLEP';
        end

        if strcmpi(options.LASSOway,'LARs')
            LARs = true;
            if ~isfield(options,'LassoCardi')
                options.LassoCardi = 10:10:50;
            end
        end
    otherwise
        error('ReguType does not exist!');
end

nSmp=size(fea,1);
% Graph construction
if isfield(options,'W')
    W = options.W;
else
    Woptions.k = k;
    if nSmp > 3000
        tmpD = EuDist2(fea(randsample(nSmp,3000),:));
    else
        tmpD = EuDist2(fea);
    end
    Woptions.t = mean(mean(tmpD));
    W = constructW(fea,Woptions);
end

% Respnse generation
Y = Eigenmap(W,ReducedDim);

% Projection learning
if bCenter
    sampleMean = mean(fea);
    fea = (fea - repmat(sampleMean,nSmp,1));
end

[model.projection, LassoCardi] = SR(options, Y, fea);

model.LARs = LARs;
model.LassoCardi = LassoCardi;

model.TYPE = 'USR';
model.options = options;

author: ZJULearning

repo: MatlabFunc


text:
function [esTheta, esPhi] = mieScatteredField(An, Bn, theta, phi, frequency)

% Compute the complex-value scattered electric far field of a sphere
% using pre-computed coefficients An and Bn. Based on the treatment in:
%
% Ruck, et. al. "Radar Cross Section Handbook", Plenum Press, 1970.
%
% The incident electric field is in the -z direction (theta = 0) and is
% theta-polarized. The time-harmonic convention exp(jwt) is assumed, and
% the Green's function is of the form exp(-jkr)/r.
%
% Inputs:
%   An: Array of Mie solution constants
%   Bn: Array of Mie solution constants
%       (An and Bn should have the same length)
%   theta: Scattered field theta angle (radians)
%   phi: Scattered field phi angle (radians)
% Outputs:
%   esTheta: Theta-polarized electric field at the given scattering angles
%   esPhi: Phi-polarized electric field at the given scattering angles
%
%   Output electric field values are normalized such that the square of the
%   magnitude is the radar cross section (RCS) in square meters.
%
%   Author: Walton C. Gibson, email: kalla@tripoint.org

% speed of light
c = 299792458.0;

% wavenumber
k = 2.0*pi*frequency/c;

sinTheta = abs(sin(theta)); % note: theta only defined from from 0 to pi
cosTheta = cos(theta); % ok for theta > pi

% first two values of the Associated Legendre Polynomial
plm(1) = -sinTheta;
plm(2) = -3.0*sinTheta*cosTheta;

S1 = 0.0;
S2 = 0.0;

p = plm(1);

% compute coefficients for scattered electric far field
for iMode = 1:length(An)

        % derivative of associated Legendre Polynomial
    if abs(cosTheta) < 0.999999
        if iMode == 1
            dp = cosTheta*plm(1)/sqrt(1.0 - cosTheta*cosTheta);
        else
            dp = (iMode*cosTheta*plm(iMode) - (iMode + 1)*plm(iMode - 1))/sqrt(1.0 - cosTheta*cosTheta);
        end
    end

    if abs(sinTheta) > 1.0e-6
        term1 = An(iMode)*p/sinTheta;
        term2 = Bn(iMode)*p/sinTheta;
    end

    if cosTheta > 0.999999
        % Ruck, et. al. (3.1-12)
        val = ((i)^(iMode-1))*(iMode*(iMode+1)/2)*(An(iMode) - i*Bn(iMode));
        S1 = S1 + val;
        S2 = S2 + val;
    elseif cosTheta < -0.999999
        % Ruck, et. al. (3.1-14)
        val = ((-i)^(iMode-1))*(iMode*(iMode+1)/2)*(An(iMode) + i*Bn(iMode));
        S1 = S1 + val;
        S2 = S2 - val;
    else
        % Ruck, et. al. (3.1-6)
        S1 = S1 + ((i)^(iMode+1))*(term1 - i*Bn(iMode)*dp);
        % Ruck, et. al. (3.1-7)
        S2 = S2 + ((i)^(iMode+1))*(An(iMode)*dp - i*term2);
    end

    % recurrence relationship for next Associated Legendre Polynomial
    if iMode > 1
        plm(iMode + 1) = (2.0*iMode + 1)*cosTheta*plm(iMode)/iMode - (iMode + 1)*plm(iMode - 1)/iMode;
    end
    p = plm(iMode + 1);
end

% complex-value scattered electric far field, Ruck, et. al. (3.1-5)
esTheta = S1*cos(phi);
esPhi = -S2*sin(phi);

% normalize electric field so square of magnitude is RCS in square meters
esTheta = esTheta*sqrt(4.0*pi)/k;
esPhi = esPhi*sqrt(4.0*pi)/k;

return


text:
function Xd = Pendulum_Cart_Dynamics_Forced(X,F,P_Dyn)
% Xd = Pendulum_Cart_Dynamics_Forced(X,F,P_Dyn)
%
% FUNCTION:
%   This function is used to compute the dynamics of the pendulum cart
%   system. The inputs are the system state, actuator force, and paramters.
%
% INPUTS:
%
%   X = [4 x N] state matrix: [x;v;th;w]
%   F = [1 x N] actuator force (applied horizontally to cart)
%   P_Dyn = a parameter struct with fields:
%       M = cart mass
%       m = bob mass
%       L = pendulum length
%       g = gravitational acceleration
%
% OUTPUTS:
%   dX = [4 x N] derivative of the state matrix [dx;dv;dth;dw]
%
% NOTES:
%   The equations of motion are found using a cartesian force balance for
%   each rigid body
%
%   The cart is modelled as a point mass, and the pendulum is modelled as a
%   point mass at the end of a slender rod. There is no friction in the
%   system.
%
%   Parts of the derivation are shown below in comments. This code has been
%   partially optimized, so some of the comments show intermediate steps
%   that were bypassed.
%

%Define Parameters
m = P_Dyn.m;    %Pendulum point mass
M = P_Dyn.M;    %Cart mass
L = P_Dyn.L;    %Pendulum length (to point mass)
g = P_Dyn.g;    %Gravity

%Define input states:
% % x = X(1);    %Horizontal position of the cart
v = X(2,:);      %Horizontal velocity of the cart
th = X(3,:);     %Angle of the pendulum: 0 = straight up, pi = straight down
w = X(4,:);      %Angular rate of the pendulum

%Evaluate sines and cosines
Sin = -sin(th);
Cos = -cos(th);
One = ones(size(Sin));

%"Constant terms"
% % App_Forces = [...
% %     -F;
% %     M*g;
% %     -m*L*w^2*Sin;
% %     m*g + m*L*w^2*Cos];

% % %Matrix for of coeffieients:
% %     % Coeff * [Tension; Normal; Trans. Acc; Rot Acc] =
% %     % Applied_Forces( X0, Y0, X1, Y1 )
% % Coeff = [...
% %     sin(th)     0   -M  0;
% %     -cos(th)    1   0   0;
% %     -sin(th)    0   -m  -m*L*cos(th);
% %     cos(th)     0   0   -m*L*sin(th)];

% Analytic Inverse of Coeff via Mathematica:
Det = m*L*(M*Cos.^2 + (m+M)*Sin.^2);
inv_Det = 1./Det;
% % Adj = [...
% %     m^2*L*Sin       0   -m*M*L*Sin          m*M*L*Cos;
% %     m^2*L*Sin*Cos   1   -m*M*L*Sin*Cos      m*M*L*Cos^2;
% %     -m*L          	0   -m*L*Sin^2          m*L*Cos*Sin;
% %     m*Cos        	0   -M*Cos              -(M+m)*Sin];

% % Inv = inv_Det*Adj;
% % UnKnowns = Inv*App_Forces; %Solve linear system of equations

% % % % Tension = Unknowns(1);
% % % % Normal = UnKnowns(2);
% % vd = UnKnowns(3);
% % wd = UnKnowns(4);

%Unknowns (By Hand)
% % UnKnowns = inv_Det*[...
% %     %-F*               M*g*       -m*L*w^2*Sin*         (m*g + m*L*w^2*Cos)*
% %     -F*m^2*L*Sin     + 0   +  m*L*w^2*Sin*m*M*L*Sin     + (m*g + m*L*w^2*Cos)*m*M*L*Cos ;    %Tension
% %     -F*m^2*L*Sin*Cos + M*g +  m*L*w^2*Sin*m*M*L*Sin*Cos + (m*g + m*L*w^2*Cos)*m*M*L*Cos^2 ;  %Normal Force
% %     F*m*L            + 0   +  m*L*w^2*Sin*m*L*Sin^2     + (m*g + m*L*w^2*Cos)*m*L*Cos*Sin ;  %Translational acceleration
% %     -F*m*Cos         + 0   +  m*L*w^2*Sin*M*Cos         +  -(m*g + m*L*w^2*Cos)*(M+m)*Sin ]; %Angular Acceleration

%Broken Out UnKnowns (By Hand)
vd = (F*m*L  + 0 +  m*L*w.^2.*Sin*m*L.*Sin.^2 + (m*g*One + m*L*w.^2.*Cos).*m*L.*Cos.*Sin).*inv_Det;%Translational acceleration
wd = (-F.*m.*Cos  + 0  + m*L*w.^2.*Sin.*M.*Cos  +  -(m*g*One + m*L*w.^2.*Cos)*(M+m).*Sin).*inv_Det;%Angular Acceleration

%Defined derivatives:
xd = v;
thd = w;

%Express as a vector
Xd = [xd;vd;thd;wd];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function ok = mrInit(varargin)
% Initialize a mrVista session. Replacement for mrInitRet, mrvInit and
% mrInit2.  Moved from mrRSTools (formerly mrVista2).
%
% USAGE: mrInit should be called when the current directory is the one
% you wish to initialize. It can be called one of three ways:
%
%	(1) Interactively: type 'mrInit' without arguments, and a series of
%		dialogues ('mrInitGUI') will collect the initialization parameters.
%
%	(2) With a params struct: mrInit(params) will initialize without user
%		feedback. params is a struct with fields described below.
%
%	(3) mrInit('param', [value], 'param', [value]...) will also initialize
%		without user feedback. The 'param' strings are fields in the param
%		structs, the values are the value for that parameters. Any
%		non-specicified parameters will use the default values.
%
%
% INPUTS:
%   functionals: path string, or cell array of path strings, to the input
%   functional files. Defaults to looking in Raw/Pfiles/P*.mag for Lucas
%   P-mags. Input functionals can be any format readable by mrLoad,
%   including ANALYZE, NIFTI, or DICOM images.
%
%   inplane: path string to the input co-registered inplane anatomy.
%   Defaults to 'Raw/Anatomy/Inplane/I*.dcm' (again consistent w/ Lucas
%   Center conventions.) Can be any format readable by mrLoad.
%
%   params: structure which specifies the organization of the 'Original'
%   data type in the new session. See mrInitDefaultParams for a full
%   description of the parameters. mrInitGUI calls a series of dialogs to
%   get these parameters interactively.
%
% OUTPUTS:
%     The main output is a new mrVista session, with the following files
%     and directory structure:
%       mrSESSION.mat: contains the variables:
%           * mrSESSION: contains info on the input inplanes and
%           functionals.
%           * dataTYPES: info on groups of scans. Initialized only to the
%           default first group, 'Original', which derives from the set of
%           functionals. Analysis parameters for traveling-wave / event
%           analyses are also stored here.
%       Inplane/ directory: contains:
%           * anat.mat: anatomy file. Soon, this should be replaced with a
%           more robust format, like anat.nii.gz (NIFTI) or anat.img
%           (ANALYZE).
%           * Original/: 'data directory' for the Original data type.
%               * Original/TSeries/: directory for time series data from the
%                 functionals. Again, this is an older format: each scan
%                 has a directory ('/Scan2/'), and contains one MATLAB file
%                 for each slice ('tSeries4.mat'), which has the variable
%                 'tSeries', format [nFrames x nVoxelsInSlice].
%           * If running motion compensation, will also create new data
%           types: 'MotionComp' for within-scan, and 'BwScansMotionComp'
%           for between scan. If running both, you can choose to do either
%           between or within-scan first, then the other on top of that.
%           (You may want to remove the intermediate within-scan data type
%           to save space, if it looks like things went well.)
%
%   Also returns the variable ok, which is 1 if everything succeeded and 0
%   otherwise.
%
% EXAMPLE:
%	Suppose we have a session with four functional files, named 'func1.img'
%	to 'func4.img', and an inplane file, 'inplane.img'. All files are
%	single ANALYZE-format images.
%
%	params.inplane = 'inplane.img';
%	for i=1:5, params.functionals{i} = sprintf('func%i.img', i); end
%	mrInit(params);
%
%
% ras, 09/2006.

% PROGRAMMING NOTES:
% This code is intended to improve upon mrInitRet, addressing a few
% concerns:
%
% (1) More input file neutral. While mrInitRet was largely built to use
%     Lucas Center P.mag files, mrInit is intended to allow any files
%     readable by the mrVista 2 mrLoad functions. Specifically,
%     ANALYZE files are supported. Also, it's no longer required that
%     you have the Inplane files be DICOM files in Raw/Anatomy/Inplane.
%     Although consistent naming is recommended, the input anatomy file can
%     also be any format or location.
%
% (2) Scriptable. mrInitRet provided a nice 'wizard' of dialogs, but was
%     extremely hard to automate. This code is intended to be more modular,
%     and not require any user dialogs to run. If the optional 'params'
%     argument is not provided, it will provide a set of dialogs (mrInitGUI)
%     to get them. But otherwise, it runs without further input from the
%     user.
%
% (3) Saves more info from the header files. Things like the corners of the
%     inplane in scanner coordinates were not extracted previously, making
%     it difficult to e.g. label left from right. These are now saved in
%     the mrSESSION file.
%
% (4) Free of several of the bugs that have been cropping up lately, e.g.
%     in EditSession.
%
% (5) Run Motion Compensation / Slice Timing Correction at the front end,
%    so you can leave a session intitializing and have it also do these
%    steps overnight.
% ras, 07/23/07: debugged / slight overhaul to be cleaner, rely more on
% mrLoad and mrSave to actually set the mrSESSION fields.
% ras, 05/08: fixed several bugs in the GUI processing stream for dealing
% with file patterns; the main code also properly tests that the specified
% session directory exists, and cleans up the workspace when it finishes.
%
% Modified by DY on 09/2008 to allow for within then between or between
% then within MC. Changed comments here and in mrInitDefaultParams.m to
% reflect this. Email DY if you would like to see an example wrapper script
% that sets all params.

%%%%% (0) ensure all input parameters are specified
if nargin==0
    % get params interactively
    [params, ok] = mrInitGUI;
    if ~ok, disp('mrInit Aborted.'); return; end

elseif length(varargin)==1 && isstruct(varargin{1})
    % params struct entered
    params = mrInitDefaultParams;
    params = mergeStructures(params, varargin{1});

else
    params = mrInitDefaultParams;

    % param/value pairs specified -- parse
    for i = 1:2:length(varargin)
        params.(varargin{i}) = varargin{i+1};
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (1) cd to session dir; init empty session %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params.startTime = datestr(now);
fprintf('***** [%s] Initializing Session %s ***** (%s)\n', mfilename, ...
    params.sessionCode, params.startTime);
ensureDirExists(params.sessionDir);
callingDir = pwd;
cd(params.sessionDir);

mrGlobals; %Declared here since we want HOMEDIR to be sessPath

initEmptySession; %Replace this save and then load of mrSESSION with that variable simply passed
% from one to the other
load mrSESSION mrSESSION dataTYPES
mrSESSION = sessionSet(mrSESSION,'description', params.description);
mrSESSION = sessionSet(mrSESSION,'sessionCode',params.sessionCode);
mrSESSION = sessionSet(mrSESSION,'subject',params.subject);
mrSESSION = sessionSet(mrSESSION,'comments',params.comments);
mrSESSION = sessionSet(mrSESSION,'Inplane Path',params.inplane); %Populates the mrSESSION inplane path var
mrSESSION = sessionSet(mrSESSION,'alignment',params.alignment);
mrSESSION = sessionSet(mrSESSION,'Version','2.1');

save mrSESSION mrSESSION -append;
save mrInit_params params   % stash the params in case we crash

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (2) Load functional data, if it exists      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isfield(params,'functionals') && ~isempty(params.functionals)
    if ischar(params.functionals), params.functionals = {params.functionals}; end

    for scan = 1:length(params.functionals)

        func.path = params.functionals{scan};

        %Read in the nifti to the tS struct, then apply the same transform as
        %the inplane data. Then, transfer the necessary components to the
        %local func struct.
        tS = niftiRead(func.path);
        tS = niftiApplyAndCreateXform(tS,'Inplane');

        func.hdr = rmfield(tS, 'data');

        %Store the orientation of the functional data. We will need the
        %anatomical Inplane to have the same orientation.
        [~, func.orientation] = niftiCreateXform(tS,'Inplane');

        %Need to move over:
        %Data
        func.data = niftiGet(tS,'Data');
        %Dims
        func.dims = niftiGet(tS, 'Dim');
        %PixDims
        func.pixdims = niftiGet(tS, 'Pix Dim', 'xyz_units', 'mm', 'time_units', 's');

        % select keepFrames if they're provided
        if isfield(params, 'keepFrames') && ~isempty(params.keepFrames)
            %Put keepFrames into func so that we can save it into mrSESSION
            func.keepFrames = params.keepFrames(scan,:);
        else
            %We will need to create keepFrames if it doesn't exist
            func.keepFrames = [0 -1];
        end

        % assign annotation if it's provided
        if length(params.annotations) >= scan && ...
                ~isempty(params.annotations{scan})
            func.name = params.annotations{scan};
        else
            func.name = [];
        end

        %This call updates dataTYPES as well
        %We will no longer call mrSave
        mrInitInplaneTseries(func,scan);

    end

end

fprintf('[%s]: Finished initializing mrVista session. \t(%s)\n', ...
    mfilename, datestr(now));

%%%%%%%%%%%%%%%%%%%%%%%
% (3)  VOLUME anatomy %
%%%%%%%%%%%%%%%%%%%%%%%
if isfield(params,'vAnatomy') && ~isempty(params.vAnatomy)
    setVAnatomyPath(params.vAnatomy);
end

%%%%%%%%%%%%%%%%%%%%%%%
% (4)  pre-processing %
%%%%%%%%%%%%%%%%%%%%%%%
% update dataTYPES as needed.

INPLANE{1} = initHiddenInplane;

% assign coherence analysis parameters
for s = cellfind(params.coParams)
    coParams = params.coParams{s};
    for f = fieldnames(coParams)'
        dataTYPES.blockedAnalysisParams(s).(f{1}) = coParams.(f{1});
    end
end

% assign GLM analysis parameters
for scan = cellfind(params.glmParams)
    er_setParams(INPLANE{1}, params.glmParams{scan}, scan);
end

% slice timing correction
if params.sliceTimingCorrection==1
    if isfield(params, 'sliceOrder')
        mrSESSION = sessionSet(mrSESSION, 'sliceorder', params.sliceOrder);
        saveSession;
    end
    INPLANE{1} = AdjustSliceTiming(INPLANE{1}, 1:length(sessionGet(mrSESSION, 'functionals')));
    INPLANE{1} = selectDataType(INPLANE{1}, 'Timed');
end

% motion compensation
if params.motionComp > 1
    switch params.motionComp
        case 1, % between scans only
            newDataType = 'BwScansMotionComp';
            if ~existDataType(newDataType), addDataType(newDataType); end
            hI = initHiddenInplane(newDataType);
            [INPLANE{1}, ~] = betweenScanMotComp(INPLANE{1}, hI, params.motionCompRefScan);
            INPLANE{1} = selectDataType(INPLANE{1}, newDataType);

        case 2, % within scans only
            INPLANE{1} = motionCompSelScan(INPLANE{1}, 'MotionComp', ...
                [], params.motionCompRefFrame);

        case 3, % both between and within scans
            motionCompNestaresFull(INPLANE{1}, [], ...
                params.motionCompRefScan, ...
                params.motionCompRefFrame);
            INPLANE{1} = selectDataType(INPLANE{1}, 'MotionComp');
        case 4, % both between and within scans
            motionCompNestaresWithin1st(INPLANE{1}, [], params.motionCompRefScan, params.motionCompRefFrame);
            mcdtName=['MotionComp_RefScan' num2str(params.motionCompRefScan)];
            INPLANE{1} = selectDataType(INPLANE{1}, mcdtName);

    end
end

% group scans, assign parfiles, apply GLM
for scan = cellfind(params.parfile)
    er_assignParfilesToScans(INPLANE{1}, scan, params.parfile(scan));
end

if ~isempty(params.scanGroups)
    for ii = 1:length(params.scanGroups)
        INPLANE{1} = er_groupScans(INPLANE{1}, params.scanGroups{ii}, 2);

        if params.applyGlm==1
            glmParams = er_getParams(INPLANE{1}, params.scanGroups{ii}(1));
            applyGlm(INPLANE{1}, [], params.scanGroups{ii}, glmParams);
        end
    end
end

% apply coherence analysis
if any(params.applyCorAnal > 0)
    INPLANE{1} = computeCorAnal(INPLANE{1}, params.applyCorAnal, 1);
end

% All tasks are now complete
saveSession;
fprintf('***** [%s] Finished Initializing Session %s (%s)*****\n', ...
    mfilename, mrSESSION.sessionCode, datestr(now));
mrvCleanWorkspace;
ok = 1;

cd(callingDir);

return
% /---------------------------------------------------------------------/ %

% /-----------------------------------------------------------------/ %
function mrInitInplaneTseries(mr, scan)

%Taken from mrSave, with the 'save' functionality removed, the rest of the
%data transfer functionality has been retained, however

mrGlobals;

% set header info in mrSESSION.functionals, dataTYPES.scanParmas
if notDefined('scan') || isempty(scan), scan = length(mrSESSION.functionals) + 1; end %#ok<NODEF>
f.PfileName = mr.path;
f.totalFrames = mr.dims(4);

f.firstName = '';  f.lastName = '';
if checkfields(mr, 'info', 'subject')
    sp = strfind(' ', mr.info.subject);
    if ~isempty(sp) % space in name
        f.firstName = mr.info.subject( 1:(sp(1)-1) );
        f.lastName = mr.info.subject( (sp(1)+1):end );
    else
        f.firstName = mr.info.subject;
    end
end

f.date = ''; f.time = '';
if checkfields(mr, 'info', 'date'), f.date = mr.info.date; end
if checkfields(mr, 'info', 'time'), f.time = mr.info.time; end

f.junkFirstFrames = 0; %This always appears to be 0. perhaps remove it?

if mr.keepFrames(2) == -1,
    % if 2nd column of keepframes is -1, keep all frames after drop frames
    nFrames = mr.dims(4) - mr.keepFrames(1);
else
    % if 2nd column of keepframes is +n, keep n frames after drop frames
    nFrames = mr.keepFrames(2);
end

f.nFrames   = nFrames;
f.slices    = 1:mr.dims(3);
f.fullSize  = mr.dims(1:2);
f.cropSize  = mr.dims(1:2);
f.crop      = [1 1; mr.dims(1:2)];
f.voxelSize = mr.pixdims(1:3);
f.effectiveResolution = mr.pixdims(1:3);
f.keepFrames = mr.keepFrames; %Keep Frames will now be udpated in both mrSESSION and dataTYPES
if checkfields(mr, 'info', 'effectiveResolution')
    f.effectiveResolution = mr.info.effectiveResolution;
end
f.framePeriod = mr.pixdims(4);
f.reconParams = mr.hdr;

if scan==1
    mrSESSION = sessionSet(mrSESSION, 'Functionals', f);
    mrSESSION = sessionSet(mrSESSION, 'Functional Orientation', mr.orientation);
else
    mrSESSION = sessionSet(mrSESSION, 'Functionals', ...
        mergeStructures(sessionGet(mrSESSION, 'Functionals', scan-1), f), scan);
end

% Default params, initializing the parameters in dataTYPES

% Copy one field at a time, so we don't get type-mismatch errors.

% scan params
srcScanParams = scanParamsDefaults(mrSESSION, scan, mr.name);
for f = fieldnames(srcScanParams)'
    dataTYPES(1).scanParams(scan).(f{1}) = srcScanParams.(f{1});
end

% blocked analysis params
srcBlockParams = blockedAnalysisDefaults;
for f = fieldnames(srcBlockParams)'
    dataTYPES(1).blockedAnalysisParams(scan).(f{1}) = ...
        srcBlockParams.(f{1});
end

% event analysis params
srcEventParams = er_defaultParams;
for f = fieldnames(srcEventParams)'
    dataTYPES(1).eventAnalysisParams(scan).(f{1}) = ...
        srcEventParams.(f{1});
end

saveSession

return
% /-----------------------------------------------------------------/ %

% /-----------------------------------------------------------------/ %
function params = scanParamsDefaults(mrSESSION, scan, annotation)
% Default scan parameters for a new scan.
params.annotation   = annotation;
params.nFrames      = mrSESSION.functionals(scan).nFrames;
params.framePeriod  = mrSESSION.functionals(scan).framePeriod;
params.slices       = mrSESSION.functionals(scan).slices;
params.cropSize     = mrSESSION.functionals(scan).cropSize;
params.PfileName    = mrSESSION.functionals(scan).PfileName;
params.inplanePath  = mrSESSION.functionals(scan).PfileName;
params.keepFrames   = mrSESSION.functionals(scan).keepFrames;
params.parfile      = '';
params.scanGroup            = sprintf('Original: %i',scan);
return
% /-----------------------------------------------------------------/ %

% /-----------------------------------------------------------------/ %
function params = blockedAnalysisDefaults
% Default values for the blocked analyses.
params.blockedAnalysis       = 1;
params.detrend               = 1;
params.inhomoCorrect         = 1;
params.temporalNormalization = 0;
params.nCycles               = 6;
return
% /-----------------------------------------------------------------/ %


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

%%
clear; close all; clc;

%% Syntax
% |Ad=dev(A);|

%% Description
% Computes the deviatoric part of the tensor A.
%
% See also |sph|.

%% Examples
%

A=rand(3,3)

Ad=dev(A)

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [M_final,shifts_g,template,options,col_shift] = normcorre_batch_even(Y,options,template)

% online motion correction through DFT subpixel registration
% Based on the dftregistration.m function from Manuel Guizar and Jim Fienup

% INPUTS
% Y:                Input data, can be already loaded in memory as a 3D
%                   tensor, a memory mapped file, or a pointer to a tiff stack
% options:          options structure for motion correction (optional, rigid registration is performed if not provided)
% template:         provide template (optional)

% OUTPUTS
% M_final:          motion corrected data
% shifts_up:        upsampled shifts
% shifts:           originally calculated shifts
% template:         calculated template

%% first determine filetype

if isa(Y,'char')
    [~,~,ext] = fileparts(Y);
    ext = ext(2:end);
    if strcmpi(ext,'tif') || strcmpi(ext,'tiff')
        tiffInfo = imfinfo(Y);
        filetype = 'tif';
        T = length(tiffInfo);
        sizY = [tiffInfo(1).Height,tiffInfo(1).Width,T];
    elseif strcmpi(ext,'mat')
        filetype = 'mem';
        Y = matfile(Y,'Writable',true);
        details = whos(Y);
        var_sizes = [details.bytes];
        [~,var_ind] = max(var_sizes);
        var_name = details(var_ind).name;
        sizY = size(Y,var_name);
        T = sizY(end);
    elseif strcmpi(ext,'hdf5') || strcmpi(ext,'h5')
        filetype = 'hdf5';
        fileinfo = hdf5info(Y);
        data_name = fileinfo.GroupHierarchy.Datasets.Name;
        sizY = fileinfo.GroupHierarchy.Datasets.Dims;
        T = sizY(end);
    elseif strcmpi(ext,'raw')
        filetype = 'raw';
        fid = fopen(Y);
        FOV = [options.d1,options.d2];
        bitsize = options.bitsize;
        imsize = FOV(1)*FOV(2)*bitsize;                                                   % Bit size of single frame
        current_seek = ftell(fid);
        fseek(fid, 0, 1);
        file_length = ftell(fid);
        fseek(fid, current_seek, -1);
        T = file_length/imsize;
        sizY = [FOV,T];
        fclose(fid);
    elseif strcmpi(ext,'avi')
        filetype = 'avi';
        sizY = size(read_file(Y));
        FOV = sizY(1:2);
        T = sizY(end);
    end
elseif isobject(Y)
    filetype = 'mem';
    var_name = 'Y';
    sizY = size(Y,var_name);
    T = sizY(end);
else % array loaded in memory
    filetype = 'mat';
    sizY = size(Y);
    T = sizY(end);
end

nd = length(sizY)-1;                          % determine whether imaging is 2d or 3d
sizY = sizY(1:nd);
otherdims = repmat({':'},1,nd);
%% set default parameters if not present

if ~exist('options','var') || isempty(options)
    options = NoRMCorreSetParms('d1',sizY(1),'d2',sizY(2));
    if nd > 2; options.d3 = sizY(3); end
end

memmap = options.memmap;
grid_size = options.grid_size;
mot_uf = options.mot_uf;
min_patch_size = options.min_patch_size;
min_diff = options.min_diff;
overlap_pre = options.overlap_pre;
overlap_post = options.overlap_post;
upd_template = options.upd_template;
bin_width = options.bin_width;
buffer_width = options.buffer_width;
max_dev_g = options.max_dev;
init_batch = options.init_batch;
us_fac = options.us_fac;
method = options.method;
filename = options.mem_filename;
iter = options.iter;
add_value = options.add_value;
max_shift = options.max_shift;
if strcmpi(options.boundary,'nan')
    fill_value = NaN;
else
    fill_value = add_value;
end
print_msg = options.print_msg;

while mod(T,bin_width) == 1
    if T == 1
        error('Movie appears to have only one frame. Use the function normcorre instead')
    end
    bin_width = bin_width + 1;
end

%% first check for offset due to bi-directional scanning

if options.correct_bidir && isempty(options.col_shift)
    col_shift = correct_bidirectional_offset(Y,options.nFrames,options.bidir_us);
elseif ~isempty(options.col_shift)
    col_shift = options.col_shift;
else
    col_shift = 0;
end
options.col_shift = col_shift;
if col_shift
    if print_msg; fprintf('Offset %1.1d pixels due to bidirectional scanning detected. \n',col_shift); end
    if strcmpi(options.shifts_method,'fft')
        options.shifts_method = 'cubic';
        if print_msg; fprintf('Cubic shifts will be applied. \n'); end
    end
end
%% read initial batch and compute template

init_batch = min(T,init_batch);
interval = ceil(T/2-init_batch/2+1):floor(T/2+init_batch/2);
switch filetype
    case 'tif'
        Y_temp = read_file(Y,interval(1),init_batch,[],tiffInfo);
    case 'hdf5'
        Y_temp = read_file(Y,interval(1),init_batch);
    case 'avi'
        Y_temp = read_file(Y,interval(1),init_batch);
    case 'mem'
        Y_temp = Y.(var_name)(otherdims{:},interval);
    case 'mat'
        Y_temp = Y(otherdims{:},interval);
    case 'raw'
        Y_temp = read_raw_file(Y,interval(1),init_batch,FOV,bitsize);
end
data_type = class(Y_temp);
Y_temp = single(Y_temp);
use_proj = true;
if nargin < 3 || isempty(template)
    if print_msg; fprintf('Registering the first %i frames just to obtain a good template....',init_batch); end
    template_in = median(Y_temp,nd+1)+add_value;
    fftTemp = fftn(template_in);
    for t = 1:size(Y_temp,nd+1)
        if nd == 2
            [~,Greg] = dftregistration_min_max(fftTemp,fftn(Y_temp(:,:,t)),us_fac,-max_shift,max_shift,options.phase_flag);
        end
        if nd == 3
            [~,Greg] = dftregistration_min_max_3d(fftTemp,fftn(Y_temp(:,:,:,t)),us_fac,-max_shift,max_shift,options.phase_flag);
        end
        M_temp = real(ifftn(Greg));
        template_in = template_in*(t-1)/t + M_temp/t;
    end
    template_in = template_in + add_value;
    if print_msg; fprintf('..done. \n'); end
else
    template_in = single(template + add_value);
end

[d1,d2,d3,~] = size(Y_temp);
if nd == 2; d3 = 1; end
%% setup grids for patches
dim = [d1,d2,d3];
patches = construct_grid_even(grid_size,overlap_pre,dim,min_diff);
shifts_g = struct('shifts',cell(T,1),'shifts_up',cell(T,1),'diff',cell(T,1));
template_patches = split_frame(template_in,patches);
patch_size = size(template_patches);
n_patches = [length(unique(patches(:,1))),length(unique(patches(:,3))),length(unique(patches(:,5)))];
%%
%maxNumCompThreads(1);
temp_mat = template_in;

fftTemp =fft(fft(fft(template_patches,[],1),[],2),[],3);
fftTempMat = fftn(temp_mat);

if ~strcmpi(options.output_type,'mat')
    options.mem_batch_size = max(min(round(options.mem_batch_size/bin_width)*bin_width,T),1);
    mem_buffer = squeeze(zeros([dim,options.mem_batch_size],'single'));
end

switch lower(options.output_type)
    case 'mat'
        M_final = zeros([sizY,T],data_type);
    case 'memmap'
        M_final = matfile(filename,'Writable',true);
        if nd == 2; M_final.Y(d1,d2,T) = zeros(1,data_type); end
        if nd == 3; M_final.Y(d1,d2,d3,T) = zeros(1,data_type); end
        M_final.Yr(d1*d2*d3,T) = zeros(1,data_type);
    case {'hdf5','h5'}
         if exist(options.h5_filename,'file')
            [pathstr,fname,ext] = fileparts(options.h5_filename);
            new_filename = fullfile(pathstr,[fname,'_',datestr(now,30),ext]);
            warning_msg = ['File ',options.h5_filename,'already exists. Saving motion corrected file as',new_filename];
            warning('%s',warning_msg);
            options.h5_filename = new_filename;
        end
        M_final = options.h5_filename;
        if nd == 2
            h5create(options.h5_filename,['/',options.h5_groupname],[d1,d2,Inf],'Chunksize',[d1,d2,options.mem_batch_size],'Datatype',data_type);
        elseif nd == 3
            h5create(options.h5_filename,['/',options.h5_groupname],[d1,d2,d3,Inf],'Chunksize',[d1,d2,d3,options.mem_batch_size],'Datatype',data_type);
        end
    case {'tif','tiff'}
        M_final = ['motion corrected file has been saved as ', options.tiff_filename];
        opts_tiff.append = true;
        opts_tiff.big = true;
        if nd == 3
            error('Saving volumetric tiff stacks is currently not supported. Use a different filetype');
        end
    otherwise
        error('This filetype is currently not supported')
end

cnt_buf = 0;
if print_msg; fprintf('Template initialization complete.  Now registering all the frames with new template. \n'); end
%%

prevstr = [];
for it = 1:iter
    for t = 1:bin_width:T
        switch filetype
            case 'tif'
                Ytm = single(read_file(Y, t, min(t+bin_width-1,T)-t+1, [], tiffInfo));
            case 'avi'
                Ytm = single(read_file(Y, t, min(t+bin_width-1,T)-t+1));
            case 'hdf5'
                Ytm = single(h5read(Y,data_name,[ones(1,nd),t],[sizY(1:nd),min(t+bin_width-1,T)-t+1]));
            case 'mem'
                Ytm = single(Y.(var_name)(otherdims{:},t:min(t+bin_width-1,T)));
            case 'mat'
                Ytm = single(Y(otherdims{:},t:min(t+bin_width-1,T)));
            case 'raw'
                Ytm = single(read_raw_file(Y,t,min(t+bin_width-1,T)-t+1,FOV,bitsize));
        end

        Mf = zeros(size(Ytm),data_type);
        lY = size(Ytm,nd+1);
        shifts = struct('shifts',cell(lY,1),'shifts_up',cell(lY,1),'diff',cell(lY,1));
        for i = lY:-1:1
            Yt = Ytm(otherdims{:},i);
            future_results(i) = parfeval(@register_frame, 2, Yt,fftTempMat,fftTemp,patches,options);
        end
        for i = 1:lY
            [idx, shifts_temp, Mf_temp] = fetchNext(future_results);
            shifts(idx).shifts = shifts_temp;
            shifts(idx).shifts_up = shifts_temp;
            Mf(otherdims{:},idx) = Mf_temp;
        end

        shifts_g(t:min(t+bin_width-1,T)) = shifts;

        if ~strcmpi(options.output_type,'mat')
            rem_mem = rem(t+lY-1,options.mem_batch_size);
            if rem_mem == 0; rem_mem = options.mem_batch_size; end
            mem_buffer(otherdims{:},rem_mem-lY+1:rem_mem) = cast(Mf,data_type);
        end
        if it == iter
            switch lower(options.output_type)
                case 'mat'
                    M_final(otherdims{:},t:min(t+bin_width-1,T)) = cast(Mf,data_type);
                case 'memmap'
                    if rem_mem == options.mem_batch_size || t+lY-1 == T
                        M_final.Y(otherdims{:},t+lY-rem_mem:t+lY-1) = mem_buffer(:,:,1:rem_mem);
                        M_final.Yr(:,t+lY-rem_mem:t+lY-1) = reshape(mem_buffer(1:d1*d2*d3*rem_mem),d1*d2*d3,rem_mem);
                    end
                case {'hdf5','h5'}
                    if rem_mem == options.mem_batch_size || t+lY-1 == T
                        if nd == 2; h5write(options.h5_filename,['/',options.h5_groupname],mem_buffer(:,:,1:rem_mem),[ones(1,nd),t+lY-rem_mem],[sizY(1:nd),rem_mem]); end
                        if nd == 3; h5write(options.h5_filename,['/',options.h5_groupname],mem_buffer(:,:,:,1:rem_mem),[ones(1,nd),t+lY-rem_mem],[sizY(1:nd),rem_mem]); end
                    end
                case {'tif','tiff'}
                    if rem_mem == options.mem_batch_size || t+lY-1 == T
                        saveastiff(cast(mem_buffer(:,:,1:rem_mem),data_type),options.tiff_filename,opts_tiff);
                    end
            end
        end
        if print_msg
            str = sprintf('%d out of %d frames registered, iteration %d out of %d..', t+lY-1, T, it, iter);
            refreshdisp(str, prevstr, t);
            prevstr=str;
        end
        % update template
        if upd_template
            cnt_buf = cnt_buf + 1;
            if strcmpi(method{2},'mean')
                new_temp = nanmean(Mf,nd+1);
            elseif strcmpi(method{2},'median')
                new_temp = nanmedian(x,nd+1);
            end
            if any(isnan(new_temp(:)))
                new_temp(isnan(new_temp)) =  template(isnan(new_temp));
            end
            if strcmpi(method{1},'mean')
                cnt = t/bin_width + 1;
                temp_mat = temp_mat*(cnt-1)/cnt + new_temp/cnt;
            elseif strcmpi(method{1},'median')
                ind_buffer = mod(cnt_buf,buffer_width) + buffer_width*(~mod(cnt_buf,buffer_width));
                buffer_med(otherdims{:},ind_buffer) = new_temp;
                temp_mat = nanmedian(buffer_med,nd+1);
            end
            template_patches = split_frame(temp_mat,patches);
            fftTemp = fft(fft(fft(template_patches,[],1),[],2),[],3);
            fftTempMat = fftn(temp_mat);
        end
    end

    if it == iter
        template = temp_mat - add_value;
    end
    if memmap
        M_final.shifts = shifts_g;
        M_final.template = template;
    end
end
if print_msg; fprintf('\n'); end
maxNumCompThreads('automatic');

author: flatironinstitute

repo: NoRMCorre


text:
function q = uminus(q)
% overload unitary minus

q.a = -q.a;
q.b = -q.b;
q.c = -q.c;
q.d = -q.d;