text:
%FITELLIPSE  Fits an ellipse around a set of 2D points
%
%     rct = cv.fitEllipse(points)
%     rct = cv.fitEllipse(points, 'OptionName',optionValue, ...)
%
% ## Input
% * __points__ Input 2D point set, stored in numeric array
%   (Nx2/Nx1x2/1xNx2) or cell array of 2-element vectors (`{[x,y], ...}`).
%   There should be at least 5 points to fit the ellipse.
%
% ## Output
% * __rct__ Output rotated rectangle struct with the following fields:
%   * __center__ The rectangle mass center `[x,y]`.
%   * __size__ Width and height of the rectangle `[w,h]`.
%   * __angle__ The rotation angle in a clockwise direction. When the angle is
%     0, 90, 180, 270 etc., the rectangle becomes an up-right rectangle.
%
% ## Options
% * __Method__ One of:
%   * __Linear__ Linear (LIN) conic fitting method. This is the default.
%   * __Direct__ Direct least square (DLS) method.
%   * __AMS__ Approximate mean square (AMS) method.
%
% ### Method = Linear
%
% The function calculates the ellipse that fits (in a least-squares sense) a
% set of 2D points best of all. It returns the rotated rectangle in which the
% ellipse is inscribed. The first algorithm described by [Fitzgibbon95] is
% used. Developer should keep in mind that it is possible that the returned
% ellipse/rotatedRect data contains negative indices, due to the data points
% being close to the border of the containing Mat element.
%
% ### Method = Direct
%
% The function calculates the ellipse that fits a set of 2D points.
% It returns the rotated rectangle in which the ellipse is inscribed.
% The Direct least square (Direct) method by [Fitzgibbon1999] is used.
%
% For an ellipse, this basis set is `chi = (x^2, x*y, y^2, x, y, 1)`, which is
% a set of six free coefficients `A^T = {A_xx, A_xy, A_yy, A_x, A_y, A_0}`.
% However, to specify an ellipse, all that is needed is five numbers; the
% major and minor axes lengths `(a,b)`, the position `(x_0,y_0)`, and the
% orientation `theta`. This is because the basis set includes lines,
% quadratics, parabolic and hyperbolic functions as well as elliptical
% functions as possible fits.
%
% The Direct method confines the fit to ellipses by ensuring that
% `4*A_xx*A_yy - A_xy^2 > 0`. The condition imposed is that
% `4*A_xx*A_yy - A_xy^2 = 1` which satisfies the inequality and as the
% coefficients can be arbitrarily scaled is not overly restrictive.
%
%     epsilon^2 = A^T * D^T * D * A
%     with A^T * C * A = 1
%     and C = [0 0 2 0 0 0; 0 -1 0 0 0 0; 2 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0]
%
% The minimum cost is found by solving the generalized eigenvalue problem.
%
%     D^T * D * A = lambda * (C) * A
%
% The system produces only one positive eigenvalue `lambda` which is chosen as
% the solution with its eigenvector `u`. These are used to find the
% coefficients:
%
%     A = sqrt(1 / (u^T * C * u)) * u
%
% The scaling factor guarantees that  `A^T * C * A = 1`.
%
% Note: If the determinant of `A` is too small, the method fallsback to
% 'Linear'.
%
% ### Method = AMS
%
% The function calculates the ellipse that fits a set of 2D points.
% It returns the rotated rectangle in which the ellipse is inscribed.
% The Approximate Mean Square (AMS) proposed by [Taubin1991] is used.
%
% For an ellipse, this basis set is `chi = (x^2, x*y, y^2, x, y, 1)`, which is
% a set of six free coefficients `A^T = {A_xx, A_xy, A_yy, A_x, A_y, A_0}`.
% However, to specify an ellipse, all that is needed is five numbers; the
% major and minor axes lengths `(a,b)`, the position `(x_0,y_0)`, and the
% orientation `theta`. This is because the basis set includes lines,
% quadratics, parabolic and hyperbolic functions as well as elliptical
% functions as possible fits.
%
% If the fit is found to be a parabolic or hyperbolic function then the
% 'Direct' method is used. The AMS method restricts the fit to parabolic,
% hyperbolic and elliptical curves by imposing the condition that
% `A^T * (D_x^T * D_x + D_y^T * D_y) * A = 1` where the matrices `Dx` and `Dy`
% are the partial derivatives of the design matrix `D` with respect to x and y.
% The matrices are formed row by row applying the following to each of the
% points in the set:
%
%     D(i,:)   = {x_i^2, x_i y_i, y_i^2, x_i, y_i, 1}
%     D_x(i,:) = {2*x_i, y_i, 0, 1, 0, 0}
%     D_y(i,:) = {0, x_i, 2*y_i, 0, 1, 0}
%
% The AMS method minimizes the cost function
%
%     epsilon^2 = (A^T * D^T * D * A) / (A^T * (D_x^T * D_x + D_y^T * D_y) * A^T)
%
% The minimum cost is found by solving the generalized eigenvalue problem.
%
%     D^T * D * A = lambda * (D_x^T * D_x + D_y^T * D_y) * A
%
% Note: If the determinant of `A` is too small, the method fallsback to
% 'Linear'.
%
% ## References
% [Fitzgibbon95]:
% > Andrew W Fitzgibbon and Robert B Fisher.
% > "A buyer's guide to conic fitting". In Proceedings of the 6th British
% > conference on Machine vision (Vol. 2), pages 513-522. BMVA Press, 1995.
% > [PDF](http://www.bmva.org/bmvc/1995/bmvc-95-050.pdf)
%
% [Fitzgibbon1999]:
% > Andrew Fitzgibbon, Maurizio Pilu, and Robert B. Fisher.
% > "Direct least square fitting of ellipses". IEEE Transactions on Pattern
% > Analysis and Machine Intelligence, 21(5):476-480, 1999.
% > [PDF](https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/ellipse-pami.pdf)
%
% [Taubin1991]:
% > Gabriel Taubin. "Estimation of planar curves, surfaces, and nonplanar
% > space curves defined by implicit equations with applications to edge and
% > range image segmentation". IEEE Transactions on Pattern Analysis and
% > Machine Intelligence, 13(11):1115-1138, 1991.
%
% See also: cv.minAreaRect
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function out = SB_BinaryStats(y,binaryMethod)
% SB_BinaryStats    Statistics on a binary symbolization of the time series
%
% Binary symbolization of the time series is a symbolic string of 0s and 1s.
%
% Provides information about the coarse-grained behavior of the time series
%
%---INPUTS:
% y, the input time series
%
% binaryMethod, the symbolization rule:
%         (i) 'diff': by whether incremental differences of the time series are
%                      positive (1), or negative (0),
%         (ii) 'mean': by whether each point is above (1) or below the mean (0)
%         (iii) 'iqr': by whether the time series is within the interquartile range
%                      (1), or not (0).
%
%---OUTPUTS:
% Include the Shannon entropy of the string, the longest stretches of 0s
% or 1s, the mean length of consecutive 0s or 1s, and the spread of consecutive
% strings of 0s or 1s.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check inputs, set defaults:
%-------------------------------------------------------------------------------
if nargin < 2 || isempty(binaryMethod)
    binaryMethod = 'diff';
end

%-------------------------------------------------------------------------------
% Binarize the time series:
%-------------------------------------------------------------------------------
yBin = BF_Binarize(y,binaryMethod);

N = length(yBin); % length of signal - 1 (difference operation)

%-------------------------------------------------------------------------------
% Stationarity of binarized time series:
%-------------------------------------------------------------------------------
% (cf. SB_MotifTwo for basic stats on binarized time series)

% Stationarity:
out.pupstat2 = sum(yBin(floor(end/2)+1:end) == 1)/sum(yBin(1:floor(end/2)) == 1);

%-------------------------------------------------------------------------------
% Consecutive string of ones / zeros (normalized by length)
%-------------------------------------------------------------------------------
difffy = diff(find([1;yBin;1]));
stretch0 = difffy(difffy ~= 1) - 1;

difffy = diff(find([0;yBin;0] == 0));
stretch1 = difffy(difffy ~= 1) - 1;

%-------------------------------------------------------------------------------
% pstretches
%-------------------------------------------------------------------------------
% Number of different stretches as proportion of the time-series length
out.pstretch1 = length(stretch1)/N;
% The following are trivially dependent on pstretch1:
% out.pstretch0 = length(stretch0)/N;
% out.pstretches = (length(stretch0)+length(stretch1))/N;

if isempty(stretch0) % all 1s (almost impossible to actually occur)
    out.longstretch0 = 0;
    out.longstretch0norm = 0;
    out.meanstretch0 = 0;
    out.meanstretch0norm = 0;
    out.stdstretch0 = NaN;
    out.stdstretch0norm = NaN;
else
    out.longstretch0 = max(stretch0); % longest consecutive stretch of zeros
    out.longstretch0norm = max(stretch0)/N; % longest consecutive stretch of zeros as proportion of time-series length
    out.meanstretch0 = mean(stretch0); % mean stretch of zeros
    out.meanstretch0norm = mean(stretch0)/N; % mean stretch of zeros
    out.stdstretch0 = std(stretch0); % standard deviation of stretch lengths of consecutive zeros
    out.stdstretch0norm = std(stretch0)/N; % standard deviation of stretch lengths of consecutive zeros
end

if isempty(stretch1) % all 0s (almost impossible to actually occur)
    out.longstretch1 = 0;
    out.longstretch1norm = 0;
    out.meanstretch1 = 0;
    out.meanstretch1norm = 0;
    out.stdstretch1 = NaN;
else
    out.longstretch1 = max(stretch1); % longest consecutive stretch of ones
    out.longstretch1norm = max(stretch1)/N; % longest consecutive stretch of ones as proportion of the time-series length
    out.meanstretch1 = mean(stretch1);
    out.meanstretch1norm = mean(stretch1)/N;
    out.stdstretch1 = std(stretch1);
    out.stdstretch1norm = std(stretch1)/N;
end

out.meanstretchdiff = (out.meanstretch1 - out.meanstretch0)/N;
out.stdstretchdiff = (out.stdstretch1 - out.stdstretch0)/N;

out.diff21stretch1 = mean(stretch1 == 2) - mean(stretch1 == 1);
out.diff21stretch0 = mean(stretch0 == 2) - mean(stretch0 == 1);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function xList=RKAdaptiveAtTimes(xInit,theTimes,f,initStepSize,order,solutionChoice,RelTol,AbsTol,maxSteps)
%%RKADAPTIVEATTIMES Perform multiple steps of Runge-Kutta propagation
%                   using an adaptive step size. Runge-Kutta methods are
%                   derivative-free techniques for solving ordinary
%                   differential equations. That is, integrating
%                   dx/dt=f(x,t) given initial conditions
%                   (xStart,tSpan(1)). More information on available
%                   algorithms for the steps is given in the comments to
%                   the function RungeKStep.
%
%INPUTS: xInit The initial value of the state (scalar or vector) over which
%              integration is being performed.
%     theTimes The times at which state estimates are desired. theTimes(1)
%              is the time of xInit.
%            f f(xVal,curT) returns the derivative of xVal taken at time
%              curT.
% initStepSize An optional initial step size (in t) to use for the
%              integration. If omitted or an empty matrix is passed, an
%              ad-hoc method is used to find an initial step size.
%        order The order of the Runge-Kutta method. If this parameter is
%              omitted, then the default order of 4 is used. Order can
%              range from 1 to 7.
%  solutionChoice When multiple formulae are implemented, this selects
%              which one to use. Otherwise, this parameter is not used.
%       RelTol The maximum relative error tolerance allowed, a positive
%              scalar. If omitted or an empty matrix is passed, the default
%              value of 1e-3 is used.
%       AbsTol The absolute error tolerance allowed, a positive scalar, of
%              the same for all components of x, or a positive NX1 vector.
%              If omitted or an empty matrix is passed, the default value
%              of 1e-6 is used.
%     maxSteps The maximum allowable number of steps to perform the
%              integration. If omitted, the default of 1024 is used.
%
%OUTPUTS: xList The state at the times given in times xList(:,1) is the
%               same as xInit.
%
%A detailed description of the adaptive step size algorithm can be found in
%the comments of RKAdaptiveOverRange.
%
%May 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<9||isempty(maxSteps))
    maxSteps=1024;
end

if(nargin<8||isempty(AbsTol))
    AbsTol=1e-6;
end

if(nargin<7||isempty(RelTol))
    RelTol=1e-3;
end
if(nargin<6||isempty(solutionChoice))
    solutionChoice=0;
end

if(nargin<5||isempty(order))
    order=5;
end

if(nargin<4||isempty(initStepSize))
    initStepSize=[];
end

xDim=length(xInit);
numTimes=length(theTimes);

xList=zeros(xDim,numTimes);
xList(:,1)=xInit;
for curTime=2:numTimes
    tSpan=theTimes(curTime-1:curTime);
    [xVals,~,~,~,~,initStepSize]=RKAdaptiveOverRange(xList(:,curTime-1),tSpan,f,initStepSize,0,order,solutionChoice,RelTol,AbsTol,maxSteps);
    xList(:,curTime)=xVals(:,end);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef ButtonBox < uix.Box
    %uix.ButtonBox  Button box base class
    %
    %  uix.ButtonBox is a base class for containers that lay out buttons.

    %  Copyright 2009-2015 The MathWorks, Inc.
    %  $Revision: 1165 $ $Date: 2015-12-06 03:09:17 -0500 (Sun, 06 Dec 2015) $

    properties( Access = public, Dependent, AbortSet )
        ButtonSize % button size, in pixels
        HorizontalAlignment % horizontal alignment [left|center|right]
        VerticalAlignment % vertical alignment [top|middle|bottom]
    end

    properties( Access = protected )
        ButtonSize_ = [60 20] % backing for ButtonSize
        HorizontalAlignment_ = 'center' % backing for HorizontalAlignment
        VerticalAlignment_ = 'middle' % backing for VerticalAlignment
    end

    methods

        function value = get.ButtonSize( obj )

            value = obj.ButtonSize_;

        end % get.ButtonSize

        function set.ButtonSize( obj, value )

            % Check
            assert( isa( value, 'double' ), 'uix:InvalidPropertyValue', ...
                'Property ''ButtonSize'' must be of type double.' )
            assert( isequal( size( value ), [1 2] ), ...
                'uix:InvalidPropertyValue', ...
                'Size of property ''ButtonSize'' must by 1-by-2.' )
            assert( all( isreal( value ) ) && ~any( isinf( value ) ) && ...
                ~any( isnan( value ) ) && ~any( value <= 0 ), ...
                'uix:InvalidPropertyValue', ...
                'Elements of property ''ButtonSize'' must be real, finite and positive.' )

            % Set
            obj.ButtonSize_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.ButtonSize

        function value = get.HorizontalAlignment( obj )

            value = obj.HorizontalAlignment_;

        end % get.HorizontalAlignment

        function set.HorizontalAlignment( obj, value )

            % Check
            assert( ischar( value ), 'uix:InvalidPropertyValue', ...
                'Property ''HorizontalAlignment'' must be a string.' )
            assert( any( strcmp( value, {'left';'center';'right'} ) ), ...
                'Property ''HorizontalAlignment'' must be ''left'', ''center'' or ''right''.' )

            % Set
            obj.HorizontalAlignment_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.HorizontalAlignment

        function value = get.VerticalAlignment( obj )

            value = obj.VerticalAlignment_;

        end % get.VerticalAlignment

        function set.VerticalAlignment( obj, value )

            % Check
            assert( ischar( value ), 'uix:InvalidPropertyValue', ...
                'Property ''VerticalAlignment'' must be a string.' )
            assert( any( strcmp( value, {'top';'middle';'bottom'} ) ), ...
                'Property ''VerticalAlignment'' must be ''top'', ''middle'' or ''bottom''.' )

            % Set
            obj.VerticalAlignment_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.VerticalAlignment

    end % accessors

end % classdef

author: goGPS-Project


text:
function [f] = spm_mc_fxa_4(x,v,a,P)
% equations of motion for the mountain car problem
% problem
% FORMAT [f] = spm_mc_fxa_4(x,v,a,P)
%
% x   - hidden states
% v   - exogenous inputs
% a   - action
% P   - parameters for mountain car
%
% returns f = dx/dt
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mc_fxa_4.m 3333 2009-08-25 16:12:44Z karl $

% physical flow
%--------------------------------------------------------------------------
dx  = spm_fx_mountaincar([x.x; x.v],v,a,P)/2;
f.x = dx(1);
f.v = dx(2);

% physiological flow
%--------------------------------------------------------------------------
A   = exp(-(x.x - 1).^2*32);
f.p = A - x.p/32;
f.d = (x.p - x.d)/64;

author: spm


text:
 function X = nufft2_bilinear(omega, x, K1, K2, n_shift)
%function X = nufft2_bilinear(omega, x, K1, K2, n_shift)
%	in:
%		omega	[M 2]		frequencies in radians
%		x	[N1 N2]		image dimensions
%		J			# of neighbors used
%		K1,K2			FFT sizes (should be > N1,N2)
%		n_shift [2]		n = 0-n_shift to N-1-n_shift
%	out:
%		X	[M 1]	DTFT2 at omega, approximated by bilinear interp
%
%	Like fft(), this expects the signals to be x(0,0), ...
%	Use n_shift = [N1/2 N2/2] for x(-N1/2,-N2/2), ...
%
% Copyright 2001-9-20, Jeff Fessler, University of Michigan

warn('This function is obsolete.  Use nufft_init with the "linear" option')

% if no arguments, then run a simple test
if nargin < 1
	N1 = 4;
	N2 = 8;
	n_shift = [2.7 3.3];
	x = [[1:N1]'*ones(1,3), ones(N1,N2-3)]; % test signal
%	x = zeros(N1,N2); x(1,1) = 1;
	if 0	% test with uniform frequency locations
		o1 = 2 * pi * [0:(N1-1)]' / N1;
		o2 = 2 * pi * [0:(N2-1)]' / N2;
		[o1, o2] = ndgrid(o1, o2);
		omega = [o1(:) o2(:)];
	else	% nonuniform frequencies
		o1 = [0 7.2 2.6 3.3]';
		o2 = [0 4.2 -1 5.5]';
		omega = [o1(:) o2(:)];
	end
	Xd = dtft2(x, omega, n_shift);
	Xb = nufft2_bilinear(omega, x, 2*N1, 2*N2, n_shift);
%	Xb = reshape(Xb, N1, N2)
%	disp([Xd Xb Xb-Xd])
	help(mfilename)
	disp(sprintf('max %% difference = %g', max_percent_diff(Xd,Xb)))
	return
end

if ~isvar('n_shift') || isempty(n_shift), n_shift = [0 0]; end
if ~isvar('useloop') || isempty(useloop), useloop = 0; end

M = size(omega,1);
[N1 N2] = size(x);

koff1 = bilinear_offset(omega(:,1), K1);	% [M 1] nearest
koff2 = bilinear_offset(omega(:,2), K2);	% [M 1] nearest

% indices into oversampled FFT components
k1 = mod(outer_sum(koff1, [0 1]), K1) + 1;	% [M 2] {1,...,K1}
k2 = mod(outer_sum(koff2, [0 1]), K2) + 1;	% [M 2] {1,...,K2}

% 1D interpolation coefficient vectors
[u1l u1r] = interp_coef(omega(:,1), koff1, K1, N1);	% [M 1]
[u2l u2r] = interp_coef(omega(:,2), koff2, K2, N2);	% [M 1]

% precorrect for bilinear interpolation filtering
if 0
	warning 'applying precorrection'
	[i1, i2] = ndgrid([0:(N1-1)]-n_shift(1), [0:(N2-1)]-n_shift(2));
	tmp = (nufft_sinc(i1 / K1) .* nufft_sinc(i2 / K2)).^2;
	printm('sinc correction min=%g', min(tmp(:)))
	x = x ./ tmp;
end

% FFT and bilinear interpolation
% can't use interp2 here because we want periodic end conditions!
% oh, maybe i could with a little padding of ends...
Xk = fft2(x, K1, K2);
kk11 = k1(:,1) + (k2(:,1) - 1) * K1;
kk21 = k1(:,2) + (k2(:,1) - 1) * K1;
kk12 = k1(:,1) + (k2(:,2) - 1) * K1;
kk22 = k1(:,2) + (k2(:,2) - 1) * K1;
t1 = u1l .* Xk(kk11) + u1r .* Xk(kk21);
t2 = u1l .* Xk(kk12) + u1r .* Xk(kk22);
X = u2l .* t1 + u2r .* t2;

% apply phase shift
phase = exp(i * (omega * n_shift(:)));	% [M 1]
X = X .* phase;

% index from 0 to K-1 (will modulo later)
function koff = bilinear_offset(om, K)
gam = 2*pi/K;
koff = floor(om / gam);

% make 1D interpolation coefficient vectors
function [ul, ur] = interp_coef(om, koff, K, N)
gam = 2*pi/K;
ul = 1 - (om/gam - koff);
ur = 1 - ul;
ul = ul .* kphase(om - koff * gam, N);
ur = ur .* kphase(om - (koff+1) * gam, N);

% natural phase function. trick: force it to be 2pi periodic
function ph = kphase(om, N)
ph = exp(-i * mod0(om,2*pi) * (N-1)/2);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [hIm, ww] = L1SR(lIm, zooming, patch_size, overlap, Dh, Dl, lambda, regres)
% Use sparse representation as the prior for image super-resolution
% Usage
%       [hIm] = L1SR(lIm, zooming, patch_size, overlap, Dh, Dl, lambda)
%
% Inputs
%   -lIm:           low resolution input image, single channel, e.g.
%   illuminance
%   -zooming:       zooming factor, e.g. 3
%   -patch_size:    patch size for the low resolution image
%   -overlap:       overlap among patches, e.g. 1
%   -Dh:            dictionary for the high resolution patches
%   -Dl:            dictionary for the low resolution patches
%   -regres:       'L1' use the sparse representation directly to high
%                   resolution dictionary;
%                   'L2' use the supports found by sparse representation
%                   and apply least square regression coefficients to high
%                   resolution dictionary.
% Ouputs
%   -hIm:           the recovered image, single channel
%
% Written by Jianchao Yang @ IFP UIUC
% April, 2009
% Webpage: http://www.ifp.illinois.edu/~jyang29/
% For any questions, please email me by jyang29@uiuc.edu
%
% Reference
% Jianchao Yang, John Wright, Thomas Huang and Yi Ma. Image superresolution
% as sparse representation of raw image patches. IEEE Computer Society
% Conference on Computer Vision and Pattern Recognition (CVPR), 2008.
%

[lhg, lwd] = size(lIm);
hhg = lhg*zooming;
hwd = lwd*zooming;

mIm = imresize(lIm, 2,'bicubic');
[mhg, mwd] = size(mIm);
hpatch_size = patch_size*zooming;
mpatch_size = patch_size*2;

% extract gradient feature from lIm
hf1 = [-1,0,1];
vf1 = [-1,0,1]';
hf2 = [1,0,-2,0,1];
vf2 = [1,0,-2,0,1]';

lImG11 = conv2(mIm,hf1,'same');
lImG12 = conv2(mIm,vf1,'same');
lImG21 = conv2(mIm,hf2,'same');
lImG22 = conv2(mIm,vf2,'same');

lImfea(:,:,1) = lImG11;
lImfea(:,:,2) = lImG12;
lImfea(:,:,3) = lImG21;
lImfea(:,:,4) = lImG22;

lgridx = 2:patch_size-overlap:lwd-patch_size;
lgridx = [lgridx, lwd-patch_size];
lgridy = 2:patch_size-overlap:lhg-patch_size;
lgridy = [lgridy, lhg-patch_size];

mgridx = (lgridx - 1)*2 + 1;
mgridy = (lgridy - 1)*2 + 1;

% using linear programming to find sparse solution
bhIm = imresize(lIm, zooming, 'bicubic');
hIm = zeros([hhg, hwd]);
nrml_mat = zeros([hhg, hwd]);

hgridx = (lgridx-1)*zooming + 1;
hgridy = (lgridy-1)*zooming + 1;

disp('Processing the patches sequentially...');
count = 0;

%ProjM = inv(Dl'*Dl+0.001*eye(size(Dl,2)))*Dl';

% loop to recover each patch
for xx = 1:length(mgridx),
    for yy = 1:length(mgridy),

        mcolx = mgridx(xx);
        mrowy = mgridy(yy);

        count = count + 1;
%         if ~mod(count, 10000),
%             fprintf('.\n');
%         else
%             fprintf('.');
%         end;
        if ~mod(count, 1000),
            fprintf('%g/%g\n',count,length(mgridx)*length(mgridy));
        end
        mpatch = mIm(mrowy:mrowy+mpatch_size-1, mcolx:mcolx+mpatch_size-1);
        mmean = mean(mpatch(:));

        mpatchfea = lImfea(mrowy:mrowy+mpatch_size-1, mcolx:mcolx+mpatch_size-1, :);
        mpatchfea = mpatchfea(:);

        mnorm = sqrt(sum(mpatchfea.^2));

        if mnorm > 1,
            y = mpatchfea./mnorm;
        else
            y = mpatchfea;
        end;
        %w = ProjM*y;
        w = SolveLasso(Dl, y, size(Dl, 2), 'nnlasso', [], lambda);
        %w = feature_sign(Dl, y, lambda*2);

        if isempty(w),
            w = zeros(size(Dl, 2), 1);
        end;
        switch regres,
            case 'L1'
                if mnorm > 1,
                    hpatch = Dh*w*mnorm;
                else
                    hpatch = Dh*w;
                end;
            case 'L2'
                idx = find(w);
                lsups = Dl(:, idx);
                hsups = Dh(:, idx);
                w = inv(lsups'*lsups)*lsups'*mpatchfea;
                hpatch = hsups*w;
            otherwise
                error('Unknown fitting!');
        end;

        hpatch = reshape(hpatch, [hpatch_size, hpatch_size]);
        hpatch = hpatch + mmean;

        hcolx = hgridx(xx);
        hrowy = hgridy(yy);

        hIm(hrowy:hrowy+hpatch_size-1, hcolx:hcolx+hpatch_size-1)...
            = hIm(hrowy:hrowy+hpatch_size-1, hcolx:hcolx+hpatch_size-1) + hpatch;
        nrml_mat(hrowy:hrowy+hpatch_size-1, hcolx:hcolx+hpatch_size-1)...
            = nrml_mat(hrowy:hrowy+hpatch_size-1, hcolx:hcolx+hpatch_size-1) + 1;
    end;
end;

fprintf('done!\n');

% fill the empty
hIm(1:3, :) = bhIm(1:3, :);
hIm(:, 1:3) = bhIm(:, 1:3);

hIm(end-2:end, :) = bhIm(end-2:end, :);
hIm(:, end-2:end) = bhIm(:, end-2:end);

nrml_mat(nrml_mat < 1) = 1;
hIm = hIm./nrml_mat;
% hIm = uint8(hIm);

author: thomas-koehler

repo: SupER


text:
function [t,x] = unrz(bits, bitrate)
% UNRZ Encode bit string using unipolar NRZ code.
%   [T, X] = UNRZ(BITS, BITRATE) encodes BITS array using unipolar NRZ
%   code with given BITRATE. Outputs are time T and encoded signal
%   values X.

% Copyright (c) 2013 Yuriy Skalko <yuriy.skalko@gmail.com>

T = length(bits)/bitrate; % full time of bit sequence
n = 200;
N = n*length(bits);
dt = T/N;
t = 0:dt:T;
x = zeros(1,length(t)); % output signal

for i = 0:length(bits)-1
  if bits(i+1) == 1
    x(i*n+1:(i+1)*n) = 1;
  else
    x(i*n+1:(i+1)*n) = 0;
  end
end

author: Sable

repo: mcbench-benchmarks


text:
function [dStates, contactForces] = dynamics_singleStanceOne(States, Actuators, Parameters)
% function [dStates, Actuators] = dynamics_singleStanceOne(States, Actuators, Parameters)
%
% Computer Generated File -- DO NOT EDIT
%
% This function was created by the function Write_ContinuousDynamics()
% 10-Dec-2013 19:40:47
%
% Dymanics Model: retractable double pendulum biped
% Motion Phase: Single Stance One
%
% Matthew Kelly
% Cornell University
%

x = States(1,:); % (m) Foot One horizontal position
y = States(2,:); % (m) Foot One vertical position
th1 = States(3,:); % (rad) Leg One absolute angle
th2 = States(4,:); % (rad) Leg Two absolute angle
L1 = States(5,:); % (m) Leg One length
L2 = States(6,:); % (m) Leg Two length
dx = States(7,:); % (m/s) Foot One horizontal velocity
dy = States(8,:); % (m/s) Foot One vertical velocity
dth1 = States(9,:); % (rad/s) Leg One absolute angular rate
dth2 = States(10,:); % (rad/s) Leg Two absolute angular rate
dL1 = States(11,:); % (m/s) Leg One extension rate
dL2 = States(12,:); % (m/s) Leg Two extensioin rate

F1 = Actuators(1,:); % (N) Compresive axial force in Leg One
F2 = Actuators(2,:); % (N) Compresive axial force in Leg Two
T1 = Actuators(3,:); % (Nm) External torque applied to Leg One
T2 = 0; %Foot Two not in contact with ground!
Thip = Actuators(5,:); % (Nm) Torque acting on Leg Two from Leg One

m1 = Parameters.m1; % (kg) Foot One mass
m2 = Parameters.m2; % (kg) Foot Two mass
M = Parameters.M; % (kg) Hip mass
g = Parameters.g; % (m/s^2) Gravity

% Constraints for this phase:
ddx = 0; %(m) Foot One horizontal position
ddy = 0; %(m) Foot One vertical position
H2 = 0; %(N) Foot Two, horizontal contact force
V2 = 0; %(N) Foot Two, vertical contact force

dStates = zeros(size(States));
dStates(1:6,:) = States((1+6):(6+6),:);
dStates(7,:) = zeros(1,size(States,2));
dStates(8,:) = zeros(1,size(States,2));
dStates(9,:) = -(L2.*Thip - L2.*T1 + L1.*T2.*cos(th1 - th2) + L1.*Thip.*cos(th1 - th2) - F2.*L1.*L2.*sin(th1 - th2) + 2.*L1.*L2.*M.*dL1.*dth1 + L1.*L2.*M.*ddy.*cos(th1) + L1.*L2.*M.*g.*cos(th1) - L1.*L2.*M.*ddx.*sin(th1))./(L1.^2.*L2.*M);
dStates(10,:) = (L1.*M.*T2 + L1.*M.*Thip - L2.*T1.*m2.*cos(th1 - th2) + L2.*Thip.*m2.*cos(th1 - th2) + L1.*T2.*m2.*cos(th1 - th2).^2 + L1.*Thip.*m2.*cos(th1 - th2).^2 + L1.*T2.*m2.*sin(th1 - th2).^2 + L1.*Thip.*m2.*sin(th1 - th2).^2 + L1.*L2.*M.*V2.*cos(th2) - H2.*L1.*L2.*M.*sin(th2) - F1.*L1.*L2.*m2.*sin(th1 - th2) - 2.*L1.*L2.*M.*dL2.*dth2.*m2 - L1.*L2.*M.*ddy.*m2.*cos(th2) - L1.*L2.*M.*g.*m2.*cos(th2) + L1.*L2.*M.*ddx.*m2.*sin(th2) - L1.*L2.*M.*ddx.*m2.*cos(th1 - th2).*sin(th1) + L1.*L2.*M.*ddx.*m2.*sin(th1 - th2).*cos(th1) + L1.*L2.*M.*ddy.*m2.*sin(th1 - th2).*sin(th1) + L1.*L2.*M.*g.*m2.*sin(th1 - th2).*sin(th1) + L1.*L2.*M.*ddy.*m2.*cos(th1 - th2).*cos(th1) + L1.*L2.*M.*g.*m2.*cos(th1 - th2).*cos(th1))./(L1.*L2.^2.*M.*m2);
dStates(11,:) = -(T2.*sin(th1 - th2) + Thip.*sin(th1 - th2) - F1.*L2 + F2.*L2.*cos(th1 - th2) - L1.*L2.*M.*dth1.^2 + L2.*M.*ddx.*cos(th1) + L2.*M.*ddy.*sin(th1) + L2.*M.*g.*sin(th1))./(L2.*M);
dStates(12,:) = (T1.*m2.*sin(th1 - th2) - Thip.*m2.*sin(th1 - th2) + F2.*L1.*M + H2.*L1.*M.*cos(th2) + L1.*M.*V2.*sin(th2) - F1.*L1.*m2.*cos(th1 - th2) + F2.*L1.*m2.*cos(th1 - th2).^2 + F2.*L1.*m2.*sin(th1 - th2).^2 + L1.*L2.*M.*dth2.^2.*m2 - L1.*M.*ddx.*m2.*cos(th2) - L1.*M.*ddy.*m2.*sin(th2) - L1.*M.*g.*m2.*sin(th2) + L1.*M.*ddx.*m2.*cos(th1 - th2).*cos(th1) + L1.*M.*ddy.*m2.*cos(th1 - th2).*sin(th1) - L1.*M.*ddy.*m2.*sin(th1 - th2).*cos(th1) + L1.*M.*g.*m2.*cos(th1 - th2).*sin(th1) - L1.*M.*g.*m2.*sin(th1 - th2).*cos(th1) + L1.*M.*ddx.*m2.*sin(th1 - th2).*sin(th1))./(L1.*M.*m2);

% contactForces(1,:) == H1 == (N) Foot One, horizontal contact force
% contactForces(2,:) == V1 == (N) Foot One, vertical contact force
% contactForces(3,:) == H2 == (N) Foot Two, horizontal contact force
% contactForces(4,:) == V2 == (N) Foot Two, vertical contact force
contactForces = zeros(4,size(States,2));
contactForces(1,:) = (Thip.*sin(th1) - T1.*sin(th1) + L1.*ddx.*m1 + F1.*L1.*cos(th1))./L1;
contactForces(2,:) = (T1.*cos(th1) - Thip.*cos(th1) + L1.*ddy.*m1 + L1.*g.*m1 + F1.*L1.*sin(th1))./L1;
contactForces(3,:) = zeros(1,size(States,2));
contactForces(4,:) = zeros(1,size(States,2));

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [ pn, dist ] = circle_imp_point_near_2d ( r, center, p )

%*****************************************************************************80
%
%% CIRCLE_IMP_POINT_NEAR_2D: nearest ( implicit circle, point ) in 2D.
%
%  Discussion:
%
%    This routine finds the distance from a point to an implicitly
%    defined circle, and returns the point on the circle that is
%    nearest to the given point.
%
%    If the given point is the center of the circle, than any point
%    on the circle is "the" nearest.
%
%    An implicit circle in 2D satisfies the equation:
%
%      ( X - CENTER(1) )**2 + ( Y - CENTER(2) )**2 = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2), the center of the circle.
%
%    Input, real P(2), the point to be checked.
%
%    Output, real PN(2), the nearest point on the circle.
%
%    Output, real DIST, the distance of the point to the circle.
%
  ndim = 2;

  if ( p(1:ndim) == center(1:ndim) )
    dist = r;
    pn(1:ndim) = center(1:ndim) + r / sqrt ( ndim );
    return
  end

  r2 = sqrt ( sum ( ( p(1:ndim) - center(1:ndim) ).^2 ) );

  dist = abs (  r2 - r );

  pn(1:ndim) = center(1:ndim) + r * ( p(1:ndim) - center(1:ndim) ) / r2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%{
 * Copyright (C) 2013-2020, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * 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 OWNER 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.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function cost = cost4Points(theta_x, theta_y, theta_z, T, X, Y, intrinsic)

%         function cost = cost4Points(app, xi, X, Y, intrinsic)
%             H = expm(hat(app, xi'));
%             R = H(1:3, 1:3);
%             T = H(1:3, 4);
    H = eye(4);
    H(1:3,1:3) = rotx(theta_x) * roty(theta_y) * rotz(theta_z);
    H(1:3,4) = T';
    P = intrinsic * [eye(3) zeros(3,1)] * H;
    cost = verifyCornerAccuracy(X, Y, P);
end

author: UMich-BipedLab


text:
function scn_spm_choose_hpfilter(spm_results_dir, varargin)
% Plots and choice of optimal high-pass filter from an SPM first-level
% model directory (with statistics and contrasts estimated.)
%
% :Usage:
% ::
%
%     scn_spm_choose_hpfilter(spm_results_dir, ['events_only'],['from_multireg',nr_regs)
%
% SPM5 compatible and SPM8.
%
% Called by: scn_spm_design_check.m
% For all regressors or events only: see scn_spm_get_events_of_interest.m
%
% ..
%    Tor Wager
%    August 2010
%    Lukas Van Oudenhove
%    March 2022: added 'from_multireg' option which was already implemented
%    in scn_spm_get_events_of_interest.m
% ..

if nargin < 1, spm_results_dir = pwd; end

spmfilename = fullfile(spm_results_dir, 'SPM.mat');
if ~exist(spmfilename, 'file')
    error('SPM.mat does not exist in %s\n', spm_results_dir);
end
load(spmfilename);

EVENTS_ONLY = false;
FROM_MULTIREG = false;
SORTBYVIF = false;

i=1;
while i<=numel(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case 'events_only'
                EVENTS_ONLY = true;
            case 'from_multireg' % added lukasvo76
                FROM_MULTIREG = true;
                i=i+1;
                nr_regs = varargin{i};
            otherwise
                error(['Unrecognized argument: ' varargin{i}])
        end
    else
        error(['Unrecognized argument: ' varargin{i}])
    end
    i=i+1;
end

% Gets events of interest: All regressors, or events only if 'events_only'
% is input as keyword, or first nr_regs if 'from_multireg' is additionally
% specified
if EVENTS_ONLY
    if FROM_MULTIREG % added lukasvo76
        wh_cols = scn_spm_get_events_of_interest(SPM, 'events_only','from_multireg',nr_regs);
    else
        wh_cols = scn_spm_get_events_of_interest(SPM, 'events_only');
    end
else
    wh_cols = scn_spm_get_events_of_interest(SPM);
end

% Prepare design columns of interest
X_interest = SPM.xX.X(:, wh_cols);
I = SPM.xX.X(:, SPM.xX.iB);
X_interest = X_interest - I * pinv(I) * X_interest;

k = size(X_interest, 2);

%% Start by plotting the regressors of interest

create_figure('regressors', 2, 2);
set(gcf,'WindowState','maximized'); % lukasvo76 edited to improve figure quality
delete(subplot(2, 2, 1:2));
axh = axes('Position', [.05 .55 .9 .4]);
set(axh, 'FontSize', 16);

nregs = length(wh_cols);

if k < 30 % edited lukasvo76
    plot_matrix_cols(X_interest, 'horiz');
else
    imagesc(X_interest); colormap gray; colorbar
end

%% Plot the regressors after HP filtering in red

spm_hplength = SPM.xX.K(1).HParam;

X_filtered = filter_X(X_interest, cat(1, SPM.xX.K(:).X0));

hold on;
if k < 30 % edited lukasvo76
    han = plot_matrix_cols(X_filtered, 'horiz');
    set(han, 'Color', 'r');
else
    imagesc(X_filtered); colormap gray; colorbar
end

title(['Task regressors: Black = before filtering, Red = after current ' num2str(spm_hplength) ' s HP filter']);
drawnow

%% Calculate the cumulative power at different frequencies
TR = SPM.xY.RT;

[cumpower, myfreq, design_var_loss, five_percent_HPlength] = ...
    cumulative_power(X_interest, TR);

currenthp = mean(cat(1, SPM.xX.K.HParam));
best_hplen = five_percent_HPlength;

% Plot it
subplot(2, 2, 3);

for i = 1:nregs
    plot(myfreq(:, i), cumpower(:, i), 'k');
end
xlabel('Frequency (Hz)') ; ylabel('Cumpow (regs)');
axis tight
title(sprintf('<5%% regressor var. lost at HPlen = %3.0f sec', five_percent_HPlength))

try
    plot_vertical_line(currenthp, 'r');
    plot_vertical_line(1./best_hplen, 'b');
catch
    disp('Bug! check me')
    keyboard
end

drawnow

%% If Contrasts are entered, then plot those too
X_contrasts = get_contrast_vectors(SPM);
% can handle t- or F-contrasts

subplot(2, 2, 4);
if isempty(X_contrasts)
    title('NO Contrasts Entered!');

else
    [cumpower, myfreq, design_var_loss, five_percent_HPlength_contrasts] = ...
        cumulative_power(X_contrasts, TR);

    best_hplen = five_percent_HPlength_contrasts;

    % Plot it
    for i = 1:size(X_contrasts, 2)
        plot(myfreq(:, i), cumpower(:, i), 'k');
    end
    xlabel('Frequency (Hz)') ; ylabel('Cumpow (contrasts)');
    axis tight
    title(sprintf('<5%% con-variance lost at HPlen = %3.0f sec', five_percent_HPlength_contrasts))

    try
        plot_vertical_line(currenthp, 'r');
        plot_vertical_line(1./best_hplen, 'b');
    catch
        disp('Bug! check me')
        keyboard
    end

end

drawnow

%% Now re-filter the regressors with the best length

[Rmtx, dummy, K] = use_spm_filter(TR, size(X_interest, 1), 'none', 'specify', best_hplen);

X_filtered = filter_X(X_interest, K);

axes(axh); hold on;

if k < 20
    han = plot_matrix_cols(X_filtered, 'horiz');
    set(han, 'Color', 'b');
else
    imagesc(X_filtered); colormap gray; colorbar
end

text(1, nregs + .7, sprintf('Blue = filtering with <=5%% variance loss at %3.0f sec', best_hplen), 'FontSize', 16);
drawnow

end

function X_filtered = filter_X(X_interest, K)

%px = pinv(K);           % pinv of the filtering matrix
Rmtx = eye(size(K, 1)) - K * pinv(K);  % Residual-inducing matrix, y' = Rmtx*y

for i = 1:size(X_interest, 2)

    y = X_interest(:, i);             % select regressor to filter
    %y = y - K * px * y;     % residuals after filtering
    y = Rmtx * y;

    X_filtered(:, i) = scale(y);

end

end

function [cumpower, myfreq, design_var_loss, five_percent_HPlength] = cumulative_power(X_interest, TR)

for i = 1:size(X_interest, 2)

    [myfft, freq] = fft_calc(X_interest(:, i), TR);

    cumpower(:, i) = cumsum(myfft);
    myfreq(:, i) = freq;

    wh = find( cumpower(:, i) < .05 );
    if isempty(wh)
        wh = 2; % first one after intercept
        design_var_loss(i) = cumpower(wh, i);
        fprintf('No cutoff with < .5%% power for reg/contrast %3.0f! Min loss is %3.2f%%\n', i, 100*design_var_loss(i))
    end

    wh = wh(end);  % highest index value with < 5% power loss

    design_var_loss(i) = cumpower(wh, i);
    five_percent_loss_freq(i) = myfreq(wh, i);

end

% HP filter length that preserves 5% design variance loss
% or less in all regressors
five_percent_HPlength = 1 ./ min(five_percent_loss_freq);

end

function X_contrasts = get_contrast_vectors(SPM)
% can handle t- or F-contrasts
X_contrasts = [];
if ~isfield(SPM, 'xCon') || length(SPM.xCon) == 0
    return
end

% Remove intercepts
X = SPM.xX.X;
I = SPM.xX.X(:, SPM.xX.iB);
X = X - I * pinv(I) * X;

for mycon = 1:length(SPM.xCon)

    d = X * SPM.xCon(mycon).c;
    X_contrasts = [X_contrasts d];

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function crf = tapas_physio_crf(t)
% cardiac response function, as described in
%
% Chang, Catie, and Gary H. Glover. �Effects of Model-based Physiological
% Noise Correction on Default Mode Network Anti-correlations and Correlations.�
% NeuroImage 47, no. 4 (October 1, 2009): 1448�1459. doi:10.1016/j.neuroimage.2009.05.012.
%
% following:
% Chang, C., Cunningham, J.P., Glover, G.H., 2009. Influence of heart rate
% on the BOLD signal: the cardiac response function. Neuroimage 44, 857�869.%
%
%   crf = tapas_physio_crf(t)
%
% IN
%   t       vector of timepoints
% OUT
%   crf     cardiac response function at sampled time points
%
% EXAMPLE
%   % just for visualization
%   t = 0:0.1:100;
%   crf = tapas_physio_crf(t);
%   figure;plot(t,crf);
%
%   See also tapas_physio_rrf

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

crf = 0.6*t.^2.7.*exp(-t/1.6) - 16/sqrt(2*pi*9).*exp(-1/2.*(t-12).^2/9);

author: translationalneuromodeling

repo: tapas


text:
function gX = whitehKernGradX(kern, X, X2)

% WHITEHKERNGRADX Gradient of WHITEH kernel with respect to input locations.
% FORMAT
% DESC computes the gradident of the whiteh noise
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO whitehKernParamInit, kernGradX, whitehKernDiagGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

if nargin<3
    X2 = X;
end

gX = zeros(size(X2, 1), size(X2, 2), size(X, 1));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [H, Z, S, psi] = sfactorization_wilson3x3(S,freq,Niterations,tol,cmbindx,fb,init,checkflag,stabilityfix)

% SFACTORIZATION_WILSON3X3 performs triplet-wise non-parametric spectral factorization on
% cross-spectra, based on Wilson's algorithm.
%
% Usage  : [H, Z, psi] = sfactorization_wilson(S,freq);
%
% Inputs : S (1-sided, 3D-spectral matrix in the form of Channel x Channel x frequency)
%        : freq (a vector of frequencies) at which S is given.
%
% Outputs: H (transfer function)
%        : Z (noise covariance)
%        : S (cross-spectral density 1-sided)
%        : psi (left spectral factor)
%
% This function is an implemention of Wilson's algorithm (Eq. 3.1)
% for spectral matrix factorization.
%
% Ref: G.T. Wilson,"The Factorization of Matricial Spectral Densities,"
% SIAM J. Appl. Math.23,420-426(1972).
% Written by M. Dhamala & G. Rangarajan, UF, Aug 3-4, 2006.
% Email addresses: mdhamala@bme.ufl.edu, rangaraj@math.iisc.ernet.in

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

if nargin<9, stabilityfix = false; end
if nargin<8, checkflag = true;   end
if nargin<7, init      = 'chol'; end
if nargin<6, fb        = 'none'; end
if nargin<5
  ft_error('FieldTrip:connectivity:sfactorization_wilson3x3', 'when requesting multiple triplet-wise spectral decomposition, ''cmbindx'' needs to be specified');
end
if nargin<4, tol        = 1e-8;   end
if nargin<3, Niterations = 1000;  end

dfreq = round(diff(freq)*1e5)./1e5; % allow for some numeric issues
if ~all(dfreq==dfreq(1))
  ft_error('FieldTrip:connectivity:sfactorization_wilson3x3', 'frequency axis is not evenly spaced');
end

if freq(1)~=0
  ft_warning('FieldTrip:connectivity:sfactorization_wilson3x3', 'when performing non-parametric spectral factorization, the frequency axis should ideally start at 0, zero padding the spectral density');
  dfreq = mean(dfreq);
  npad  = freq(1)./dfreq;

  % update the freq axis and keep track of the frequency bins that are
  % expected in the output
  selfreq  = (1:numel(freq)) + npad;
  freq     = [(0:(npad-1))./dfreq freq];
  S        = cat(3, zeros(size(S,1), size(S,1), npad), S);
else
  selfreq  = 1:numel(freq);
end

% ensure input S is double (mex-files don't work with single)
S = double(S);

% check whether the last frequency bin is strictly real-valued.
% if that's the case, then it is assumed to be the Nyquist frequency
% and the two-sided spectral density will have an even number of
% frequency bins. if not, in order to preserve hermitian symmetry,
% the number of frequency bins needs to be odd.
Send = S(:,:,end);
N    = numel(freq);
m    = size(cmbindx,1);
if all(imag(Send(:))<abs(trace(Send)./size(Send,1)*1e-9))
  hasnyq = true;
  N2     = 2*(N-1);
else
  hasnyq = false;
  N2     = 2*(N-1)+1;
end

% preallocate memory for the identity matrix
I      = repmat(eye(3),[1 1 m N2]); % Defining 3 x 3 identity matrix

% %Step 1: Forming 2-sided spectral densities for ifft routine in matlab
% Sarr   = zeros(3,3,m,N2) + 1i.*zeros(3,3,m,N2);
% for c = 1:m
%   Stmp  = S(cmbindx(c,:),cmbindx(c,:),:);
%
%   % the input cross-spectral density is assumed to be weighted with a
%   % factor of 2 in all non-DC and Nyquist bins, therefore weight the
%   % DC-bin with a factor of 2 to get a correct two-sided representation
%   Sarr(:,:,c,1) = Stmp(:,:,1).*2;
%
%   for f_ind = 2:N
%     Sarr(:,:,c,       f_ind) = Stmp(:,:,f_ind);
%     Sarr(:,:,c,(N2+2)-f_ind) = Stmp(:,:,f_ind).';
%   end
% end
% Sarr2 = Sarr;

% preallocate memory for the 2-sided spectral density
Sarr = zeros(3,3,N2,m) + 1i.*zeros(3,3,N2,m);
for c = 1:m
  Sarr(:,:,1:N,c) = S(cmbindx(c,:),cmbindx(c,:),:);
end
if hasnyq
  N1 = N;
else
  N1 = N + 1; % the highest frequency needs to be represented twice, for symmetry purposes
end
Sarr(:,:,N1:N2,:) = flip(Sarr(:,:,2:N,:),3);
Sarr(2,1,N1:N2,:) = conj(Sarr(2,1,N1:N2,:));
Sarr(1,2,N1:N2,:) = conj(Sarr(1,2,N1:N2,:));
Sarr(3,1,N1:N2,:) = conj(Sarr(3,1,N1:N2,:));
Sarr(1,3,N1:N2,:) = conj(Sarr(1,3,N1:N2,:));
Sarr(3,2,N1:N2,:) = conj(Sarr(3,2,N1:N2,:));
Sarr(2,3,N1:N2,:) = conj(Sarr(2,3,N1:N2,:));
Sarr              = permute(Sarr, [1 2 4 3]);
Sarr(:,:,:,1)     = Sarr(:,:,:,1).*2; % weight the DC-bin

% the input cross-spectral density is assumed to be weighted with a
% factor of 2 in all non-DC and Nyquist bins, therefore weight the
% Nyquist bin with a factor of 2 to get a correct two-sided representation
if hasnyq
  Sarr(:,:,:,N) = Sarr(:,:,:,N).*2;
end

%Step 2: Computing covariance matrices
gam = real(reshape(ifft(reshape(Sarr, [9*m N2]), [], 2),[3 3 m N2]));

%Step 3: Initializing for iterations
gam0 = gam(:,:,:,1);

h    = complex(zeros(size(gam0)));
for k = 1:m
  switch init
    case 'chol'
      [tmp, dum] = chol(gam0(:,:,k));
      if dum
        ft_warning('initialization with ''chol'' for iterations did not work well, using arbitrary starting condition');
        tmp = rand(3,3); %arbitrary initial condition
        tmp = triu(tmp);
      end
    case 'rand'
      tmp = rand(3,3); %arbitrary initial condition
      tmp = triu(tmp);
    otherwise
      ft_error('initialization method should be eithe ''chol'' or ''rand''');
  end
  h(:,:,k) = tmp;

  %h(:,:,k) = chol(gam0(:,:,k));
end
psi  = repmat(h, [1 1 1 N2]);

%Step 4: Iterating to get spectral factors
ft_progress('init', fb, 'computing spectral factorization');
for iter = 1:Niterations
  ft_progress(iter./Niterations, 'computing iteration %d/%d\n', iter, Niterations);
  invpsi = inv3x3(psi);
  g      = sandwich3x3(invpsi, Sarr) + I;
  gp     = PlusOperator3x3(g,m,N,stabilityfix); %gp constitutes positive and half of zero lags

  psi_old = psi;
  psi     = mtimes3x3(psi, gp);

  if checkflag
    psierr  = abs(psi-psi_old)./abs(psi);
    psierrf = mean(psierr(:));
    if(psierrf<tol)
      fprintf('reaching convergence at iteration %d\n',iter);
      break;
    end % checking convergence
  end
end
ft_progress('close');

%Step 5: Getting covariance matrix from spectral factors
gamtmp = reshape(real(ifft(transpose(reshape(psi, [9*m N2]))))', [3 3 m N2]);

%Step 6: Getting noise covariance & transfer function (see Example pp. 424)
%
A0    = gamtmp(:,:,:,1);
A0inv = inv3x3(A0);

Z = zeros(3,3,m);
for k = 1:m
  %Z     = A0*A0.'*fs; %Noise covariance matrix
  Z(:,:,k) = A0(:,:,k)*A0(:,:,k).'; %Noise covariance matrix not multiplied by sampling frequency
  %FIXME check this; at least not multiplying it removes the need to correct later on
  %this also makes it more equivalent to the noisecov estimated by biosig's mvar-function
end

% H = complex(zeros(3,3,m,N));
% S = complex(zeros(3,3,m,N));
% for k = 1:N
%   for kk = 1:m
%     H(:,:,kk,k) = psi(:,:,kk,k)*A0inv(:,:,kk);  % Transfer function
%     S(:,:,kk,k) = psi(:,:,kk,k)*psi(:,:,kk,k)'; % Cross-spectral density
%   end
% end
H = mtimes3x3(psi,A0inv(:,:,:,ones(1,size(psi,4))));
S = mtimes3x3(psi,ctranspose3x3(psi));

siz = [size(H) 1 1];
H   = reshape(H, [9*siz(3) siz(4:end)]);
siz = [size(S) 1 1];
S   = reshape(S, [9*siz(3) siz(4:end)]);
siz = [size(Z) 1 1];
Z   = reshape(Z, [9*siz(3) siz(4:end)]);
siz = [size(psi) 1 1];
psi = reshape(psi, [9*siz(3) siz(4:end)]);

%if numel(selfreq)~=numel(freq)
  % return only the frequency bins that were in the input
  H   =   H(:,selfreq,:,:);
  S   =   S(:,selfreq,:,:);
  psi = psi(:,selfreq,:,:);
%end

%---------------------------------------------------------------------
function gp = PlusOperator3x3(g,ncmb,nfreq, stabilityfix)

% This function is for [ ]+operation:
% to take the positive lags & half of the zero lag and reconstitute
% M. Dhamala, UF, August 2006

g   = transpose(reshape(g, 9*ncmb, []));
gam = ifft(g);

% taking only the positive lags and half of the zero lag
gamp  = gam;
beta0 = 0.5*gam(1,:);

%for k = 1:ncmb
%  gamp(1,(k-1)*4+1:k*4) = reshape(triu(reshape(beta0(1,(k-1)*4+1:k*4),[2 2])),[1 4]);
%end
beta0(2:9:9*ncmb)   = 0;
beta0(3:9:9*ncmb)   = 0;
beta0(6:9:9*ncmb)   = 0;
gamp(1,:)           = beta0;
gamp(nfreq+1:end,:) = 0;

% smooth with a window, only for the long latency boundary: this is a
% stabilityfix proposed by Martin Vinck
if stabilityfix
  w = tukeywin(nfreq*2, 0.5);
  gamp(1:nfreq,:) = gamp(1:nfreq,:).*repmat(w(nfreq+1:end),[1 nchan^2]);
else
  % nothing to be done here
end

% reconstituting
gp = fft(gamp);
gp = reshape(transpose(gp), [3 3 ncmb numel(gp)/(9*ncmb)]);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% RING_POLAR_MAP_DER: Jacobian of the parameterization given in ring_polar_map.

function varargout = ring_polar_map_der (pts)

  if (iscell (pts))
    u = reshape (repmat (pts{1}(:), 1, numel (pts{2})), 1, []);
    v = reshape (repmat (pts{2}(:)', numel (pts{1}), 1), 1, []);
  else
    u = pts(1,:); v = pts(2,:);
  end

  jac = zeros (2, 2, numel (u));
  jac(1,1,:) =  cos (pi.*v/2);
  jac(2,1,:) =  sin (pi.*v/2);
  jac(1,2,:) = -pi*(u+1) .* sin (pi.*v/2)/2;
  jac(2,2,:) =  pi*(u+1) .* cos (pi.*v/2)/2;

  if (nargout == 1)
    varargout{1} = jac;
  elseif (nargout == 2)
    F = ring_polar_map (pts);
    varargout{1} = F;
    varargout{2} = jac;
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
% GRADIENT CALLBACK function.
% This function is used as a callback function for fminunc and it aggregates
% cost and gradient values.
function [cost, gradients] = gradient_callback(X, y, theta, lambda)
    % X - training set.
    % y - training output values.
    % theta - model parameters.
    % lambda - regularization parameter.

    % Calculate cost function.
    cost = cost_function(X, y, theta, lambda);

    % Do one gradient step.
    gradients = gradient_step(X, y, theta, lambda);
end

author: trekhleb

repo: machine-learning-octave


text:
function FlyBowlJAABADetect(expdir,varargin)

version = '0.2';

[analysis_protocol,settingsdir,datalocparamsfilestr,...
  forcecompute,DEBUG] = ...
  myparse(varargin,'analysis_protocol','current',...
  'settingsdir','/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings',...
  'datalocparamsfilestr','dataloc_params.txt',...
  'forcecompute',false,...
  'debug',false);

if isunix && settingsdir(1) ~= '/',
  warning('settingsdir path must be the global path');
end

if ischar(forcecompute),
  forcecompute = str2double(forcecompute);
end
if ischar(DEBUG),
  DEBUG = str2double(DEBUG);
end

datalocparamsfile = fullfile(settingsdir,analysis_protocol,datalocparamsfilestr);
dataloc_params = ReadParams(datalocparamsfile);
paramsfile = fullfile(settingsdir,analysis_protocol,dataloc_params.jaabadetectparamsfilestr);
params = ReadParams(paramsfile);

%% log file

if isfield(dataloc_params,'jaabadetect_logfilestr') && ~DEBUG,
  logfile = fullfile(expdir,dataloc_params.jaabadetect_logfilestr);
  logfid = fopen(logfile,'a');
  if logfid < 1,
    warning('Could not open log file %s\n',logfile);
    logfid = 1;
  end
else
  logfid = 1;
end

timestamp = datestr(now,'yyyymmddTHHMMSS');
real_analysis_protocol = GetRealAnalysisProtocol(analysis_protocol,settingsdir);

jaabapath = fileparts(mfilename('fullpath'));
versionfile = fullfile(jaabapath,'version.txt');
jaaba_version = 'unknown';
if exist(versionfile,'file'),
  versionfid = fopen(versionfile,'r');
  if versionfid > 0,
    while true,
      s = fgetl(versionfid);
      if ~ischar(s),
        break;
      end
      s = strtrim(s);
      if ~isempty(s),
        jaaba_version = s;
        break;
      end
    end
    fclose(versionfid);
  end
end

fprintf(logfid,'\n\n***\nRunning FlyBowlJAABADetect version %s (JAABA version %s), analysis_protocol %s (linked to %s) at %s\n',version,jaaba_version,analysis_protocol,real_analysis_protocol,timestamp);

%% main loop

%jaabadir = fullfile(settingsdir,analysis_protocol,params.jaabadir);

%
% oldpwd = pwd;
% cd(params.jaabadir);

% loop through the classifier params files to process in order
if ~iscell(params.classifierparamsfiles),
  params.classifierparamsfiles = {params.classifierparamsfiles};
end

classifierinfo = [];
for i = 1:numel(params.classifierparamsfiles),

  classifierparamsfile = fullfile(settingsdir,analysis_protocol,params.classifierparamsfiles{i});
  classifierinfocurr = JAABADetect(expdir,'jablistfile',classifierparamsfile,...
    'forcecompute',forcecompute,...
    'debug',DEBUG);
  %'isrelativepath',true,...
  %'fnsrelative',{'featureparamfilename'});
  classifierinfo = structappend(classifierinfo,classifierinfocurr);

end
%
% cd(oldpwd);

%% save info

if ~DEBUG,

  savefile = fullfile(expdir,dataloc_params.jaabadetectinfomatfilestr);
  if exist(savefile,'file'),
    try %#ok<TRYNC>
      delete(savefile);
    end
  end
  try
    save(savefile,'classifierinfo','timestamp','analysis_protocol','real_analysis_protocol','version');
  catch ME
    warning('Could not save jaabadetect info to file %s: %s',savefile,getReport(ME));
  end

end

%% close log

fprintf(logfid,'Finished running FlyBowlJAABADetect at %s.\n',datestr(now,'yyyymmddTHHMMSS'));

if logfid > 1,
  fclose(logfid);
end

author: kristinbranson

repo: JAABA


text:
function [paramList, paramPath, addCur] = getParamList(param, bottomFlag)
% Get all the end parameters and their paths for matching CPLEX parameters appropriately
% (e.g., if param.simplex.display is a parameter, then we will have 'display'
% in paramList and 'param.simplex' in paramPath)
%
% USAGE:
%
%    [paramList, paramPath, addCur] = getParamList(param, bottomFlag)
%
% INPUTS:
%    param:           existing structure with parameters for Cplex
%    bottomFlag:      boolean switch to extract parameters based on 'Cur' at the
%                     bottom level in the user-supplied parameter structure
%
% OUTPUTS:
%    paramList:       Cell with proper parameter names
%    paramPath:       Cell with superseeding parameter structure
%    addCur:          Structure with booleans whether or not `.Cur` has been
%                     added to the `paramPath` for user-supplied parameters
%

structCur = param;  % current structure
lv = 1;  % current level of the structure
% use a stack structure to store the found parameters
lvFieldN = zeros(10,1);
lvFieldN(1) = 1;  % the number of fields being checked at the current level
lvField = cell(10, 1);
lvField{lv} = fieldnames(structCur);  % all fields in the current level
paramPath = {};
paramList = {};
addCur = [];
while lv > 0
    if isstruct(structCur.(lvField{lv}{lvFieldN(lv)}))
        % if the current level is still a structure
        if ~isempty(fieldnames(structCur.(lvField{lv}{lvFieldN(lv)})))
            % if it is not an empty structure, go into that structure
            structCur = structCur.(lvField{lv}{lvFieldN(lv)});
            lv = lv + 1;  % level + 1
            % start checking from the first field in the next level of structure
            lvFieldN(lv) = 1;
            % get all fields for the next structure
            lvField{lv} = fieldnames(structCur);
        else
            % an empty structure, go back to the closest level with any unchecked structures
            while lvFieldN(lv) == numel(lvField{lv})
                lv = lv - 1;
                if lv == 0
                    % finish if level = 0
                    break
                end
            end
            if lv > 0
                % check the next structure in the current level
                lvFieldN(lv) = lvFieldN(lv) + 1;
                structCur = param;
                for j = 1:lv-1
                    structCur = structCur.(lvField{j}{lvFieldN(j)});
                end
            end
        end
    else
        addCurLv = false;
        if ~bottomFlag || strcmpi(lvField{lv}{lvFieldN(lv)}, 'Cur')
            % if the current level is not a structure, with the bottomFlag off
            % (not going to the bottom level, for LP.Param), or if the current level
            % is named 'Cur' (may happen in the user-supplied parameter structure),
            % decrease one level since the bottom level of all parameters is either
            % 'Min', 'Max', 'Cur', 'Def', 'Name', 'Help' containing the information
            % for the parameter, which is the name for the structure immediately on
            % top of the current level
            lv = lv - 1;
            if bottomFlag  % 'Cur' at the bottom level in the user-supplied parameters
                addCurLv = true;
            end
        end
        if lv > 0
            % get the path for the current parameter
            c = {};
            for j = 1:lv
                c = [c lvField{j}(lvFieldN(j))];
            end
            paramPath = [paramPath; strjoin(c,'.')];
            % whether or not need to add back .Cur to the path for user-supplied parameters
            if addCurLv
                addCur = [addCur; true];
            else
                addCur = [addCur; false];
            end
            paramList = [paramList; c(end)];
            % go back to the closest level with any unchecked structures
            while lvFieldN(lv) == numel(lvField{lv})
                lv = lv - 1;
                if lv == 0
                    % finish if level = 0
                    break
                end
            end
            if lv > 0
                % check the next structure in the current level
                lvFieldN(lv) = lvFieldN(lv) + 1;
                structCur = param;
                for j = 1:lv-1
                    structCur = structCur.(lvField{j}{lvFieldN(j)});
                end
            end
        else
            % if lv = 0 because of bottomFlag, continue at lv = 1
            lv = 1;
            if lvFieldN(lv) == numel(lvField{lv})
                break
            else
                lvFieldN(1) = lvFieldN(1) + 1;
            end
        end
    end
end
end

author: opencobra

repo: cobratoolbox


text:
function aedat = ResampleSpace(aedat, newX, newY)

%{
This function take a structure which has been imported and squeezes or
stretches the events (polarity only for now ...)
into the params newX and newY.
parameters give the array size and one-based, but the address space is zero based.
%}

dbstop if error

if isfield(aedat, 'data') && isfield(aedat.data, 'polarity')
    % There should be a more principled way of finding out the array size,
    % but for now just look at the max addresses ...
    mx = uint64(max(aedat.data.polarity.x));
	temp = uint64(aedat.data.polarity.x) * (newX - 1);
	aedat.data.polarity.x = uint16(temp / mx);

    mx = uint64(max(aedat.data.polarity.y));
	temp = uint64(aedat.data.polarity.y) * (newY - 1);
	aedat.data.polarity.y = uint16(temp / mx);
end

author: panpanfei


text:
% StackExchange Signal Processing Q79314
% https://dsp.stackexchange.com/questions/79314
% Image Segmentation Using Deep Learning
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     27/11/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Constants

SIG_TYPE_RAND       = 1;
SIG_TYPE_GAUSSIAN   = 2;

%% Simulation Parameters

% Signal
signalType      = SIG_TYPE_GAUSSIAN;
numSamplesSig   = 201; %<! Make sure it is odd
gaussianStd     = 2;

% Support (Recieved Signal)
numSamples  = 300;

% Noise
minNoiseStd = 0;
maxNoiseStd = 1;
numNoisePts = 50;
numRealizations = 200; %<! Per noise level

%% Generate Data

signalRadius = ceil(4 * gaussianStd);
vSignalSupport = linspace(-signalRadius, signalRadius, numSamplesSig);
vSignal = exp(-(vSignalSupport .^ 2) ./ (2 * gaussianStd));
vSignal = single(vSignal(:));

figure(); plot(vSignalSupport, vSignal);

vNoiseLevel = linspace(minNoiseStd, maxNoiseStd, numNoisePts);
vNoiseLevel = single(vNoiseLevel);

signalRadius    = (numSamplesSig - 1) / 2;
numClasses      = numSamples - numSamplesSig + 1;

mData   = zeros(numSamples, numClasses * numNoisePts * numRealizations, 'single');
vLabels = zeros(numClasses * numNoisePts * numRealizations, 1);

vNoise  = zeros(numSamples, 1, 'single');
vX      = zeros(numSamples, 1, 'single');
sigIdx = 1;

for ii = 1:numClasses
    startIdx = ii;
    endIdx = startIdx + numSamplesSig - 1;
    vX(startIdx:endIdx) = vSignal;
    for jj = 1:numNoisePts
        noiseStd = vNoiseLevel(jj);
        for kk = 1:numRealizations
            vNoise(:) = noiseStd * randn(numSamples, 1, 'single');
            mData(:, sigIdx) = vX + vNoise;
            vLabels(sigIdx) = ii;
            sigIdx = sigIdx + 1;
        end
    end
    vX(startIdx:endIdx) = 0;
end

% ii = 20001; figure(); plot(mData(:, ii)); title(vLabels(ii));

% Verify data by generating it without noise.
% Then:
% [vA, vB] = max(mData, [], 1);
% vC = vB - signalRadius;
% isequal(vC(:), vLabels)

%% Save Data

save('Data', 'vSignal', 'mData', 'vLabels', 'subStreamNumber');

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function result = cvRsimpls(x,y,kmax,rmsecv,h,k)

%CVRIMPLS calculates the robust RMSECV (root mean squared error of cross-validation) curve
% for RSIMPLS or the robust RMSEP(root mean squared error of prediction) value in a fast way.
% The R-RMSECV curve can be used to make a selection of the optimal number of
% components to include in the regression model. The function is used in rsimpls.m.
%
% Input arguments:
%    x          : the explanatory variables
%    y          : the rsponse variables
%    kmax       : maximal number of variables to include in the model.
%    rmsecv     : Optional. If equal to 1 (default), the rmsecv is computed.
%                 Else, rmsecv = 0 and then the rss and R2 are computed.
%    h          : optional input argument.
%    k          : optional input argument. If k = 0 (default) then RMSECV is calculated. Else the RMSEP wil be computed.
%
% Output:
%  if RMSECV is computed:
%   result.rmsecv     : the R-RMSECV values (obtained with the minimum weights).
%  if RMSEP is computed:
%   result.rmsep      : the R-RMSEP values
%   result.rss        : the RSS values for every k.
%   result.R2         : the coefficient of determination for every k.
%   result.residu     : the residuals for every k = 1,...,kmax
%   result.outWeights : the weights used to compute the robust R-RMSEP values.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Sanne Engelen
% Last Update: 08/10/2004, 03/07/2006
% Last Revision: 03/07/2006

% some initialisations
n = size(x,1);
p = size(x,2);
q = size(y,2);
r = rank(x);
rz = rank([x,y]);
teller_if_lus = 0;
cutoffWeights = sqrt(chi2inv(0.975,q));

if nargin < 4
    alfa = 0.75;
    kmaxr=min([kmax+q,rz]);
    h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*alfa);
    k = 0;
    rmsecv = 1;
elseif nargin == 4
    alfa = 0.75;
    kmaxr=min([kmax+q,rz]);
    h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*alfa);
    k = 0;
elseif nargin == 5
    k = 0;
end

outWeights = weightcvRsimpls(x,y,kmax,h,k,cutoffWeights);

if rmsecv
    if k == 0
        w_min = outWeights.w_min;
    end
    rob = outWeights.ResRob;

    % Assigning the input variables
    if size(rob.Hsubsets.H0,2)==1
        rob.Hsubsets.H0=rob.Hsubsets.H0';
    end
    Hsets = [rob.Hsubsets.H0;rob.Hsubsets.H1;rob.Hsubsets.Hfreq];
    same.value = 0;
    data = [x,y];

    for i = 1:n
        disp(['observation ',num2str(i),' is left out'])
        X_min_i = removal(x,i,0);
        Y_min_i = removal(y,i,0);
        data_min_i = removal(data,i,0);

        same.value = 0;
        if isempty(find(rob.Hsubsets.H0 == i))
            if teller_if_lus >= 1
                same.value = 1;
            end
            teller_if_lus = teller_if_lus + 1;
        end

        % constructing Hsets of right size: h - 1.
        Hsets_min_i = RemoveObsHsets(Hsets,i);

        if k == 0
            res = removeObsRobpca(data,i,kmax + q,Hsets_min_i,same,1);
        else
            res = removeObsRobpca(data,i,k + q,Hsets_min_i,same);
        end

        if isempty(find(rob.Hsubsets.H0 == i))
            same.res = res;
        end

        Prob_min_i = res.Pk_min_i;
        Lrob_min_i = res.Lk_min_i;
        murob_min_i = res.muk_min_i;
        Trob_min_i = (data_min_i - repmat(murob_min_i,n-1,1))*Prob_min_i;

        % Computing weights corresponding with the ROBPCA results.
        sdrob_min_i = sqrt(libra_mahalanobis(Trob_min_i,zeros(1,size(Trob_min_i,2)),'invcov',1./Lrob_min_i))';

        if k == 0
            cutoff.sd=sqrt(chi2inv(0.975,kmax));
        else
            cutoff.sd = sqrt(chi2inv(0.975,k));
        end

        % Orthogonal distances to robust PCA subspace
        XRc=data_min_i-repmat(murob_min_i,n-1,1);
        Xtilde=Trob_min_i*Prob_min_i';
        Rdiff=XRc-Xtilde;
        for j=1:(n-1)
            odrob_min_i(j,:)=norm(Rdiff(j,:));
        end
        % Robust cutoff-value for the orthogonal distance
        if k == 0
            test_k = kmax;
        else
            test_k = k;
        end

        if test_k~=r
            [m,s]=unimcd(odrob_min_i.^(2/3),h);
            cutoff.od = sqrt(norminv(0.975,m,s).^3);
            wrob_min_i = (odrob_min_i<=cutoff.od)&(sdrob_min_i<=cutoff.sd);
        else
            cutoff.od=0;
            wrob_min_i = (sdrob_min_i<=cutoff.sd);
        end

        % start the deflation:
        xycentr_min_i = [];
        sigmax_min_i = [];
        xcentr_min_i = [];
        sigmaxy_min_i = [];
        sigmayx_min_i = [];

        xycentr_min_i = murob_min_i;
        sigmax_min_i = Prob_min_i(1:p,:)*diag(Lrob_min_i)*Prob_min_i(1:p,:)';
        xcentr_min_i = X_min_i - repmat(murob_min_i(1:p),n-1,1);
        sigmaxy_min_i = Prob_min_i(1:p,:)*diag(Lrob_min_i)*Prob_min_i(p+1:p+q,:)';
        sigmayx_min_i = sigmaxy_min_i';

        % calculation of the scores.
        nScores = 1;
        R_min_i = [];
        T_min_i = [];
        P_min_i = [];
        V_min_i = [];

        if k == 0
            countScores = kmax;
        else
            countScores = k;
        end

        while nScores <= countScores
            if q == 1
                qq_min_i = 1;
            else
                [QQ,LL] = eig(sigmayx_min_i*sigmaxy_min_i);
                [LL,I] = greatsort(diag(LL));
                qq_min_i = QQ(:,I(1));
            end

            rr_min_i = sigmaxy_min_i*qq_min_i;
            rr_min_i = rr_min_i/norm(rr_min_i);
            tt_min_i = xcentr_min_i*rr_min_i;
            pp_min_i = sigmax_min_i*rr_min_i/(rr_min_i'*sigmax_min_i*rr_min_i);
            vv_min_i = pp_min_i;

            if nScores > 1
                vv_min_i = vv_min_i - V_min_i*(V_min_i'*pp_min_i);
            end
            vv_min_i = vv_min_i./norm(vv_min_i);
            sigmaxy_min_i = sigmaxy_min_i - vv_min_i*(vv_min_i'*sigmaxy_min_i);

            V_min_i(:,nScores) = vv_min_i;
            T_min_i(:,nScores) = tt_min_i;
            R_min_i(:,nScores) = rr_min_i;
            P_min_i(:,nScores) = pp_min_i;

            nScores = nScores + 1;
        end

        if k == 0
            outRegr = runRegr(x,y,i,T_min_i,Y_min_i,R_min_i,murob_min_i,wrob_min_i,kmax,cutoffWeights);
            for j = 1:kmax
                Tk_min_i = T_min_i(:,1:j);
                geg = [Tk_min_i,Y_min_i];

                outRegr.Mu = outRegr.center';
                outRegr.Sigma = outRegr.sigma;
                [Bk,intk,sigmayykmaxrew_k,sigmattkmaxrew_k] = extractmcdregres(outRegr,Tk_min_i,Y_min_i,kmax,n-1,q,j,h-1,cutoffWeights);
                coeffk = [Bk;intk];
                b_min_i = R_min_i(:,1:j)*coeffk(1:j,:);
                int_min_i = coeffk(j+1,:)  - murob_min_i(1:p)*R_min_i(:,1:j)*coeffk(1:j,:);
                Yhat_min_i = x(i,:)*b_min_i + int_min_i;
                resid_min_i(i,(j-1)*q + 1:j*q) = y(i,:) - Yhat_min_i;

                % calculation of the resd:
                rewE2=sigmayykmaxrew_k- coeffk(1:j,1:q)'*sigmattkmaxrew_k*coeffk(1:j,1:q);

                if q > 1
                    cov = rewE2;
                    cen=zeros(q,1);
                    resd(i,j)=sqrt(libra_mahalanobis(resid_min_i(i,(j-1)*q + 1:j*q),cen','cov',cov))'; %robust distances of residuals
                else
                    scale = sqrt(rewE2);
                    resd(i,j) = resid_min_i(i,(j-1)*q + 1:j*q)/scale;
                end
            end
        else
            outRegr = runRegr(x,y,i,T_min_i,Y_min_i,R_min_i,murob_min_i,wrob_min_i,k,cutoffWeights);
            resid_min_i(i,:) = outRegr.resid_min_i;
            resd(i,:) = outRegr.resd';
        end
    end

    if k == 0
        for j = 1:kmax
            resk = resid_min_i(:,(j-1)*q + 1:j*q);
            if q == 1
                rmsecv(j) = sqrt(1/sum(w_min)*w_min*(resk).^2);
            else
                rmsecv(j) = sqrt(1/sum(w_min)*w_min*(mean((resk').^2))');
            end
        end
        result.rmsecv = rmsecv;
        result.residu = resid_min_i;
    else
        weights = outWeights.weightsk;
        if q == 1
            rmsep = sqrt(1/sum(weights)*weights'*(resid_min_i).^2);
        else
            rmsep = sqrt(1/sum(weights)*weights'*(mean((resid_min_i').^2))');
        end
        result.rmsep = rmsep;
        result.residu = resid_min_i;
    end
end

result.outWeights = outWeights;
result.rss = outWeights.rss;
result.R2 = outWeights.R2;

%------------------------------------------------------------------
function out = runRegr(x,y,i,T_min_i,Y_min_i,R_min_i,mukmax_min_i,wkmax_min_i,k,cutoffWeights)

[n,p] = size(x);
[n,q] = size(y);

% perform the robpca regression:
breg = [];
b_min_i = [];
int_min_i= [];
Yhat_min_i = [];
robpcareg = robpcaregres(T_min_i(:,1:k),Y_min_i,wkmax_min_i',cutoffWeights);
out.center = robpcareg.center;
out.sigma = robpcareg.sigma;
breg = robpcareg.coeffs(1:k,:);
b_min_i = R_min_i(:,1:k)*breg;
int_min_i = robpcareg.coeffs(k+1,:) - mukmax_min_i(1:p)*R_min_i(:,1:k)*breg;
Yhat_min_i = x(i,:)*b_min_i + int_min_i;
resid_min_i = y(i,:) - Yhat_min_i;

% calculation of the resd:
if q > 1
    cov = robpcareg.cov;
    cen=zeros(q,1);
    resd=sqrt(libra_mahalanobis(resid_min_i,cen','cov',cov))'; %robust distances of residuals
else
    scale = sqrt(robpcareg.cov);
    resd = resid_min_i/scale;
end

out.resid_min_i = resid_min_i;
out.resd = resd;
%----------------------------------------------------------------------------------------
function out = weightcvRsimpls(x,y,kmax,h,k,cutoffWeights)

% Computes the weights for the robust RMSECV/RMSEP values.
%
% input:
%   x    : the independent variables.
%   y    : the response variables.
%   kmax : the maximal number of components to be considered.
%   h    : the number of observations on which the calculations are based.
%   k    : if equal to zero, robpca is performed on kmax components (case RMSECV). (default).
%          Else, robpca is performed for a certain number of components (case RMSEP)
%
% output:
%   out.w_min    : the weights obtained by taking the minimum over all k
%   out.weightsk : the weights for all observations and all k (n x kmax)
%   out.resrob   : the results of robpca on [x,y].
%   out.R2       : the weighted Rsquared for each value of k
%   out.rss      : the weighted rss for each value of k

n = size(x,1);
p = size(x,2);
q = size(y,2);
r = rank(x);

if nargin < 5
    k = 0;
end

if k == 0
    ResRob = robpca([x,y],'plots',0,'k',kmax + q,'kmax',kmax + q,'h',h);
else
    ResRob = robpca([x,y],'plots',0,'k',k + q,'kmax',kmax + q,'h',h);
end

Trob = ResRob.T;
Prob = ResRob.P;
Lrob = ResRob.L;
murob = ResRob.M;
wrob = ResRob.flag.all;

%deflation

xycentr = [];
sigmax = [];
xcentr = [];
sigmaxy = [];
sigmayx = [];

xycentr = murob;
sigmax = Prob(1:p,:)*diag(Lrob)*Prob(1:p,:)';
xcentr = x - repmat(murob(1:p),n,1);
sigmaxy = Prob(1:p,:)*diag(Lrob)*Prob(p+1:p+q,:)';
sigmayx = sigmaxy';

% calculation of the scores.
nScores = 1;
R = [];
T = [];
P = [];
V = [];

if k == 0
    countScores = kmax;
else
    countScores = k;
end

while nScores <= countScores
    if q == 1
        qq = 1;
    else
        [QQ,LL] = eig(sigmayx*sigmaxy);
        [LL,I] = greatsort(diag(LL));
        qq = QQ(:,I(1));
    end

    rr = sigmaxy*qq;
    rr = rr/norm(rr);
    tt = xcentr*rr;
    pp = sigmax*rr/(rr'*sigmax*rr);
    vv = pp;

    if nScores > 1
        vv = vv - V*(V'*pp);
    end
    vv = vv./norm(vv);
    sigmaxy = sigmaxy - vv*(vv'*sigmaxy);

    V(:,nScores) = vv;
    T(:,nScores) = tt;
    R(:,nScores) = rr;
    P(:,nScores) = pp;

    nScores = nScores + 1;
end

if k == 0
    outRobRegr = robpcaregres(T(:,1:kmax),y,wrob');

    for j = 1:kmax
        outRobRegr.Mu = outRobRegr.center';
        outRobRegr.Sigma = outRobRegr.sigma;
        [Bk,intk,sigmayykmaxrew_k,sigmattkmaxrew_k] = extractmcdregres(outRobRegr,T(:,1:j),y,kmax,n,q,j,h,cutoffWeights);
        coeffk = [Bk;intk];
        finalB = R(:,1:j)*coeffk(1:j,:);
        finalInt = coeffk(j+1,:)  - murob(1:p)*R(:,1:j)*coeffk(1:j,:);
        Yhat = x*finalB + repmat(finalInt,n,1);
        resid(:,(j-1)*q+1:j*q) = y - Yhat;

        % calculation of the rd:
        rewE2=sigmayykmaxrew_k- coeffk(1:j,1:q)'*sigmattkmaxrew_k*coeffk(1:j,1:q);

        if q > 1
            cov = rewE2;
            cen=zeros(q,1);
            resd = sqrt(libra_mahalanobis(resid(:,(j-1)*q+1:j*q),cen','cov',cov))'; %robust distances of residuals
            weightsk(:,j) = (abs(resd)<=cutoffWeights);
        else
            scale = sqrt(rewE2);
            resd = resid(:,(j-1)*q+1:j*q)/scale;
            weightsk(:,j) = (abs(resd)<=cutoffWeights);
        end
    end
else
    outRobRegr = robpcaregres(T(:,1:k),y,wrob');

    % robust residual distance:
    if q==1
        resd=outRobRegr.resids/sqrt(outRobRegr.cov);
    else
        resd=sqrt(libra_mahalanobis(outRobRegr.resids,zeros(1,q),'cov',outRobRegr.cov))';
    end

    weightsk = (abs(resd)<=cutoffWeights);
end

if k == 0
    if kmax == 1
        w_min = weightsk';
    else
        w_min = min(weightsk');
    end

    out.w_min = w_min;
    out.weightsk = weightsk;

    yw = mean(y(w_min == 1,:));
    y2=(y-repmat(yw,n,1)).^2;
    R = resid.^2;
    D=sum(y2(w_min==1,:));
    for j = 1:kmax
        R1=R(w_min==1,(j-1)*q+1:j*q);
        rss(j) = sum(sum(R1));
        R2(j)=1-rss(j)/sum(D);
    end

    out.rss = 1/(q*sum(w_min))*rss;
    out.R2 = R2;
else
    out.weightsk = weightsk;

    s=sum(weightsk);
    yw=sum(y(weightsk==1,:))/s;
    y2=(y-repmat(yw,n,1)).^2;
    R = outRobRegr.resids.^2;
    D=sum(y2(weightsk==1,:));
    R1=R(weightsk==1,:);
    rss = sum(sum(R1));
    R2=1-rss/sum(D);

    out.rss = 1/(q*sum(weightsk))*rss;
    out.R2 = R2;
end
out.ResRob = ResRob;
%---------------------------------------------------------------------------------------
function Hsets_min_i = RemoveObsHsets(Hsets,i)

% removes the right index from the $h$-subsets in Hsets to
% obtain (h - 1)-subsets.
% every h-set is put as a row in Hsets.
% i is the index of the observation that is removed from the whole data.

for r = 1:size(Hsets,1)
    if ~isempty(find(Hsets(r,:)== i))
        Hsets_min_i(r,:) = removal(Hsets(r,:),0,find(Hsets(r,:) == i));
    else
        Hsets_min_i(r,:) = Hsets(r,1:(end-1));
    end

    for j = 1:length(Hsets_min_i(r,:))
        if Hsets_min_i(r,j) > i
            Hsets_min_i(r,j) = Hsets_min_i(r,j) - 1;
        end
    end
end

author: marianux

repo: ecg-kit

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


text:
function [mesh, cfg] = ft_prepare_mesh(cfg, data)

% FT_PREPARE_MESH creates a triangulated surface mesh or tetrahedral/hexahedral
% volume mesh that can be used as geometrical description for a volume conduction
% model. The mesh can either be created manually from anatomical MRI data or can be
% generated starting from a segmented MRI. This function can also be used to create a
% cortex hull, i.e. the smoothed envelope around the pial surface created by
% freesurfer.
%
% Use as
%   mesh = ft_prepare_mesh(cfg)
%   mesh = ft_prepare_mesh(cfg, mri)
%   mesh = ft_prepare_mesh(cfg, seg)
% where the mri input argument is the result from FT_READ_MRI, FT_VOLUMEREALIGN or
% FT_VOLUMERESLICE and the seg input argument is from FT_VOLUMESEGMENT. If you
% specify an anatomical MRI, it will be segmented on the fly.
%
% The cfg argument is a structure that can contain:
%   cfg.method      = string, can be 'interactive', 'projectmesh', 'iso2mesh', 'isosurface',
%                     'headshape', 'hexahedral', 'tetrahedral', 'cortexhull' or 'fittemplate'
%   cfg.tissue      = cell-array with strings representing the tissue types, or numeric vector with integer values
%   cfg.numvertices = numeric vector, should have same number of elements as the number of tissues
%
% When providing an anatomical MRI or a segmentation, you should specify
%   cfg.downsample  = integer number (default = 1, i.e. no downsampling), see FT_VOLUMEDOWNSAMPLE
%   cfg.spmversion  = string, 'spm2', 'spm8', 'spm12' (default = 'spm12')
%
% For method 'headshape' you should specify
%   cfg.headshape   = a filename containing headshape, a Nx3 matrix with surface
%                     points, or a structure with a single or multiple boundaries
%
% For method 'cortexhull' you should not give input data, but specify
%   cfg.headshape   = string, filename containing the pial surface computed by freesurfer recon-all
%
% For method 'fittemplate' you should specify
%   cfg.headshape   = a filename containing headshape
%   cfg.template    = a filename containing headshape
% With this method you are fitting the headshape from the configuration to the template;
% the resulting affine transformation is applied to the input mesh (or set of meshes),
% which is subsequently returned as output variable.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% Example
%   mri             = ft_read_mri('Subject01.mri');
%
%   cfg             = [];
%   cfg.output      = {'scalp', 'skull', 'brain'};
%   segmentation    = ft_volumesegment(cfg, mri);
%
%   cfg             = [];
%   cfg.tissue      = {'scalp', 'skull', 'brain'};
%   cfg.numvertices = [800, 1600, 2400];
%   mesh            = ft_prepare_mesh(cfg, segmentation);
%
%   cfg             = [];
%   cfg.method      = 'cortexhull';
%   cfg.headshape   = '/path/to/surf/lh.pial';
%   cfg.fshome      = '/path/to/freesurfer dir';
%   cortex_hull     = ft_prepare_mesh(cfg);
%
% See also FT_VOLUMESEGMENT, FT_PREPARE_HEADMODEL, FT_PLOT_MESH

% Undocumented functionality: at this moment it allows for either
%   mesh = ft_prepare_mesh(cfg)             or
%   mesh = ft_prepare_mesh(cfg, headmodel)
% but more consistent would be to specify a volume conduction model with
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%   cfg.headshape     = name of file containing the volume conduction model, see FT_READ_HEADMODEL
%
% Undocumented options, I have no clue why they exist
%   cfg.method = {'singlesphere' 'concentricspheres' 'localspheres'}

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% we cannot use nargin, because the data might have been loaded from cfg.inputfile
hasdata = exist('data', 'var');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden', {'numcompartments', 'outputfile', 'sourceunits', 'mriunits'});

% get the options
cfg.downsample  = ft_getopt(cfg, 'downsample', 1);            % default is no downsampling
cfg.numvertices = ft_getopt(cfg, 'numvertices', 3000, true);  % set the default, [] is also a meaningful value
cfg.smooth      = ft_getopt(cfg, 'smooth');                   % no default
cfg.spmversion  = ft_getopt(cfg, 'spmversion', 'spm12');

% Translate the input options in the appropriate default for cfg.method
if isfield(cfg, 'headshape') && ~isempty(cfg.headshape)
  cfg.method = ft_getopt(cfg, 'method', 'headshape');
elseif hasdata && ~strcmp(ft_headmodeltype(data), 'unknown')
  cfg.method = ft_getopt(cfg, 'method', ft_headmodeltype(data));
elseif hasdata
  cfg.method = ft_getopt(cfg, 'method', 'projectmesh');
else
  cfg.method = ft_getopt(cfg, 'method', []);
end

if hasdata && cfg.downsample~=1
  % optionally downsample the anatomical volume and/or tissue segmentations
  tmpcfg = keepfields(cfg, {'downsample', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  data = ft_volumedownsample(tmpcfg, data);
  % restore the provenance information and put back cfg.smooth
  tmpsmooth = cfg.smooth;
  [cfg, data] = rollback_provenance(cfg, data);
  cfg.smooth = tmpsmooth;
end

switch cfg.method
  case 'interactive'
    % this makes sense with a non-segmented MRI as input
    % call the corresponding helper function
    mesh = prepare_mesh_manual(cfg, data);

  case {'projectmesh', 'iso2mesh', 'isosurface'}
    % this makes sense with a segmented MRI as input
    % call the corresponding helper function
    mesh = prepare_mesh_segmentation(cfg, data);

  case 'headshape'
    % call the corresponding helper function
    mesh = prepare_mesh_headshape(cfg);

  case 'hexahedral'
    % the MRI is assumed to contain a segmentation
    % call the corresponding helper function
    mesh = prepare_mesh_hexahedral(cfg, data);

  case 'tetrahedral'
    % the input data is assumed to contain a segmentation
    % call the corresponding helper function
    mesh = prepare_mesh_tetrahedral(cfg, data);

  case {'singlesphere' 'concentricspheres'}
    % the input data is assumed to contain a spherical head model
    data = ft_datatype_headmodel(data);   % ensure that it is consistent and up-to-date
    data = ft_determine_units(data);      % ensure that it has units

    if length(data.r)>1 && numel(cfg.numvertices)==1
      % use the same number of vertices for each tissue
      cfg.numvertices = repmat(cfg.numvertices, length(data.r));
    end

    mesh = [];
    for i=1:length(data.r)
      ft_info('triangulating sphere %d in the volume conductor\n', i);
      [pos, tri] = mesh_sphere(cfg.numvertices(i));
      mesh(i).pos(:,1) = pos(:,1)*data.r(i) + data.o(1);
      mesh(i).pos(:,2) = pos(:,2)*data.r(i) + data.o(2);
      mesh(i).pos(:,3) = pos(:,3)*data.r(i) + data.o(3);
      mesh(i).tri = tri;
    end

  case 'cortexhull'
    mesh = prepare_mesh_cortexhull(cfg);

  case 'fittemplate'
    M    = prepare_mesh_fittemplate(cfg.headshape.pos, cfg.template.pos);
    mesh = ft_transform_geometry(M, data);

  otherwise
    ft_error('unsupported cfg.method')
end

% copy the geometrical units from the input to the output
if ~isfield(mesh, 'unit') && hasdata && isfield(data, 'unit')
  for i=1:numel(mesh)
    mesh(i).unit = data.unit;
  end
elseif ~isfield(mesh, 'unit')
  mesh = ft_determine_units(mesh);
end

% copy the coordinate system from the input to the output
if ~isfield(mesh, 'coordsys') && hasdata && isfield(data, 'coordsys')
  for i=1:numel(mesh)
    mesh(i).coordsys = data.coordsys;
  end
end

% smooth the mesh
if ~isempty(cfg.smooth)
  tmpcfg = keepfields(cfg, {'smooth', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.numvertices = []; % the number of vertices should not be changed
  tmpcfg.headshape = mesh;
  mesh = prepare_mesh_headshape(tmpcfg);
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data
ft_postamble provenance mesh
ft_postamble history    mesh


text:
function determ = wilk04_determinant ( )

%*****************************************************************************80
%
%% WILK04_DETERMINANT returns the determinant of the WILK04 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real DETERM, the determinant.
%
  determ = 0.9143E-04 * 0.7156E-04 * 0.9504E-04 * 0.7123E-04;

  return
end


text:
function anonStructS = anonymizeCERRdataStruct(structS,anonDefinS)
% function anonStructS = anonymizeCERRdataStruct(structS,anonDefinS)
%
% Anonymizes the passed CERR/Matlab data structure according to the
% anonymization definitions in anonDefinS.
%
% The anonymization is applied to the fields of structS structure array if
% the field names are anything other than Item_*. If the fields of structS
% are Item_* (Item_1, Item_2, Item_3), the anonymization is applied to the
% elements of structS.item_*.
%
% APA, 1/29/2018

% initialize the anonStructS data structure to be same as structS
anonStructS = structS;

% Get the field names of structS
fieldC = fieldnames(anonStructS);

% If the fields are Item_*, then anonymize its sub-fields
itemIteratorFlag = 0;
if length(anonStructS) == 1 && isequal(unique(strtok(fieldC,'_')),{'Item'})
    itemIteratorFlag = 1;
    itemC = fieldC;
    fieldC = fieldnames(anonStructS.Item_1);
end

% Loop over all the fields and anonymize
for i = 1:length(fieldC)
    % remove field if it is not defined in the anonymization
    rmFieldFlag = 0;
    if ~itemIteratorFlag && ~isfield(anonDefinS,fieldC{i})
        %anonStructS = rmfield(structS,fieldC{i});
        rmFieldFlag = 1;
    elseif itemIteratorFlag && ~isfield(anonDefinS,fieldC{i})
        rmFieldFlag = 1;
    end

    if itemIteratorFlag
        % check string values in case of pre-defined strings
        for elemNum = 1:length(itemC)
            if ~rmFieldFlag
                valC = anonDefinS.(fieldC{i});
                if isfield(anonStructS.(itemC{elemNum}),(fieldC{i}))
                    fieldVal = anonStructS.(itemC{elemNum}).(fieldC{i});
                else
                    fieldVal = '';
                end
                fieldCalss = class(fieldVal);
            else
                fieldCalss = 'remove';
            end
            switch fieldCalss
                case {'double','logical','single','uint8','uint16','uint32' ...
                        ,'uint64','int8','int16','int32','int64'}
                    % do nothing to the numeric fields
                case 'char'
                    if isequal(valC,'keep')
                        % keep the field as is (UID fields)
                    elseif isequal(valC,'date')
                        % assign a dummy date
                        disp(strcat('Anonymizing date from "', fieldC{i}, '" field'))
                        anonStructS.(itemC{elemNum}).(fieldC{i}) = '11111111';
                    elseif ~isempty(fieldVal) && ~ismember(fieldVal,valC)
                        disp(strcat('Anonymizing "', fieldC{i}, '" field'))
                        anonStructS.(itemC{elemNum}).(fieldC{i}) = 'CERR anonymized';
                    end
                case 'struct'
                    if ~isstruct(anonDefinS.(fieldC{i})) && isequal(anonDefinS.(fieldC{i}),'keep')
                        % keep all the sub-fields of this data-structure
                    else
                        % anonymize sub-fields as per the definitions
                        anonStructS.(itemC{elemNum}).(fieldC{i}) = anonymizeCERRdataStruct...
                            (anonStructS.(itemC{elemNum}).(fieldC{i}), anonDefinS.(fieldC{i}));
                    end
                case 'remove'
                    anonStructS.(itemC{elemNum}) = ...
                        rmfield(anonStructS.(itemC{elemNum}),(fieldC{i}));
                case 'cell'

            end
        end
    else
        % check string values in case of pre-defined strings
        for elemNum = 1:length(structS)
            if ~rmFieldFlag
                valC = anonDefinS.(fieldC{i});
                if isfield(anonStructS(elemNum),(fieldC{i}))
                    fieldVal = anonStructS(elemNum).(fieldC{i});
                else
                    fieldVal = '';
                end
                fieldCalss = class(fieldVal);
            elseif isfield(anonStructS,fieldC{i})
                disp(strcat('Removing "', fieldC{i}, '" field'))
                anonStructS = rmfield(anonStructS,fieldC{i});
                continue;
            else
                continue;
            end
            switch fieldCalss
                case {'double','logical','single','uint8','uint16','uint32' ...
                        ,'uint64','int8','int16','int32','int64'}
                    % do nothing to the numeric fields
                case 'char'
                    if isequal(valC,'keep')
                        % keep the field as is (UID fields)
                    elseif isequal(valC,'date')
                        % assign a dummy date
                        disp(strcat('Anonymizing date from "', fieldC{i}, '" field'))
                        anonStructS(elemNum).(fieldC{i}) = '11111111';
                    elseif ~isempty(fieldVal) && ~ismember(fieldVal,valC)
                        disp(strcat('Anonymizing "', fieldC{i}, '" field'))
                        anonStructS(elemNum).(fieldC{i}) = 'CERR anonymized';
                    end
                case 'struct'
                    if ~isstruct(anonDefinS.(fieldC{i})) && isequal(anonDefinS.(fieldC{i}),'keep')
                        % keep all the sub-fields of this data-structure
                    else
                        % anonymize sub-fields as per the definitions
                        anonStructS(elemNum).(fieldC{i}) = anonymizeCERRdataStruct...
                        (anonStructS(elemNum).(fieldC{i}), anonDefinS.(fieldC{i}));
                    end

                case 'cell'

            end
        end
    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [  ] = convert2vid( vid_data, varargin )
%convert2vid - converts the video structure, as created by sie2mat.exe,
%and converts it to a video.
% Program is not optimised, but requires no video processing toolbox
%   Syntax
%       convert2vid( vid_data )  - returns a SIE_DATA structure with data from
%                               the sie file as defined in sie_flnme
%       convert2vid( vid_data, 'path', path_str )  -  path_str defines
%                               the location where the video is saved.
%                               Default - current folder
%       convert2vid( vid_data, 'video_name', video_name_str )  -  video_name
%                               changes the name of the output video
%                               Default - 'SIE_VIDEO'
%       convert2vid( vid_data, 'profile', profile_str )  -  profile
%                               String enclosed in single quotation marks that describes
%                               the type of file to create. Specifying a profile sets default
%                               values for video properties such as VideoCompressionMethod.
%                               Possible values:
%                                     'Archival' - Motion JPEG 2000 file with lossless compression
%                                     'Motion JPEG AVI' - Compressed AVI file using Motion JPEG codec
%                                     'Motion JPEG 2000' - Compressed Motion JPEG 2000 file
%                                     'MPEG-4' - Compressed MPEG-4 file with H.264 encoding (Windows 7 systems only)
%                                     'Uncompressed AVI' - Uncompressed AVI file with RGB24 video
%                               Default: 'Motion JPEG AVI'
%       convert2vid( vid_data, 'framerate', FrameRate_int )  - FrameRate
%                               Rate of playback for the video in frames per second.
%                               Default: 30
%
% convert2vid.m
% Created by: MJ Stallmann
% Last modified: 02/08/2013

%% Start sie2mat
% get options
[options] = parseOptions( varargin{:});
[pathstr, name, ext] = fileparts(options.video_name) ;
vid_name=fullfile(options.path,name,'.avi');
writerObj = VideoWriter(vid_name,options.profile);
writerObj.FrameRate=options.framerate;
open(writerObj);
nframes=length(vid_data);
disp('Processing frames.... Please wait...')
for k = 1:nframes
    fid = fopen('temp.jpg', 'w');
    fwrite(fid, vid_data(k).Video_data, 'uint8');
    fclose(fid);
    frame = imread('temp.jpg');
    writeVideo(writerObj,frame);
end

close(writerObj);
delete('temp.jpg');
end

%% Set options
function options = parseOptions(varargin)
IP = inputParser;
IP.addParamValue('path', pwd, @ischar)
IP.addParamValue('video_name','SIE_VIDEO', @ischar)
IP.addParamValue('profile','Motion JPEG AVI', @ischar)
IP.addParamValue('framerate',30, @isnumeric)
IP.parse(varargin{:});
options = IP.Results;
end