text:
function [f,nxe_train,w] = train_qfuser(scores,scrQ,key,quiet,w0,cstepHessian,niters,obj_func,prior)
% Trains a linear fusion of scores with quality measures starting
% from the linear fusion of the only the scores.
% Inputs:
%   scores: A matrix of dev scores where the number of rows is the
%     number of systems to be fused and the number of columns
%     is the number of scores.
%   scrQ: An object of type Quality containing the quality measures
%     for the dev trials.
%   key: A Key object indicating target and non-target trials.
%   quiet: Boolean.  If true, the training algorithm outputs fewer
%     messages describing its progress.
%   w0: The fusion weights from the linear (without quality) fusion
%     step.
%   cstepHessian: Boolean.  If true, the training algorithm will
%     calculate the Hessian matrix using complex step
%     differentiation which should make training faster.
%   niters: The maximum number of training iterations.
%   obj_func: The objective function for the fusion training.  If [],
%     cllr objective is used.
%   prior: The effective target prior.
% Outputs:
%   f: A function handle for the trained fusion function (of type qfuser_backoff).
%   nxe_train: The normalized cross-entropy of the solution found
%     when training the fusion weights.
%   w: The fusion weights found by the trainer.

assert(nargin==9)
assert(isa(key,'Key'))

trials = key.tar(:)' | key.non(:)';
n1 = size(scores,2);
n2 = sum(trials);
assert(n1==n2,'# score trials = %i, but # key trials = %i',n1,n2)

% dump bad trials
validscores = scrQ.scoremask(trials);
if sum(trials)-sum(validscores)>0
    log_warning('Quality fusion training has %i missing trials.\n',sum(trials)-sum(validscores));
end
scores = scores(:,validscores);

% remove models and segs without quality measures from key (and the
% corresponding trials)
key.modelset = key.modelset(scrQ.hasmodel);
key.segset = key.segset(scrQ.hasseg);
key.tar = key.tar(scrQ.hasmodel,scrQ.hasseg);
key.non = key.non(scrQ.hasmodel,scrQ.hasseg);
trials = key.tar(:)' | key.non(:)';

% also remove them from scrQ
scrQ.modelset = scrQ.modelset(scrQ.hasmodel);
scrQ.segset = scrQ.segset(scrQ.hasseg);
scrQ.scoremask = scrQ.scoremask(scrQ.hasmodel,scrQ.hasseg);
scrQ.modelQ = scrQ.modelQ(:,scrQ.hasmodel);
scrQ.segQ = scrQ.segQ(:,scrQ.hasseg);

scrQ.hasmodel = scrQ.hasmodel(scrQ.hasmodel);
scrQ.hasseg = scrQ.hasseg(scrQ.hasseg);

assert(size(scores,2)==sum(trials));

% set up function that must be trained
[fusion,w0] = qfuser_linear([],scores,scrQ,key,w0);

% Create vector that classifies scores as target/non-target
tar = key.tar(trials);
tar = tar(:)';
non = key.non(trials);
non = non(:)';
classf = zeros(1,size(scores,2));
classf(tar) = 1;
classf(non) = -1;
% rest is zero and will be ignored

% train the fusion weights
[w,nxe_train] = train_binary_classifier(fusion,classf,w0,obj_func,prior,[],0,niters,[],[],quiet,cstepHessian);

% create function handle to do fusion using trained weights
f = @(scores,scrQ,key,backoff_scores) qfuser_backoff(w,scores,scrQ,key,backoff_scores);


text:
function filteredOutS = sitkWrapper(sitkLibPath, scan3M, filterType, paramS)
%function filteredOutS = sitkWrapper(sitkLibPath, scan3M, filterType, paramS, planC)
%Calculate image filters using the Simple ITK Python Library
%--------------------------------------------------------------------------
%sitkLibPath - location of sitk python wrappscan to be filtered
%description - Short description string if desired.
%filterType - name of sitk filter
%             Supported options: 'GradientImageFilter','HistogramMatchingImageFilter',...
%             'LaplacianRecursiveGaussianImageFilter',...
%             'N4BiasFieldCorrectionImageFilter',
%             'N4BiasAndHistogramCorrectionImageFilter'.
%paramS: parameters required to calculate the filter
%planC: to convert scan back to cerr
%--------------------------------------------------------------------------
% example usage:
% filterType = 'GradientImageFilter';
% paramS.useImageSpacing = false;
% paramS.useImageDirection = true;
% sitkLibPath = 'C:\Python34\Lib\site-packages\SimpleITK\';
% planC = loadPlanC(cerrFileName,tempdir);
% planC = updatePlanFields(planC);
% % Quality assure
% planC = quality_assure_planC(cerrFileName,planC);
% indexS = planC{end};
%  % Get Scan
%  scanNum  =1;
%  scanM = double(planC{indexS.scan}(scanNum).scanArray) ...
%      - planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;%
%
% sitkWrapper(sitkLibPath, scanM, filterType, paramS)
%
%
% Rutu Pandya, Dec, 16 2019.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.

% import python module SimpleITK
sitkModule = 'SimpleITK';
numpyModule = 'numpy';
%p = pyenv;

P = py.sys.path;
currentPath = pwd;
cd(sitkLibPath);
sitkFileName = fullfile(sitkLibPath,sitkModule);
numpyFileName = fullfile(sitkLibPath,numpyModule);

try
    if count(P,numpyFileName) == 0
        insert(P,int32(0),numpyFileName);
    end
    np = py.importlib.import_module(numpyModule);
    if count(P,sitkFileName) == 0
        insert(P,int32(0),sitkFileName);
    end
    sitk = py.importlib.import_module(sitkModule);
catch
    disp('SimpleITK module could not be imported, check the path');
end

cd(currentPath);

%origScanSize = size(scanM);
% visualize original scan
%slc = 50;
%figure, imagesc(scanM(:,:,50)), title('orig Image')

% % convert scan to numpy array and integer
% scanPy = py.numpy.array(scanM(:).');
% scanPy = scanPy.astype(py.numpy.int64);
%
% % get original shape of the scan
% origShape = py.numpy.array(size(scanM));
% origShape = origShape.astype(py.numpy.int64);

switch filterType
    case 'GradientImageFilter'
        % paramS inputs needed:
        % useImageSpacing (bool), true by default
        % useImageDirection (bool), true by default

        % convert scan to numpy array and integer
        scanPy = np.array(scan3M);
        scanPy = scanPy.astype(py.numpy.float32);

        % get original shape of the scan
        %origShape = py.numpy.array(size(scanM));
        %origShape = origShape.astype(py.numpy.int64);

        % reshape numpy array to original shape
        %scanPy = reshape(scanPy,origShape);

        % Get image from the array
        itkimg = sitk.GetImageFromArray(scanPy);

        % calculate gradient
        gradient = sitk.GradientImageFilter();

        if(paramS.useImageSpacing == false)
            gradient.SetUseImageSpacing(paramS.useImageSpacing);
        end
        if(paramS.useImageDirection == false)
            gradient.SetUseImageDirection(paramS.useImageDirection);
        end
        gradImg = gradient.Execute(itkimg);

        % extract numpy array from resulting image
        npGradImg = sitk.GetArrayFromImage(gradImg);

        % convert resulting numpy array to matlab array in required shape
        %dblGradResultM = double(py.array.array('d',py.numpy.nditer(npGradImg)));
        %gradMatM = reshape(dblGradResultM,[3,origScanSize]);
        gradMatM = double(npGradImg);
        gradMatM = permute(gradMatM,[2,3,4,1]);

        %             %visualize
        %             size(gradMatM)
        %             figure, imagesc(gradMatM(:,:,50,1))
        %             figure, imagesc(gradMatM(:,:,50,2))
        %             figure, imagesc(gradMatM(:,:,50,3))

        filteredOutS.xGradient = gradMatM(:,:,:,1);
        filteredOutS.yGradient = gradMatM(:,:,:,2);
        filteredOutS.zGradient = gradMatM(:,:,:,3);

    case 'LaplacianRecursiveGaussianImageFilter'

        % convert scan to numpy array and integer
        % scanPy = py.numpy.array(scanM(:).');
        scanPy = np.array(scan3M);
        scanPy = scanPy.astype(py.numpy.float32);

        % get original shape of the scan
        %origShape = py.numpy.array(size(scanM));
        %origShape = origShape.astype(py.numpy.int64);

        % reshape numpy array to original shape
        %scanPy = reshape(scanPy,origShape);

        % Get image from the array
        itkimg = sitk.GetImageFromArray(scanPy);

        % Set Image spacing
        zSpacing = paramS.VoxelSize_mm.val(3);
        xSpacing = paramS.VoxelSize_mm.val(1);
        ySpacing = paramS.VoxelSize_mm.val(2);
        itkimg.SetSpacing([zSpacing, ySpacing, xSpacing]);

        % calculate gradient
        logRecursiveFilt = sitk.LaplacianRecursiveGaussianImageFilter();

        logRecursiveFilt.SetNormalizeAcrossScale(true)

        sigmaVal = paramS.Sigma_mm.val;
        logRecursiveFilt.SetSigma(sigmaVal);

        % Execute the filter
        logImg = logRecursiveFilt.Execute(itkimg);

        % extract numpy array from resulting image
        npLogImg = sitk.GetArrayFromImage(logImg);

        % convert resulting numpy array to matlab array in required shape
        %dblLogResultM = double(py.array.array('d',py.numpy.nditer(npLogImg)));
        %logMatM = reshape(dblLogResultM,origScanSize);

        logMatM = double(npLogImg);

        filteredOutS.logImg3M = logMatM;

    case 'HistogramMatchingImageFilter'
        % paramS inputs needed:
        % numHistLevel (int), paramS.numMatchPts (int),
        % ThresholdAtMeanIntensityOn (bool),
        % refImgPath (char vector)

        % convert scan to numpy array and integer
        scan3M = permute(scan3M,[3,1,2]);
        scanPy = np.array(scan3M);
        scanPy = scanPy.astype(py.numpy.float32);

        % get original shape of the scan
        %origShape = py.numpy.array(size(scanM(:)'));
        %origShape = origShape.astype(py.numpy.int64);

        % reshape numpy array to original shape
        %scanPy = reshape(scanPy,origShape);

        % Get image from the array

        srcItkImg = sitk.GetImageFromArray(scanPy);

        %             srcItkImg = py.SimpleITK.ReadImage('E:\data\TumorAware_MR\nrrdScanFormat.nrrd');

        % get ref image
        maskImgPath = fullfile(getCERRPath,'ModelImplementationLibrary/SegmentationModels/referenceImages');
        if isfield(paramS,'refImg') && ~isempty(paramS.refImg.val)
            [refPath,~,~,] = fileparts(paramS.refImg.val);
            if isempty(refPath)
                maskImgPath = fullfile(maskImgPath,paramS.refImg.val);
            else
                maskImgPath = paramS.refImg.val;
            end
            refItkImg = sitk.ReadImage(maskImgPath);
            refItkImg = sitk.Cast(refItkImg,py.SimpleITK.sitkFloat32);

            %Adjust to RTOG-compliant orientation
            %maskPy = sitk.GetArrayFromImage(refItkImg);
            %mask3M = single(maskPy);
            %mask3M = flip(flip(permute(mask3M,[2,3,1]),1),3); % flip y and z dims to convert from RAS to CERR coordinate system
            %mask3M = permute(mask3M,[3,1,2]);
            %scanRefPy = np.array(mask3M);
            %scanRefPy = scanRefPy.astype(py.numpy.float32);
            %refItkImg1 = sitk.GetImageFromArray(scanRefPy);

            %mask3M = flip(flip(mask3M,1),2);
            %mask3M = flip(mask3M,3);
            %maskPy = np.array(mask3M);
            %refItkImg = py.extra.GetImageFromArray(maskPy);
        elseif isfield(paramS,'refImgMat') && ~isempty(paramS.refImgMat.val)
            maskPy = np.array(paramS.refImgMat.val);
            maskPy = maskPy.astype(py.numpy.float32);
            % Get image from the array
            refItkImg = sitk.GetImageFromArray(maskPy);
        else
            error('Reference image not specified for histogram matching')
        end

        %refNumElems = int64(py.numpy.prod(refScanPy.shape));
        %refShape = py.numpy.array([1,refNumElems]);
        %refScanPy = refScanPy.astype(py.numpy.int64);
        % reshape numpy array to original shape
        %refScanPy = reshape(refScanPy,refShape);

        %refItkImg = py.SimpleITK.reshape(refItkImg, scanPy);
        % execute Histogram Matching
        matcher = sitk.HistogramMatchingImageFilter();
        matcher.SetNumberOfHistogramLevels(uint32(paramS.numHistLevel.val));
        matcher.SetNumberOfMatchPoints(uint32(paramS.numMatchPts.val));
        matcher.SetThresholdAtMeanIntensity(false)
        if paramS.thresholdAtMeanIntensityOn.val
            matcher.SetThresholdAtMeanIntensity(true);
        end
        matchedImg = matcher.Execute(srcItkImg,refItkImg);

        % extract numpy array from resulting image
        npHistImg = sitk.GetArrayFromImage(matchedImg);

        % convert resulting numpy array to matlab array in required shape
        %dblHistResultM = double(py.array.array('d',py.numpy.nditer(npHistImg)));
        %histMatM = reshape(dblHistResultM,[origScanSize(1),origScanSize(2),origScanSize(3)]);
        histMatM = double(npHistImg);

        filteredOutS.histMatchedImage = permute(histMatM,[2,3,1]);

        %             %visualize
        %             size(histMatM)
        %             figure, imagesc(histMatM(:,:,slc,1))
        %             figure, imagesc(histMatM(:,:,slc,2))
        %             figure, imagesc(histMatM(:,:,slc,3))

    case 'N4BiasFieldCorrectionImageFilter'

        % Convert to sitk image
        scan3M = permute(scan3M,[3,1,2]);
        scanPy = np.array(scan3M);
        scanPy = scanPy.astype(py.numpy.float32);
        srcItkImg = sitk.GetImageFromArray(scanPy);

        % Get mask
        if isfield(paramS,'maskImg') && ~isempty(paramS.maskImg.val)
            maskImgPath = fullfile(paramS.maskImg.val);
            maskImg = sitk.ReadImage(maskImgPath);
            maskImg = sitk.Cast(maskImg,py.SimpleITK.sitkFloat32);
            %Adjust to RTOG-compliant orientation
            maskPy = sitk.GetArrayFromImage(maskImg);
            mask3M = single(maskPy);
            mask3M = permute(mask3M,[2,3,1]);
            mask3M = flip(flip(mask3M,1),2);
            mask3M = flip(mask3M,3);
            maskPy = np.array(mask3M);
            maskImg = sitk.GetImageFromArray(maskPy);
        else
            %Use defaults
            %maskImg = py.SimpleITK.OtsuThreshold(srcItkImg,...
            %    uint8(0),uint8(1),uint32(200));
            maskImg = [];
        end

        corrector = sitk.N4BiasFieldCorrectionImageFilter();

        % Get optional parameters
        if isfield(paramS,'shrinkFactor') && paramS.shrinkFactor.val > 1
            shrink = paramS.shrinkFactor.val;
            srcItkImg = sitk.Shrink(srcItkImg,...
                shrink*srcItkImg.GetDimension());
            maskImg = sitk.Shrink(maskImg,...
                shrink*srcItkImg.GetDimension());
        end

        if isfield(paramS,'numFittingLevels')
            numFitLevels = paramS.numFittingLevels.val;
        else
            numFitLevels = 4;
        end

        if isfield(paramS,'numIterations')
            numIterationsV = paramS.numIterations.val; %Vector of values per fit level
            maxIterations = int32(numIterationsV);
            corrector.SetMaximumNumberOfIterations(maxIterations);
        end

        % Apply bias correction
        if ~isempty(maskImg)
            outImg = corrector.Execute(srcItkImg, maskImg);
        else
            outImg = corrector.Execute(srcItkImg);
        end

        logBiasFieldImg = corrector.GetLogBiasFieldAsImage(srcItkImg);

        % Convert to matlab array
        npLogBiasFieldImg = sitk.GetArrayFromImage(logBiasFieldImg);
        logBiasFieldImg3M = double(npLogBiasFieldImg);

        correctedImageFullResolution = scan3M ./ exp(logBiasFieldImg3M);

        filteredOutS.biasCorrectedImage = permute(correctedImageFullResolution,[2,3,1]);

        %if isfield(paramS,'shrinkFactor') && paramS.shrinkFactor.val > 1
        %  npCorrectedImg = py.extra.GetArrayFromImage(outImg);
        %  correctedImgShrunk3M = double(npCorrectedImg);
        %end

    case 'N4BiasAndHistogramCorrectionImageFilter'
        filteredOutS = sitkWrapper(sitkLibPath, scan3M, ...
            'N4BiasFieldCorrectionImageFilter', paramS);
        correctedScan3M = filteredOutS.biasCorrectedImage;
        filteredOutS = sitkWrapper(sitkLibPath, correctedScan3M, ...
            'HistogramMatchingImageFilter', paramS);

    otherwise

        msgStr = [filterType,' not defined. Add it to sitkWrapper.m'];
        error(msgStr)
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function score = grade(board,goal,solutionVector)
modSolutionVector = validateSolutionVector(solutionVector,board);
boardToScore = runSolutionVector(modSolutionVector, board, goal, false);
score = sum(boardToScore(:)) + sum(modSolutionVector);

author: Sable

repo: mcbench-benchmarks


text:
function a = wilk04 ( )

%*****************************************************************************80
%
%% WILK04 returns the WILK04 matrix.
%
%  Formula:
%
%    0.9143D-04  0.0         0.0         0.0
%    0.8762      0.7156D-04  0.0         0.0
%    0.7943      0.8143      0.9504D-04  0.0
%    0.8017      0.6123      0.7165      0.7123D-04
%
%  Properties:
%
%    A is lower triangular.
%
%    LAMBDA = ( 0.9143D-04, 0.7156D-04, 0.9504D-04, 0.7123D-04 ).
%
%  Discussion:
%
%    Since the matrix is already in lower triangular form, errors can
%    occur only in the backsubstitution.  However, even a double
%    precision calculation will show a significant degradation in the
%    solution.  It is also instructive to compare the actual error in
%    the solution to the residual error, A*x-b.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Wilkinson,
%    Rounding Errors in Algebraic Processes,
%    Prentice Hall, 1963, page 105.
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%

%
%  Note that the matrix entries are listed by row.
%
  a(1:4,1:4) = [ 0.9143E-04, 0.8762,     0.7943,     0.8017; ...
                 0.0000,     0.7156E-04, 0.8143,     0.6123; ...
                 0.0000,     0.0000,     0.9504E-04, 0.7165; ...
                 0.0000,     0.0000,     0.0000,     0.7123E-04 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%SETLABLIST Set names of classes or targets
%
%    A = SETLABLIST(A,LABLIST)
%
% LABLIST should be a column vector of C elements (integers, characters,
% or strings), in which C is the number of classes or targets of A.
% In case of multiple label lists this resets the current label list.
%
%    A = SETLABLIST(A)
%
% Removes entries in the lablist of A to which no objects are assigned,
% i.e. it remove empty classes. This command also removes duplicates in the
% lablist. An example of the merge of two classes in a 3-class dataset X
% with class names 'A', 'B' and 'C' can be realized by
%   X = setlablist(X,char('A','B','A'));
%   X = setlablist(X);
%
% SEE ALSO (PRTools Guide)
% MULTI_LABELING

author: marianux

repo: ecg-kit

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


text:
 function ob = Gtranslate(mask, varargin)
%function ob = Gtranslate(mask, options)
%|
%| Construct Gtranslate object for image registration.
%| (Internally stores image-sized arrays so not very memory efficient.)
%|
%| See Gtranslate_test() below for example usage.
%|
%| in
%|	mask	size(image)	logical array of object support.
%|
%| options
%|	'shift'	[ndim]		shift amount (default 0)
%|	'type'	char		translation method:
%|				'circshift' (integer shifts only) (default)
%|				'fft' (terrible for non-integer shifts)
%|				'interpn,linear,circ' linear interpolation
%|				with circulant boundary conditions
%|
%|	these methods do not have a properly matched adjoint (todo):
%|				'bspline3'
%|				'interpn,linear' linear interpolation
%|
%| out
%|	ob [nd np]	np = sum(mask(:)), so it is already "masked"
%|			nd = np for 'conv,same' type
%|
%| Copyright 2010-3-22, Jeff Fessler, University of Michigan
%|
%| 2012-8-1 'interpn,linear,circ' option added by F. Zhao and M. Muckley

if nargin == 1 && streq(mask, 'test'), Gtranslate_test, return, end
if nargin < 1, ir_usage, end

arg.mask = mask;

% option defaults
arg.shift = 0; % identity matrix
arg.type = 'circshift';
arg.class = 'fatrix2';

% options specified by name/value pairs
arg = vararg_pair(arg, varargin);

arg.shift = arg.shift(:)'; % row vector

arg.ndim = ndims(mask);
if arg.ndim == 2 && size(mask,2) == 1
	arg.ndim = 1;
end
if numel(arg.shift) ~= arg.ndim
	fail 'ndim mismatch'
end

switch arg.type

case 'bspline3'
	if ndims(mask) ~= 2, fail 'only 2d done', end

	ig = image_geom('nx', size(mask,1), 'dx', 1, ...
			'ny', size(mask,2), 'dy', 1);
	kg = knot_geom( 'nx', 1, 'mx', 8*ig.nx, 'offset_x', ig.nx/2, ...
			'ny', 1, 'my', 8*ig.ny, 'offset_y', ig.ny/2);
	Bx = makeB(ig, kg);
	By = Bx;
	alphax = zeros(kg.dim, 'single');
	alphay = zeros(kg.dim, 'single');
	alphax(1) = -2.25 * arg.shift(1); % trick: empirical
	if arg.ndim > 1
		alphay(1) = -2.25 * arg.shift(2);
	else
		alphay(1) = 0; % 1d
	end
	arg.W = makeW({Bx, By}, {alphax, alphay});

	arg.fun_forw = @(arg, x) ...
		reshape(arg.W * BsplVal2CoMirr(single(x)), size(x));
	arg.fun_back = @(arg, y) arg.mask .* ...
		reshape(arg.W' * y, size(y)); % todo: wrong
%	arg.fun_back = @(arg, y) BsplCo2ValTranMirr(arg.W' * y);
%	arg.fun_back = @(arg, y) fail 'bspline3 adjoint not done';
	warn 'bspline3 adjoint not done to match!';

case 'circshift'
	if any(round(arg.shift) ~= arg.shift) % need to use fft?
		fail 'circshift needs integer shifts'
	end
	arg.fun_forw = @(arg, x) circshift(x, arg.shift);
	arg.fun_back = @(arg, y) arg.mask .* circshift(y, -arg.shift);

case 'fft'
	Nd = size(mask);
	for id = 1:arg.ndim
		N = Nd(id);
		kk = single([0:N-1]);
		phase{id} = exp(-2i * pi / N * kk * arg.shift(id));
	end
	tmp = ndgrid_jf('mat', phase);
	tmp = prod(tmp, 1 + arg.ndim);
	arg.phase = single(tmp);

%	arg.fun_forw = @(x) ifftn(fftn(x) .* phase);
	arg.fun_forw = @(arg, x) translate_fft(x, arg.phase);
	arg.fun_back = @(arg, y) arg.mask .* translate_fft(y, conj(arg.phase));

case 'interpn,linear'
	Nd = size(mask);
	for id = 1:arg.ndim
		N = Nd(id);
		xf{id} = single([1:N]) - arg.shift(id);
		xb{id} = single([1:N]) + arg.shift(id);
	end
	if ndims(mask) == 2 && size(mask,2) == 1 % 1d
		arg.xf = col(xf{1});
		arg.xb = col(xb{1});
		arg.fun_forw = @(arg, x) interp1(x, arg.xf, 'linear', 0);
		arg.fun_back = @(arg, y) interp1(y, arg.xb, 'linear', 0);
	else
		arg.xf = ndgrid_jf('cell', xf);
		arg.xb = ndgrid_jf('cell', xb);
		arg.fun_forw = @(arg, x) interpn(x, arg.xf{:}, 'linear', 0);
		arg.fun_back = @(arg, y) arg.mask .* ...
			interpn(y, arg.xb{:}, 'linear', 0); % todo: wrong
	end
	warn 'interpn adjoint not done to match!';
%	arg.fun_back = @(arg, y) fail 'interp adjoint not done';

case 'interpn,linear,circ'
	Nd = size(mask) + 2 * ceil(abs(arg.shift));
	for id = 1:arg.ndim
		N = Nd(id);
		xf{id} = single([1:N]) - arg.shift(id);
		xb{id} = single([1:N]) + arg.shift(id);
	end
	arg.xf = ndgrid_jf('cell', xf);
	arg.xb = ndgrid_jf('cell', xb);
	arg.fun_forw = @(arg, x) translate_interp_circ_forw(arg, x);
	arg.fun_back = @(arg, y) arg.mask .* translate_interp_circ_back(arg, y);

otherwise
	fail('type "%s" unknown', arg.type)
end

% build object
switch arg.class

case 'Fatrix'
	arg.np = sum(mask(:));
	arg.nd = prod(size(mask));
	dim = [arg.nd arg.np];
	ob = Fatrix(dim, arg, 'caller', 'Gtranslate', ...
		'forw', @Gtranslate_forw_Fatrix, ...
		'back', @Gtranslate_back_Fatrix);

case 'fatrix2'
	idim = size(mask);
	if numel(idim) == 2 && idim(2) == 1
		idim = idim(1); % 1d
	end
	ob = fatrix2('mask', mask, 'arg', arg, ...
		'idim', idim, 'odim', idim, ...
		'forw', arg.fun_forw, 'back', arg.fun_back);

otherwise
	fail('class "%s" unknown', arg.class)
end

% Gtranslate_forw_Fatrix()
% y = A * x
function y = Gtranslate_forw_Fatrix(arg, x)

[x ei] = embed_in(x, arg.mask, arg.np);

if ndims(x) > arg.ndim
	y = zeros(arg.nd, 1, class(x));
	for ll=1:size(x, ndims(x))
		y(:,ll) = col(arg.fun_forw(arg, stackpick(x, ll)));
	end
	if ndims(arg.mask) > 2 || size(arg.mask, ndims(arg.mask)) > 1 % > 1d
		y = reshapee(y, size(arg.mask), []);
	end
else
	y = arg.fun_forw(arg, x);
end

y = ei.shape(y);

% Gtranslate_back_Fatrix()
% x = A' * y (adjoint)
function x = Gtranslate_back_Fatrix(arg, y)

[y eo] = embed_out(y, size(arg.mask));

if ndims(y) > ndims(arg.mask)
	x = zeros(arg.nd, 1, class(y)); % trick!
	for ll=1:size(y, ndims(y))
		x(:,ll) = col(arg.fun_back(arg, stackpick(y, ll)));
	end
	x = reshapee(x, size(arg.mask), []);
else
	x = arg.fun_back(arg, y);
end

x = eo.shape(x, arg.mask, arg.np);

% translate_fft()
function y = translate_fft(x, phase)
y = ifftn(fftn(x) .* phase);
if isreal(x)
	y = real(y);
end

% translate_interp_circ_forw()
function y = translate_interp_circ_forw(arg, x)
padsize = ceil(abs(arg.shift));
%x1 = padarray(x, padsize, 'circular', 'both'); % needs image toolbox
if numel(arg.shift) == 1
	npad = 2*padsize + size(x,1);
else
	npad = 2*padsize + size(x);
end
x1 = ir_pad_into_center(x, npad, 'circ', 1);
if numel(arg.shift) == 1
	y = interp1(x1, arg.xf{1}, 'linear', 0);
else
	y = interpn(x1, arg.xf{:}, 'linear', 0);
end
%msk = logical(padarray(true(size(x)), padsize)); % image toolbox
msk = ir_pad_into_center(true(size(x)), npad);
y = reshape(y(msk), size(x));

% translate_interp_circ_back()
function x = translate_interp_circ_back(arg, y)
padsize = ceil(abs(arg.shift));
%y1 = padarray(y, padsize, 'circular', 'both'); % needs image toolbox
if numel(arg.shift) == 1
	npad = 2*padsize + size(y,1);
else
	npad = 2*padsize + size(y);
end
y1 = ir_pad_into_center(y, npad, 'circ', 1);
if numel(arg.shift) == 1
	x = interp1(y1, arg.xb{1}, 'linear', 0);
else
	x = interpn(y1, arg.xb{:}, 'linear', 0);
end
%msk = logical(padarray(true(size(y)), padsize)); % image toolbox
msk = ir_pad_into_center(true(size(y)), npad);
x = reshape(x(msk), size(y));

% Gtranslate_test()
function Gtranslate_test

rng(0)
classes = {'fatrix2', 'Fatrix'};
for ic=1:2
for id=1:2
	if id == 1
		mask = true(8,1); ishift = 2; fshift = [2.3]; % 1d
	else
		mask = true(10,8); ishift = [7 3]; fshift = [2.3 1.7]; % 2d
	end
	mask(1) = false;
	args = {mask, 'class', classes{ic}};

	x = mask .* rand(size(mask));

	btypes = {'bspline3', 'interpn,linear'}; % non-matched adjoint cases
	ftypes = {'fft', 'interpn,linear,circ'}; % arbitrary shifts with adjoint
	itypes = {'circshift'}; % integer shifts (orthonormal operation)

	for it = 1:numel(btypes)
		btype = btypes{it};
		A = Gtranslate(args{:}, 'type', btype, 'shift', fshift);
		if 0 % todo
			fatrix2_tests(A, 'complex', 0, 'halt', 0, ...
			'check1', false, 'full', false) % because of bad adjoint
		end

		if 0
			a1 = full(A)
			a2 = full(A')'
		end
%		test_adjoint(A, 'complex', 1); % todo: fails due to bad adjoint
	end

	for it = 1:numel(ftypes)
		ftype = ftypes{it};
		A = Gtranslate(args{:}, 'type', ftype, 'shift', fshift);
		fatrix2_tests(A, 'complex', 1)
		test_adjoint(A, 'complex', 1);
	end

	for it = 1:numel(itypes)
		itype = itypes{it};
		A = Gtranslate(args{:}, 'type', itype, 'shift', ishift);

		if id > 1
			jf_equal(x, A' * (A * x)) % orthonormal
		end
		jf_equal(x(mask), A' * (A * x(mask))) % orthonormal
		fatrix2_tests(A, 'complex', 1)
		test_adjoint(A, 'complex', 1);
	end
end
end

if 0
	tmp = 0 * mask; tmp(ceil(end/2), ceil(end/2)) = 1;
	im plc 2 3
	im(1, tmp)
	tmp = A * tmp;
	im(2, abs(tmp), 'abs'), cbar
	im(3, real(tmp), 'real'), cbar
	im(4, imag(tmp), 'imag'), cbar
return
end

% try it for image registration
if 1
%	f.type = 'fft'; % bad results!
%	f.type = 'bspline3';
%	f.type = 'interpn,linear';
	f.type = 'interpn,linear,circ';

	f1 = ellipse_im(128, [], 'oversample', 2);

	im plc 2 2
	im(1, f1), cbar

	mask = true(size(f1));
	A = Gtranslate(mask, 'type', f.type, 'shift', fshift);
	y1 = A * f1;
	im(2, y1), cbar

	if 0 % examine integer shifts
		Ac = Gtranslate(mask, 'type', 'circshift', 'shift', ishift);
		yc = Ac * f1;
		Ab = Gtranslate(mask, 'type', f.type, 'shift', ishift);
		yb = Ab * f1;
		yb = reshape(yb, size(mask));
		im clf, im_toggle(yb, yc, [0 255])
	return
	end

	if 0 % examine shifts
		Af = Gtranslate(mask, 'type', f.type, 'shift', fshift);
		yf = Af * f1;
		Ab = Gtranslate(mask, 'type', 'bspline3', 'shift', fshift);
		yb = Ab * f1;
		yb = reshape(yb, size(mask));
		im clf, im_toggle(yf, yb, [0 255])
	return
	end

	shifts = linspace(0, 4, 21);
	cost = zeros(size(shifts));
	for is=1:length(shifts)
		tmp = [shifts(is) fshift(2)]; % in dim1 only
		A2 = Gtranslate(mask, 'type', f.type, 'shift', tmp);
		y2 = A2 * f1;
		cost(is) = norm(y2(:) - y1(:))^2;
	end

	if im
		im subplot 3
		plot(shifts, cost, '-o', fshift(1), 0, 'rx')
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Clear Matlab command window.
clc;

% Example #1: Basic usage.
% Note: Filters are computed for each run.
disp('Example #1: Basic Usage');
theta = [0:15:360];
for i = [1:length(theta)]
   [J,H] = steerGauss([],theta(i),3,true);
   filters{i} = H;
   pause(0.1);
end
disp('   Press any key to continue.'); pause;

% Load "mandrill" test image.
I = imread('mandrill.jpg');

% Example #2: Using pre-computed filters.
disp('Example #2: Using pre-computed filters.');
for i = [1:length(filters)]
   [J,H] = steerGauss(I,filters{i},true);
   pause(0.1);
end


text:
%
% Author: Marius Drulea
% http://www.cv.utcluj.ro/optical-flow.html
%
% References
% M. Drulea and S. Nedevschi, "Total variation regularization of
% local-global optical flow," in Intelligent Transportation Systems (ITSC),
% 2011 14th International IEEE Conference on, 2011, pp. 318-323.
%
% Copyright (C) 2011 Technical University of Cluj-Napoca
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

function [u] = tv_min(u0, lambda)
%   The output u approximately minimizes the Rudin-Osher-Fatemi (ROF)
%   denoising model
%
%       Min  TV(u) + 1/(2*lambda)* || u - u0 ||^2_2,
%        u

u = zeros(size(u0));

% initialization of the dual variable
p1 = zeros(size(u0));
p2 = zeros(size(u0));

tau = 1/4;
max_iters = 100;

der_mask = [0 -1 1];
adjoint_der_mask = [-1 1 0];

for i=1:max_iters

    % the divergence
    div_p = imfilter(p1, adjoint_der_mask, 'replicate') + ...
        imfilter(p2, adjoint_der_mask', 'replicate');

    t = div_p - u0/lambda;

    % the derivatives
    tx = imfilter(t, der_mask, 'replicate');
    ty = imfilter(t, der_mask', 'replicate');

    denominator = 1 + tau * sqrt(tx.^2 + ty.^2);

    % update dual variable; gradient ascent and reprojection
    p1 = (p1 + tau*tx)./denominator;
    p2 = (p2 + tau*ty)./denominator;

    % update variable; gradient descent
    u = u0 - lambda*div_p;
end

end

author: thomas-koehler

repo: SupER


text:
function [err, out_fn] = Norm (InBrik, N_runs, run, run_norm, run_mask, foutid_log)
%
%   [err, norm_fn] = Norm(InBrik, N_runs, run)
%
%Purpose:
%
%   Normalize the original dataset so that the beta values in 3dDeconvolve would be automatically percent signal change.
%
%Input Parameters:
%
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%
%
%
%Key Terms:
%
%More Info :
%
%
%
%
%     Author : Gang Chen
%     Date : Mon Oct 20 11:26:38 EDT 2003
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'Norm';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

[InBrikPrefix, InBrikView, InBrikExt] =  AfniPrefix (InBrik);
rm_cmmd = 'rm -f ';

%Originally 3dClipLevel was used to generate a cliplevel for each run,
%then a step function chops off those voxel outside of the brain.
%But due to the concern of losing the sensitive amygdala and orbitofrontal cortex,
%it might be better to run 3dAutomask, which is roughly equivalent to
%3dClipLevel but would protect amygdala and orbitofrontal cortex regions
%from being masked.

%create command for executing 3dAutomask
if (run_mask.do == 1),
   automask_fn = sprintf('%s_msk', InBrikPrefix);
   automask = sprintf('3dAutomask -dilate %i -prefix %s %s%s', run_mask.N_dilate, automask_fn, InBrikPrefix, InBrikView);
   fprintf(2, '\nRunning: %s\n', automask);
   fprintf(foutid_log, '%s\n', automask);
   [s,w] = system(automask);
   if (s),
      fprintf(2,'Error during 3dAutomask: \n\n%s\n\n',w);
      while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
   end
	rm_cmmd = sprintf('%s %s%s%s', rm_cmmd, automask_fn, InBrikView, InBrikExt);

	%If no normalization, do masking
   if (run_norm == 0),
	   out_fn = sprintf('%s_mask', InBrikPrefix);
      mask_signal = sprintf('3dcalc -fscale -a  %s%s  \\\n', InBrikPrefix, InBrikView);
 	   mask_signal = sprintf('%s    -b %s%s \\\n', mask_signal, automask_fn, InBrikView);
	   mask_signal = sprintf('%s    -expr  "a*b" \\\n', mask_signal);
      mask_signal = sprintf('%s    -prefix %s', mask_signal, out_fn);
	   fprintf(2, '\nRunning: %s\n', mask_signal);
	   fprintf(foutid_log, '%s\n', mask_signal);
	   [s,w] = system(mask_signal);
      if (s),
         fprintf(2,'Error during 3dcalc: \n\n%s\n\n',w);
         while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
      end
   end
	rm_cmmd = sprintf('%s %s%s%s', rm_cmmd, automask_fn, InBrikView, InBrikExt);
end

if (run_norm == 1),   %for normalization
   out_fn = sprintf('%s_norm', InBrikPrefix);
   concat = sprintf('3dTcat -prefix %s ', out_fn);

   for (i=1:1:N_runs),

   %store each mean dataset in a temporary file, and delete it before leaving this function.
      mean_fn = sprintf('%s_mean%i', InBrikPrefix, i);
	   get_mean = sprintf('3dTstat -mean -prefix %s %s%s''[%i..%i]''', mean_fn, InBrikPrefix, InBrikView, run(i).first, run(i).last);
      fprintf(2, '\nRunning: %s\n', get_mean);
	   fprintf(foutid_log, '%s\n', get_mean);
	   [s, w] = system(get_mean);
	   if (s),
         fprintf(2,'Error at running 3dTstat: \n\n\n%s\n', w);
         while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
      end

	%normalizing
%	norm_signal = sprintf('3dcalc -datum float -a  %s%s''[%i..%i]''  \\\n', InBrikPrefix, InBrikView, run(i).first, run(i).last);
	   norm_signal = sprintf('3dcalc -fscale -a  %s%s''[%i..%i]''  \\\n', InBrikPrefix, InBrikView, run(i).first, run(i).last);
   	norm_signal = sprintf('%s    -b %s%s \\\n', norm_signal, mean_fn, InBrikView);
		if (run_mask.do == 1),
		   norm_signal = sprintf('%s    -c %s%s \\\n', norm_signal, automask_fn,InBrikView);
		%	norm_signal = sprintf('%s    -expr  "(a/b*100)*step(b-clip_level)" \\\n', norm_signal);
   	   norm_signal = sprintf('%s    -expr  "(a/b*100)*c" \\\n', norm_signal);
		else norm_signal = sprintf('%s    -expr  "a/b*100" \\\n', norm_signal);
		end
	   norm_out_fn = sprintf('%s_norm%i', InBrikPrefix, i);
   	norm_signal = sprintf('%s    -prefix %s', norm_signal, norm_out_fn);
		fprintf(2, '\nRunning: %s\n', norm_signal);
	   fprintf(foutid_log, '%s\n', norm_signal);
	   [s,w] = system(norm_signal);
      if (s),
         fprintf(2,'Error during 3dcalc: \n\n%s\n\n',w);
         while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
      end

	%Concatenate the above files together
	   concat = sprintf('%s %s%s%s', concat, norm_out_fn, InBrikView, InBrikExt);

	%create remove command for deleting those individual percent signal change and mean file for each run.
	   rm_cmmd = sprintf('%s %s%s%s %s%s%s', rm_cmmd, norm_out_fn, InBrikView, InBrikExt, mean_fn, InBrikView, InBrikExt);

   end

%Concatenate all runs back together
   fprintf(2, '\nRunning: %s\n', concat);
   fprintf(foutid_log, '%s\n', concat);
   [s,w] = system(concat);
   if (s),
      fprintf(2,'Error during running 3dTcat: \n\n%s\n\n',w);
      while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
   end
end

%Don't forget to remove those temporary files above.
fprintf(2, '\nRemoving: %s\n', rm_cmmd);
fprintf(foutid_log, '%s\n', rm_cmmd);
[s,w] = system(rm_cmmd);
if (s),
    fprintf(2,'Error during rm: \n\n%s\n\n',w);
    while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
end

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function prob_test135 ( )

%*****************************************************************************80
%
%% TEST135 tests SECH_CDF, SECH_CDF_INV, SECH_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST135\n' );
  fprintf ( 1, '  For the Hyperbolic Secant PDF:\n' );
  fprintf ( 1, '  SECH_CDF evaluates the CDF.\n' );
  fprintf ( 1, '  SECH_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  SECH_PDF evaluates the PDF.\n' );

  a = 3.0;
  b = 2.0;

  check = sech_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST135 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =         %14f\n', a );
  fprintf ( 1, '  PDF parameter B =         %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = sech_sample ( a, b, seed );

    pdf = sech_pdf ( x, a, b );

    cdf = sech_cdf ( x, a, b );

    x2 = sech_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
function xy = spike_position(dd, wPCA, wTEMP, yc, xc)

wTPC = gpuArray(wPCA' * wTEMP);

dT = gpuArray(dd);
for j = 1:size(dd,3)
    dT(:,:,j) = dT(:,:,j) * wTPC;
end

[nspikes, nPC, nchan] = size(dT);
[~, imax] = max(max(dT.^2, [], 3), [], 2);
dBest = gpuArray.zeros(nspikes, nchan, 'single');
for j = 1:nPC
    iX = imax==j;
    dBest(iX, :) = dT(iX, j, :);
end

dBest = max(0, dBest);
dBest = dBest ./ sum(dBest,2);

ysp = dBest * yc;
xsp = dBest * xc;

xy = [xsp, ysp];

xy = gather(xy);
end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
clear all;clc;
%  close all;
outputDir = GetOutputDataDir;

%% Init load
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% params
isStimLockRanking = 0; % init
isMotorRanking = 0; % init

for caseflag = 5
    switch caseflag
        case 1 % seeds
            stimrange = [];
            M_clus_range = {[11,2],[12,1]};
            M_clus_name = {'motorseed2','ABN'};
            M_fishrange = {[1:3,5:18],[1:12,14:18]};
            n_reg = length(M_clus_name);

        case 2 % testing
            stimrange = [];
            M_clus_range = {[11,2],[11,3],[11,4]};
            M_clus_name = {'motorseed_v2','motorseed_LRmanual','motorseed_Fwmanual'};
            M_fishrange = {[5:13],[5:13],[5:13]};
            n_reg = length(M_clus_name);

        case 3 % conserved clusters
            stimrange = 1;
            M_clus_range = {[10,1]};
            M_clus_name = {'conservedclus_SLranked'};
            M_fishrange = {1:18};
            n_reg = length(M_clus_name);
            isStimLockRanking = 1;

        case 4 % all Auto0.7 clusters
            stimrange = [];
            M_clus_range = {[6,1]};
            M_clus_name = {'A0.7_def_stimlockrank'};
            M_fishrange = {1:18};
            n_reg = length(M_clus_name);
            isStimLockRanking = 1;

        case 5 % all Auto0.7 clusters
            stimrange = [];
            M_clus_range = {[6,1]};
            M_clus_name = {'A0.7_def_motorrank'};
            M_fishrange = {[1:3,5:18]};
            n_reg = length(M_clus_name);
            isStimLockRanking = 0;
            isMotorRanking = 1;

        case 6 % HBO stripes
            stimrange = 2;
            M_clus_range = {[10,4]};
            M_clus_name = {'HBO_anat'};
%             M_clus_name = {'HBO_OMR_stimlockrank'};
            M_fishrange = {[1:3,5:18]};%{8:18};%{1:18};%
            n_reg = length(M_clus_name);
            isStimLockRanking = 0;%1;
            isMotorRanking = 0;

            P_stimscore = zeros(18,4);

        case 7 % sparse network
            stimrange = 2;
            M_clus_range = {[12,2]};
            M_clus_name = {'sparse_nw'};
            %             M_clus_name = {'HBO_OMR_stimlockrank'};
            M_fishrange = {[1:8]};%{8:18};%{1:18};%
            n_reg = length(M_clus_name);

%                     case 6 % HBO stripes
%             stimrange = 1;
%             M_clus_range = {[10,4]};
%             M_clus_name = {'HBO_PT_stimlockrank'};
%             M_fishrange = {1:18};%{1:18};%{[1:3,5:18]};
%             n_reg = length(M_clus_name);
%             isStimLockRanking = 1;
%
%             P_stimscore = zeros(18,4);
%
%         case 3 % Looming
%             stimrange = 5;
%             M_stimmotorflag = [1,0]; % 1 for stim and 0 for motor
%             M_reg_name = {'Loom_LR','Loom_SwimLR'};
%             M_reg_range =  {[11,12],[1,3]};
%             M_fishrange = {[9:15,17:18],[9:15,17:18]};
%             n_reg = length(M_reg_name);
%
%         case 4 % Dark Flash (black-white)
%             stimrange = 3;
%             M_stimmotorflag = [1,0]; % 1 for stim and 0 for motor
%             M_reg_name = {'DF_BW','DF_SwimLR'};
%             M_reg_range =  {[1,4],[1,3]};
%             M_fishrange = {[12:15,17:18],[12:15,17:18]};
%             n_reg = length(M_reg_name);
    end

    M_fishrange_im = M_fishrange;

    %% Load fish
    range_fish = 1:18;
    IM_full = cell(n_reg,max(range_fish));
    %%
    for i_fish = range_fish
        %% load fish
        if ismember(i_fish,cell2mat(M_fishrange))
            ClusterIDs = [1,1];
            if isempty(stimrange)
                LoadSingleFishDefault(i_fish,hfig,ClusterIDs);
            else
                LoadSingleFishDefault(i_fish,hfig,ClusterIDs,stimrange);
            end
            %         [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);
        end

        %% load clusters
        for i_set = 1:n_reg
            range_set = M_fishrange{i_set};
            if ismember(i_fish,range_set)
                ClusterIDs = M_clus_range{i_set};%[12,1];% GetClusterIDs('all');
                [cIX,gIX] = LoadCluster_Direct(i_fish,ClusterIDs(1),ClusterIDs(2));

                % override:
                if caseflag==7
                   gIX = ones(size(cIX));
                end
%                 if caseflag==6
%                     [cIX1,gIX1] = LoadCluster_Direct(i_fish,10,2);
%                     [cIX2,gIX2] = LoadCluster_Direct(i_fish,10,3);
%                     cIX = [cIX1;cIX2];
%                     gIX = [gIX1;5-gIX2];

%                     absIX = getappdata(hfig,'absIX');
%                     name = 'HBO_4stripes';
%                     clusgroupID = 10;
%                     clusIDoverride = 4;
%                     SaveCluster_Direct(cIX,gIX,absIX,i_fish,name,clusgroupID,clusIDoverride);
%                 end

                if isempty(cIX)
                    M_fishrange_im{i_set} = setdiff(M_fishrange_im{i_set} ,i_fish);
                    continue;
                end

                %% [option: ranking]
                if caseflag>=3
                    if isStimLockRanking
                        M = UpdateIndices_Manual(hfig,cIX,gIX);
                        C = FindClustermeans(gIX,M);
                        [~,~,H] = GetTrialAvrLongTrace(hfig,C);
                        [gIX,rankscore] = SortGroupIXbyScore(H,gIX);

                        if caseflag==6
                            P_stimscore(i_fish,:) = H;
                        end
                    elseif isMotorRanking
                        M = UpdateIndices_Manual(hfig,cIX,gIX);
                        C = FindClustermeans(gIX,M);
                        numU = max(gIX);
                       [gIX,rankscore] = RankByMotorReg_Direct(hfig,gIX,numU,C,1);
                    end
                end

                %% make figure

                if caseflag==3 || caseflag==4 || caseflag==6
                    setappdata(hfig,'clrmap_name','hsv_new');
%                 elseif caseflag==6
%                     setappdata(hfig,'clrmap_name','jet');
                else
                    setappdata(hfig,'clrmap_name','hsv_old');
                end

                I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
                [h,im_full] = DrawCellsOnAnat(I);

                %% save figure
                close(h);
                IM_full{i_set,i_fish} = im_full;
            end
        end
    end

%     if caseflag==6
%         % save VAR
%         SaveVARtoMat(hfig);
%     end

    %% save as tiff stack
    for i_set = 1:n_reg
        range_im = M_fishrange_im{i_set};
        tiffdir = fullfile(outputDir,[M_clus_name{i_set},'_allfish.tiff']);
        IM = IM_full(i_set,range_im);

        SaveImToTiffStack(IM,tiffdir);
    end

    %% [for later] plot from tiff stack
    isPlotfromtiffstack = 0;

    if isPlotfromtiffstack
        IM_full = cell(n_reg,18);
        for i_set = 1:n_reg
            %% get tiff-stack path
            tiffdir = fullfile(outputDir,[M_clus_name{i_set},'_allfish.tiff']);

            %% load

            for i_fish = 1:18
                im = double(imread(tiffdir,i_fish))./255;
                IM_full{i_set,i_fish} = im(317:1236,1:621,:);
            end
        end
    end

    %% save average tiff image
    for i_set = 1:n_reg;
        range_im = M_fishrange_im{i_set};%[1:3,5:7];%[1:3,5:18];
        cellarray = IM_full(i_set,range_im);

        % adjust params for visualization
        k_scale = 0.5;
        k_contrast = 1;% 1.2;

        [h_anat,im_avr] = AverageAnatPlot(cellarray,k_contrast,k_scale);

        tiffdir = fullfile(outputDir,[M_clus_name{i_set},'_avr.tiff']);
        imwrite(im_avr, tiffdir, 'compression','none','writemode','overwrite');
    end

end

author: xiuyechen

repo: FishExplorer


text:
function [Y_F,U] = SimulateFactorsResiduals(Corr_Y_F, volU, J)
%[Y_F,U] = SimulateFactorsResiduals(Corr_Y_F, volU, J)
%  Generates simulations of unit normal factors Y_F and residuals U. The
%  first two sample joint moments of Y_F are forced to match the given
%  inputs and assumptions. Sample covariance of each residual U with each
%  factor is enforced to be zero, but the sample covariances between
%  elements of U are not forced to be zero (for complexity reasons).
%  U is assumed to be zero-mean, uncorrelated and jointly normal.
%
% IMPORTANT: The order of computations in the following expressions is
% important for efficiency. Please use the same order, or have a good reason
% to change it!
% Code by S. Gollamudi. This version December 2009.

K = size(Corr_Y_F,1);
N = length(volU);
if size(volU,2)==1, volU = volU'; end

% Simulate Y_F with sample mean and covariance exactly matching the
% desired values
Y_F = MvnRndMatchCrossCov(Corr_Y_F, zeros(0,K), zeros(J,0));

% Simulate residuals that are orthogonal to columns of F
U = randn(J/2,N);
U = [U
     -U];
U = U - (Y_F/Corr_Y_F)*((Y_F'*U)/J);

% Scale residuals to have the correct sample variance
norm_factors = volU./sqrt(mean(U.^2,1));
U = U .* repmat(norm_factors,J,1);


text:
function a = full(a)
%FULL         Convert slope to full
%

% written  04/06/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  a.r = full(a.r);
  a.s = full(a.s);

repo: OpenMAS


text:
function obj = subsasgn(obj,subs,dat)
% Overloaded subsasgn function for file_array objects.
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

if isempty(subs)
    return;
end;
if ~strcmp(subs(1).type,'()'),
    if strcmp(subs(1).type,'.'),
        %error('Attempt to reference field of non-structure array.');
        if numel(struct(obj))~=1,
            error('Can only change the fields of simple file_array objects.');
        end;
        switch(subs(1).subs)
        case 'fname',      obj = asgn(obj,@fname,     subs(2:end),dat); %fname(obj,dat);
        case 'dtype',      obj = asgn(obj,@dtype,     subs(2:end),dat); %dtype(obj,dat);
        case 'offset',     obj = asgn(obj,@offset,    subs(2:end),dat); %offset(obj,dat);
        case 'dim',        obj = asgn(obj,@dim,       subs(2:end),dat); %obj = dim(obj,dat);
        case 'scl_slope',  obj = asgn(obj,@scl_slope, subs(2:end),dat); %scl_slope(obj,dat);
        case 'scl_inter',  obj = asgn(obj,@scl_inter, subs(2:end),dat); %scl_inter(obj,dat);
        case 'permission', obj = asgn(obj,@permission,subs(2:end),dat); %permission(obj,dat);
        otherwise, error(['Reference to non-existent field "' subs.subs '".']);
        end;
        return;
    end;
    if strcmp(subs(1).type,'{}'), error('Cell contents reference from a non-cell array object.'); end;
end;

if numel(subs)~=1, error('Expression too complicated');end;
dm   = size(obj);
sobj = struct(obj);

if length(subs.subs) < length(dm),
    l   = length(subs.subs);
    dm  = [dm(1:(l-1)) prod(dm(l:end))];
    if numel(sobj) ~= 1,
        error('Can only reshape simple file_array objects.');
    end;
    if numel(sobj.scl_slope)>1 || numel(sobj.scl_inter)>1,
        error('Can not reshape file_array objects with multiple slopes and intercepts.');
    end;
end;

dm   = [dm ones(1,16)];
do   = ones(1,16);
args = {};
for i=1:length(subs.subs),
    if ischar(subs.subs{i}),
        if ~strcmp(subs.subs{i},':'), error('This shouldn''t happen....'); end;
        args{i} = int32(1:dm(i));
    else
        args{i} = int32(subs.subs{i});
    end;
    do(i) = length(args{i});
end;
for j=1:length(sobj),
    if strcmp(sobj(j).permission,'ro'),
        error('Array is read-only.');
    end
end

if length(sobj)==1
    sobj.dim = dm;
    if numel(dat)~=1,
        subfun(sobj,double(dat),args{:});
    else
        dat1 = double(dat) + zeros(do);
        subfun(sobj,dat1,args{:});
    end;
else
    for j=1:length(sobj),
        ps  = [sobj(j).pos ones(1,length(args))];
        dm  = [sobj(j).dim ones(1,length(args))];
        siz = ones(1,16);
        for i=1:length(args),
            msk      = args{i}>=ps(i) & args{i}<(ps(i)+dm(i));
            args2{i} = find(msk);
            args3{i} = int32(double(args{i}(msk))-ps(i)+1);
            siz(i)   = numel(args2{i});
        end;
        if numel(dat)~=1,
            dat1 = double(subsref(dat,struct('type','()','subs',{args2})));
        else
            dat1 = double(dat) + zeros(siz);
        end;
        subfun(sobj(j),dat1,args3{:});
    end
end
return

function sobj = subfun(sobj,dat,varargin)
va = varargin;

dt  = datatypes;
ind = find(cat(1,dt.code)==sobj.dtype);
if isempty(ind), error('Unknown datatype'); end;
if dt(ind).isint, dat(~isfinite(dat)) = 0; end;

if ~isempty(sobj.scl_inter),
    inter = sobj.scl_inter;
    if numel(inter)>1,
        inter = resize_scales(inter,sobj.dim,varargin);
    end;
    dat = double(dat) - inter;
end;

if ~isempty(sobj.scl_slope),
    slope = sobj.scl_slope;
    if numel(slope)>1,
        slope = resize_scales(slope,sobj.dim,varargin);
        dat   = double(dat)./slope;
    else
        dat   = double(dat)/slope;
    end;
end;

if dt(ind).isint, dat = round(dat); end;

% Avoid warning messages in R14 SP3
wrn = warning;
warning('off');
dat   = feval(dt(ind).conv,dat);
warning(wrn);

nelem = dt(ind).nelem;
if nelem==1,
    mat2file(sobj,dat,va{:});
elseif nelem==2,
    sobj1       = sobj;
    sobj1.dim   = [2 sobj.dim];
    sobj1.dtype = dt(find(strcmp(dt(ind).prec,{dt.prec}) & (cat(2,dt.nelem)==1))).code;
    dat         = reshape(dat,[1 size(dat)]);
    dat         = [real(dat) ; imag(dat)];
    mat2file(sobj1,dat,int32([1 2]),va{:});
else
    error('Inappropriate number of elements per voxel.');
end;
return

function obj = asgn(obj,fun,subs,dat)
if ~isempty(subs),
    tmp = feval(fun,obj);
    tmp = subsasgn(tmp,subs,dat);
    obj = feval(fun,obj,tmp);
else
    obj = feval(fun,obj,dat);
end;


text:
function[varargout]=ecconv(varargin)
%ECCONV  Converts between eccentricity measures.
%
%   Y=ECCONV(X,STR) converts the eccentricity measure X into eccentricity
%   measure Y as specified by the string STR.
%
%   STR is of the form 'in2out' where the names 'in' and 'out' may be
%   chosen from the following table.
%
%       STR        Name                    Definition
%       ----------------------------------------------------------------
%      'ecc'       Eccentricity            sgn(b)*sqrt(1-b^2/a^2)
%      'ell'       Ellipticity             b/a
%      'lin'       Linearity               sgn(b)*(a^2-b^2)/(a^2+b^2)
%      'circ'      Circularity             2ab/(a^2+b^2)
%      'nu'        Eccentricity angle      asin(b/sqrt(a^2 + b^2))
%      'rot'       Rotary ratio            sgn(b)*(a-abs(b))/(a+abs(b))
%
%   where "a" and "b" are the semi-major and semi-minor axes respective.
%   For example, ECC=ECCONV(L,'lin2ecc') converts the ellipse linearity L
%   into eccentricity ECC.
%
%   All quantities are positive for mathematically positive rotation
%   and negative for mathematically negative rotation.  Note the linearity
%   is undefined for the case of purely circular rotation a=abs(b);
%
%   For further details, see Lilly and Gascard (2006).
%   ____________________________________________________________________
%
%   Cell array input/output
%
%   If ECCONV is given cell array input, it returns cell array output.
%
%   Thus X may be a cell arrays of numerical arrays, and the output Y will
%   also be a cell array having an identical size as X.
%   ____________________________________________________________________
%
%   See also ELLCONV, ELLDIFF, TRANSCONV.
%
%   Usage:   lin=ecconv(b./a,'ell2lin');
%            ecc=ecconv(lin,'lin2ecc');
%
%   'ecconv --t' runs a test.
%   'ecconv --f' generates a sample figure.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2005--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
  ecconv_test,return
end
if strcmpi(varargin{1}, '--f')
  type makefigs_ecconv
  makefigs_ecconv;
  return
end

%if nargin==2

x=varargin{1};
str=varargin{2};
ecconv_checkstr(str);

if ~iscell(x)
   y=ecconv_one(x,str);
   y(isinf(x))=inf;
else
    for i=1:length(x)
         y{i,1}=ecconv_one(x{i},str);
         y{i}(isinf(x{i}))=inf;
    end
end

varargout{1}=y;

function[y]=ecconv_one(x,str)

sizex=size(x);
x=x(:);

str1=str(1:strfind(str,'2')-1);
str2=str(strfind(str,'2')+1:end);
if aresame(str1,'circ')&&~aresame(str2,'circ')
    lin=sign(x).*sqrt(1-x.^2);
    eval(['y=lin2' str2 '(x);'])
elseif ~aresame(str1,'circ')&&aresame(str2,'circ')
    eval(['lin=' str1 '2lin(x);'])
    y=sign(lin).*sqrt(1-lin.^2);
else
    eval(['y=' str '(x);'])
end
y=reshape(y,sizex);

% elseif nargin==3
%   x1=varargin{1};
%   x2=varargin{2};
%   str=varargin{3};
%   ecconv_checkstr(str);
%   eval(['[y1,y2]=' str '(x1,x2);'])
%   varargout{1}=y1;
%   varargout{2}=y2;
% end

function[]=ecconv_checkstr(str)
ii=strfind(str,'2');
if isempty(ii)
  error('STR should be of the form ''ecc2lin''.')
end
str1=str(1:ii-1);
str2=str(ii+1:end);

if ~aresame(str1,'nu') && ~aresame(str1,'ecc') && ~aresame(str1,'lin') && ~ ...
      aresame(str1,'rot') && ~aresame(str1,'ell') && ~aresame(str1,'circ')
  error(['Input parameter name ' str1 ' is not supported.'])
end
if ~aresame(str2,'nu') && ~aresame(str2,'ecc') && ~aresame(str2,'lin') && ~ ...
      aresame(str2,'rot') && ~aresame(str2,'ell') &&~aresame(str2,'circ')
  error(['Output parameter name ' str2 ' is not supported.'])
end

function[nu]=ell2nu(ell)
nu=atan(ell);
function[ell]=nu2ell(nu)
ell=tan(nu);
function[ecc]=ell2ecc(ell)
ecc=nu2ecc(ell2nu(ell));
function[ell]=ecc2ell(ecc)
ell=nu2ell(ecc2nu(ecc));
function[rot]=ell2rot(ell)
rot=nu2rot(ell2nu(ell));
function[ell]=rot2ell(rot)
ell=nu2ell(rot2nu(rot));
function[lin]=ell2lin(ell)
lin=nu2lin(ell2nu(ell));
function[ell]=lin2ell(lin)
ell=nu2ell(lin2nu(lin));

function[ecc]=nu2ecc(nu)
ecc=sign(nu).*sqrt(1-squared(tan(nu)));
vindexinto(ecc,1,find(nu==0),1);
vindexinto(ecc,0,find(abs(nu)==pi/4),1);

function[nu]=ecc2nu(ecc)
nu=sign(ecc).*atan(sqrt(1-squared(ecc)));
vindexinto(nu,pi/4,find(ecc==0),1);

function[lambda]=nu2lin(nu)
lambda=sign(nu).*cos(2*nu);
vindexinto(lambda,1,find(nu==0),1);
vindexinto(lambda,0,find(abs(nu)==pi/4),1);

function[nu]=lin2nu(lambda)
nu=sign(lambda).*atan(sqrt(frac(1-abs(lambda),1+abs(lambda))));
vindexinto(nu,pi/4,find(lambda==0),1);

function[alpha]=nu2rot(nu)
alpha=sign(nu).*frac(1-abs(tan(nu)),1+abs(tan(nu)));
vindexinto(alpha,1,find(nu==0),1);
vindexinto(alpha,0,find(abs(nu)==pi/4),1);

function[nu]=rot2nu(alpha)
nu=sign(alpha).*atan(frac(1-abs(alpha),1+abs(alpha)));
vindexinto(nu,pi/4,find(alpha==0),1);

function[ecc]=lin2ecc(lambda)
ecc=sign(lambda).*sqrt(frac(2*abs(lambda),1+abs(lambda)));
vindexinto(ecc,0,find(lambda==0),1);

function[lambda]=ecc2lin(ecc)
lambda=sign(ecc).*frac(squared(ecc),2-squared(ecc));
vindexinto(lambda,0,find(ecc==0),1);

function[alpha]=lin2rot(lambda)
alpha=nu2rot(lin2nu(lambda));
vindexinto(alpha,0,find(lambda==0),1);

function[lambda]=rot2lin(alpha)
lambda=nu2lin(rot2nu(alpha));
vindexinto(lambda,0,find(alpha==0),1);

function[alpha]=ecc2rot(ecc)
alpha=nu2rot(ecc2nu(ecc));
vindexinto(alpha,0,find(ecc==0),1);

function[ecc]=rot2ecc(alpha)
ecc=nu2ecc(rot2nu(alpha));
vindexinto(ecc,0,find(alpha==0),1);

function[x]=rot2rot(x)
function[x]=nu2nu(x)
function[x]=ecc2ecc(x)
function[x]=lin2lin(x)
function[x]=ell2ell(x)
function[x]=circ2circ(x)

%/********************************************************
function[]=ecconv_test
str{1}='nu';
str{2}='ecc';
str{3}='lin';
str{4}='rot';
str{5}='ell';

x=rand(100,1)*2-1;
x(end+1,1)=1;  %Add an exact one
x(end+1,1)=0;  %Add an exact zero

clear y z x2 bool

%Turn lambda into all others
for i=1:length(str)
  y(:,i)=ecconv(x,['lin2' str{i}]);
end

%Turn all others into everything else
for i=1:length(str)
  for j=1:length(str)
     z(:,i,j)=ecconv(y(:,i),[str{i} '2' str{j}]);
  end
end

%Turn everything back into lambda
for i=1:length(str)
  for j=1:length(str)
     x2(:,i,j)=ecconv(z(:,i,j),[str{j} '2lin']);
  end
end

for i=1:size(x2,2)
  for j=1:size(x2,2)
    bool(i,j)=aresame(squeeze(x2(:,i,j)),x,1e-14);
  end
end

reporttest('ECCONV', allall(bool))
%\********************************************************

% %Not currently used
% figure
% lambda=(0:.001:1)';
%
% ecc=ecconv(lambda,'lin2ecc');
% alpha=ecconv(lambda,'lin2rot');
% ecc1=sqrt(2)*sqrt(lambda);
% alpha1=lambda/2;
%
% plot(lambda,[ecc,alpha,ecc1,alpha1]);
% linestyle k 2k k-- 2k--
% axis([0 1 0 1]),axis square
% text(0.4,0.7,'\epsilon')
% text(0.7,0.5,'\alpha')
% title('Eccentricity measures')
% xlabel('Ellipse parameter \lambda')
% ylabel('Eccentricity \epsilon or rotary ratio \alpha')
% set(gcf,'paperposition', [2 2 3.5 3.5])
% xtick(.1),ytick(.1),fixlabels(-1)
% fontsize 14 14 14 14

% function[kappa,nu]=ab2kn(a,b)
% kappa=frac(1,sqrt(2))*sqrt(squared(a)+squared(b));
% nu=atan(frac(b,a));
%
% function[a,b]=kn2ab(kappa,nu)
% a=cos(nu).*kappa*sqrt(2);
% b=sin(nu).*kappa*sqrt(2);
%
% function[nu]=ab2nu(a,b)
% nu=atan(frac(b,a));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function writeGLTF(filename,V,F,varargin)
  % writeGLTF(filename,V,F)
  %
  % Inputs:
  %   filename  path to .gltf file
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of mesh triangle indices into rows of V
  %   Optional:
  %   'SkinningTransforms' followed by
  %     T  3 by 4 by #T by #A list of animated global skinning transformations
  %   'SkinningWeights' followed by
  %     W  #V by #T list of skinning weights
  %   'MorphTargets' followed by
  %     MV  #V by 3 by #M list of morph targets
  %   'MorphWeights' followed by
  %     MW  #frames by #M list of morph weights
  %   'TextureCoordinates' followed by
  %     UV  #V by 2 list of texture coordinates
  %
  % Examples:
  %   !python -m json.tool octopus.gltf > octopus-pretty.gltf
  NEAREST = 9728;
  LINEAR = 9729;
  NEAREST_MIPMAP_NEAREST = 9984;
  LINEAR_MIPMAP_NEAREST = 9985;
  NEAREST_MIPMAP_LINEAR = 9986;
  LINEAR_MIPMAP_LINEAR = 9987;
  CLAMP_TO_EDGE = 33071;
  MIRRORED_REPEAT = 33648;
  REPEAT = 10497;

  T = [];
  W = [];
  MV = [];
  MN = [];
  MW = [];
  N = [];
  fps = 30;
  nf = 0;
  UV = [];
  tex_im = [];
  mag_filter = LINEAR;
  min_filter = LINEAR_MIPMAP_LINEAR;
  wrap_s = MIRRORED_REPEAT;
  wrap_t = MIRRORED_REPEAT;

  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'FPS','MagFilter','MinFilter','MorphTargets','MorphNormals','MorphWeights','Normals','SkinningTransforms','SkinningWeights','TextureCoordinates','TextureImage'}, ...
    {'fps','mag_filter','min_filter','MV','MN','MW','N','T','W','UV','tex_im'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  mag_filter = string_to_magic_value(mag_filter);
  min_filter = string_to_magic_value(min_filter);
  function value = string_to_magic_value(string)
    if ischar(string)
      switch(string)
      case 'LINEAR'
        value = LINEAR;
      case 'NEAREST'
        value = NEAREST;
      case 'LINEAR_MIPMAP_NEAREST'
        value = LINEAR_MIPMAP_NEAREST;
      case 'LINEAR_MIPMAP_LINEAR'
        value = LINEAR_MIPMAP_LINEAR;
      otherwise
        error([string ' not found/implemented.']);
      end
    else
      value = string;
    end
  end

  % set flags
  has_normals = ~isempty(N);
  has_skinning_weights = ~isempty(W);
  has_skinning_transforms = ~isempty(T);
  has_morph_target_positions = ~isempty(MV);
  has_morph_target_weights = ~isempty(MW);
  has_morph_target_normals = ~isempty(MN);
  has_texture_coordinates = ~isempty(UV);
  has_texture_image = ~isempty(tex_im);

  % Determine what type to use for faces
  BYTE = 5120;
  UNSIGNED_BYTE = 5121;
  SHORT = 5122;
  UNSIGNED_SHORT = 5123;
  UNSIGNED_INT = 5125;
  FLOAT = 5126;
  ARRAY_BUFFER = 34962;
  ELEMENT_ARRAY_BUFFER = 34963;
  if max(F(:)) < 2^16
    F_size = 2;
    F_type = UNSIGNED_SHORT;
    F_matlab_type = 'uint16';
  else
    F_size = 4;
    F_type = UNSIGNED_INT;
    F_matlab_type = 'uint32';
  end

  % prepare things to write
  buffered_data = {};
  buffered_data{end+1} = struct('Name','F','Data',F,'Type','SCALAR','ComponentType',F_type,'Size',F_size,'MatlabType',F_matlab_type,'Target',ELEMENT_ARRAY_BUFFER);
  buffered_data{end+1} = struct('Name','V','Data',V,'Target',ARRAY_BUFFER);
  if has_normals
    buffered_data{end+1} = struct('Name','N','Data',normalizerow(N),'Target',ARRAY_BUFFER);
  end
  if has_morph_target_weights
    assert(size(MW,2) == size(MV,3));
    buffered_data{end+1} = struct('Name','MW','Data',MW,'Type','SCALAR');
    MWt = (0:size(MW,1)-1)*1/fps;
    buffered_data{end+1} = struct('Name','MWt','Data',MWt,'Type','SCALAR');
  end
  if has_morph_target_positions
    buffered_data{end+1} = struct('Name','MV','Data',MV,'Target',ARRAY_BUFFER);
  end
  if has_morph_target_normals
    buffered_data{end+1} = struct('Name','MN','Data',MN,'Target',ARRAY_BUFFER);
  end
  if has_skinning_weights
    [WIk,Wk] = prepare_skinning_weights(W);
    if max(WIk(:)) < 2^8
      WIk_size = 1;
      WIk_type = UNSIGNED_BYTE;
      WIk_matlab_type = 'uint8';
    else
      WIk_size = 2;
      WIk_type = UNSIGNED_SHORT;
      WIk_matlab_type = 'uint16';
    end
    buffered_data{end+1} = struct('Name','WIk','Data',WIk,'Size',WIk_size,'MatlabType',WIk_matlab_type,'ComponentType',WIk_type,'Target',ARRAY_BUFFER);
    buffered_data{end+1} = struct('Name','Wk','Data',Wk,'Target',ARRAY_BUFFER);
  end
  if has_skinning_transforms
    assert(size(W,2) == size(T,3));
    [tT,qU,sS,qVT] = prepare_skinning_transforms(T);
    buffered_data{end+1} = struct('Name','tT','Data',tT);
    buffered_data{end+1} = struct('Name','qU','Data',qU);
    buffered_data{end+1} = struct('Name','sS','Data',sS);
    buffered_data{end+1} = struct('Name','qVT','Data',qVT);
    Tt = (0:size(T,4)-1)*1/fps;
    buffered_data{end+1} = struct('Name','Tt','Data',Tt,'Type','SCALAR');
  end
  if has_texture_coordinates
    buffered_data{end+1} = struct('Name','UV','Data',UV,'Target',ARRAY_BUFFER);
  end

  [buffer,bufferViews,accessors,accessor_hash] = prepare_buffer(buffered_data);

  % First node is the mesh
  attributes = struct("POSITION",accessor_hash.V);

  if has_normals
    attributes.NORMAL = accessor_hash.N;
  end
  if has_texture_coordinates
    attributes.TEXCOORD_0 = accessor_hash.UV;
  end
  if has_skinning_weights
    attributes.JOINTS_0 = accessor_hash.WIk;
    attributes.WEIGHTS_0 = accessor_hash.Wk;
  end
  mesh = struct( ...
    "primitives",{{struct("attributes",attributes,"indices",accessor_hash.F)}});
  if has_morph_target_positions
    if ~isfield(mesh.primitives{1},'targets')
      mesh.primitives{1}.targets = num2cell(repmat(struct(),size(MV,3),1));
    end
    for i = 1:size(MV,3)
      mesh.primitives{1}.targets{i} = setfield(mesh.primitives{1}.targets{i},'POSITION',accessor_hash.MV(i));
    end
    mesh.weights = zeros(1,size(MV,3));
  end
  if has_morph_target_normals
    if ~isfield(mesh.primitives{1},'targets')
      mesh.primitives{1}.targets = num2cell(repmat(struct(),size(MN,3),1));
    end
    for i = 1:size(MN,3)
      mesh.primitives{1}.targets{i} = setfield(mesh.primitives{1}.targets{i},'NORMAL',accessor_hash.MN(i));
    end
  end

  nodes_list = {0};
  nodes = {struct("mesh",0)};
  if has_skinning_weights
    m = size(W,2);
    nodes{end}.skin = 0;
    % First put leaf nodes that will animate V-rotation
    for i = 1:m
      nodes{end+1} = struct( ...
        "name",sprintf('V-%d',i), ...
        "rotation",[0 0 0 1]);
    end
    % Then put parent auxiliary nodes that will animate translation, U-rotation, scale
    for i = 1:m
      nodes{end+1} = struct( ...
        "name",sprintf('TUS-%d',i), ...
        "children",{{i}}, ...
        "translation",[0 0 0], ...
        "rotation",[0 0 0 1], ...
        "scale",[1 1 1]);
    end
    %Put fake root at scene level and attach parents
    nodes_list{end+1} = 2*m+1;
    common_root = struct("name","root","children",{num2cell(m+(1:m))});
    nodes = {nodes{:} common_root};
  end

  % These {{double-brackets}} are needed when using the `struct` function.
  gltf = struct( ...
    "scene",0, ...
    "scenes",{{struct("nodes",{nodes_list})}}, ...
    "nodes",{nodes}, ...
    "meshes",{{mesh}}, ...
    "buffers",{{buffer}}, ...
    "bufferViews",{bufferViews}, ...
    "accessors",{accessors}, ...
    "asset",struct("version","2.0") ...
    );

  % Are these just for animation or do they overlap with texture samplers?
  samplers = {};
  channels = {};
  if has_skinning_weights
    gltf.skins = {struct('joints',{num2cell(0+1:m)})};
    if has_skinning_transforms
      for i = 1:m
        channels{end+1} = struct('sampler',numel(samplers),'target',struct('node',m+i,'path','translation'));
        samplers{end+1} = struct('input',accessor_hash.Tt,'output',accessor_hash.tT(i),'interpolation', 'STEP');
        channels{end+1} = struct('sampler',numel(samplers),'target',struct('node',m+i,'path','rotation'));
        samplers{end+1} = struct('input',accessor_hash.Tt,'output',accessor_hash.qU(i),'interpolation', 'STEP');
        channels{end+1} = struct('sampler',numel(samplers),'target',struct('node',m+i,'path','scale'));
        samplers{end+1} = struct('input',accessor_hash.Tt,'output',accessor_hash.sS(i),'interpolation', 'STEP');
        channels{end+1} = struct('sampler',numel(samplers),'target',struct('node',i,'path','rotation'));
        samplers{end+1} = struct('input',accessor_hash.Tt,'output',accessor_hash.qVT(i),'interpolation','STEP');
      end
    end
  end
  if has_morph_target_weights
    channels{end+1} = struct('sampler',numel(samplers),'target',struct('node',0,'path','weights'));
    samplers{end+1} = struct('input',accessor_hash.MWt,'output',accessor_hash.MW,'interpolation','LINEAR');
  end
  if ~isempty(samplers)
    gltf.animations = {struct('samplers',{samplers},'channels',{channels})};
  end
  if has_texture_image

    if mag_filter == LINEAR
      tmpfile = 'writeGLTF.jpg';
      imwrite(flipud(tex_im),tmpfile,'Quality',100);
    else
      tmpfile = 'writeGLTF.png';
      imwrite(flipud(tex_im),tmpfile);
    end
    tex_image_uri = imdata(tmpfile);
    delete(tmpfile);

    gltf.images = {struct('uri',tex_image_uri)};
    tex_image_index = numel(gltf.images)-1;
    if ~isfield(gltf,'samplers')
      gltf.samplers = {};
    end
    gltf.samplers{end+1} = struct('magFilter',mag_filter,'minFilter',min_filter,'wrapS',wrap_s,'wrapT',wrap_t);
    gltf.textures = {struct('sampler',numel(gltf.samplers)-1,'source',tex_image_index)};

    if has_texture_coordinates
      gltf.materials = {struct( ...
        'name','default', ...
        'pbrMetallicRoughness', ...
          struct( ...
            'baseColorTexture',struct('index',0), ...
            'metallicFactor',0, ...
            'roughnessFactor',1))};
      gltf.meshes{1}.primitives{1}.material = numel(gltf.materials)-1;
    end
  end

  fid = fopen(filename,'w');
  fprintf(fid,'%s',jsonencode(gltf));
  fclose(fid);

  function [buffer,bufferViews,accessors,accessor_hash] = prepare_buffer(buffered_data)
    buffer_index = 0;
    bufferViews = {};
    accessors = {};
    accessor_hash = struct();
    byte_count = 0;
    % write all data, tracking offsets
    tmpfile = 'writeGLTF.tmp';
    fp = fopen(tmpfile,'wb');
    for bi = 1:numel(buffered_data)
      item = buffered_data{bi};
      if isfield(item,'Size')
        adjusted_data = item.Data - 1;
        pad_amount = 4;
      else
        % default to single
        item.ComponentType = FLOAT;
        item.Size = 4;
        item.MatlabType = 'single';
        pad_amount = 8;
        adjusted_data = item.Data;
      end
      if ~isfield(item,'Type')
        switch(size(adjusted_data,2))
        case 2
          item.Type = 'VEC2';
        case 3
          item.Type = 'VEC3';
        case 4
          item.Type = 'VEC4';
        otherwise
          item.Type = 'SCALAR';
        end
      end

      accessor_ids = numel(accessors)-1 + (1:size(item.Data,3));
      for i = 1:size(item.Data,3)
        adjusted_data_i = adjusted_data(:,:,i);
        switch(item.Type)
        case 'SCALAR'
          accessor_count = numel(adjusted_data_i);
          accessor_min = {{min(  adjusted_data_i,[],[1 2])}};
          accessor_max = {{max(  adjusted_data_i,[],[1 2])}};
        case {'VEC2','VEC3','VEC4'}
          accessor_count = size( adjusted_data_i,1);
          accessor_min = min(    adjusted_data_i,[],1);
          accessor_max = max(    adjusted_data_i,[],1);
        end

        fwrite(fp,adjusted_data_i',item.MatlabType);
        % pad so that next block starts aligned at 4-byte
        bi_pad = mod(pad_amount-mod(numel(adjusted_data_i),pad_amount),pad_amount);
        fwrite(fp,zeros(bi_pad,1),item.MatlabType);
        % prepare bufferView and accessor
        % Always use 'buffer',0  → caller MUST adjust all if needed
        bufferViews{end+1} = struct( ...
          'name',item.Name, ...
          'buffer',buffer_index, ...
          'byteOffset',byte_count, ...
          'byteLength',numel(adjusted_data_i)*item.Size);
        if isfield(item,'Target')
          bufferViews{end} = setfield(bufferViews{end},'target',item.Target);
        end
        accessors{end+1} = struct( ...
          'name',item.Name, ...
          'bufferView',numel(bufferViews)-1, ...
          'byteOffset',0, ...
          'componentType',item.ComponentType, ...
          'type',item.Type, ...
          'count',accessor_count, ...
          'min',accessor_min, ...
          'max',accessor_max);
        % Increment byte_count
        byte_count = byte_count + (numel(adjusted_data_i)+bi_pad)*item.Size;
      end
      accessor_hash = setfield(accessor_hash,item.Name,accessor_ids);
    end
    fclose(fp);
    fp = fopen(tmpfile,'r');
    B = fread(fp,byte_count);
    fclose(fp);
    delete(tmpfile);
    Bstr = base64encode(B);
    %assert(isequal(Bstr,matlab.net.base64encode(uint8(B))))
    buffer = struct( ...
      'uri',['data:application/octet-stream;base64,' Bstr], ...
      'byteLength', numel(B));
  end

  function [WIk,Wk] = prepare_skinning_weights(W)
    % Number of skinning weights/handles/bones
    m = size(W,2) * ~isempty(W);
    % Convert weight matrix into multiple of 4 index-value pairs
    assert(size(W,2) == m);
    k = max(sum(W~=0,2));
    k4 = ceil(k/4)*4;
    % I'm not sure how >4 is handled.
    assert(k4 == 4);
    [Wk,WIk] = maxk(W,k,2);
    % Pad with null bones
    Wk(:,k+1:k4)  = 0;
    WIk(:,k+1:k4) = 0+1;%repmat(m+(1:k4-k),size(WIk,1),1);
    WIk(Wk==0) = 0+1;
  end

  function [tT,qU,sS,qVT] = prepare_skinning_transforms(T)
    m = size(T,3);
    nf = size(T,4) * ~isempty(T);
    T = permute(T,[1 2 4 3]);
    assert(size(T,2) == 4);
    tT = T(:,4,:,:);
    [sU,sS,sV] = pagesvd(T(1:3,1:3,:,:));
    % Push all of the reflection into sS
    % matlabFunction(det(sym('m',[3 3])))
    pagedet = @(m1_1,m1_2,m1_3,m2_1,m2_2,m2_3,m3_1,m3_2,m3_3) m1_1.*m2_2.*m3_3-m1_1.*m2_3.*m3_2-m1_2.*m2_1.*m3_3+m1_2.*m2_3.*m3_1+m1_3.*m2_1.*m3_2-m1_3.*m2_2.*m3_1;
    pagedet = @(m) pagedet(m(1,1,:,:),m(1,2,:,:),m(1,3,:,:),m(2,1,:,:),m(2,2,:,:),m(2,3,:,:),m(3,1,:,:),m(3,2,:,:),m(3,3,:,:));
    sUdet = sign(pagedet(sU));
    sVdet = sign(pagedet(sV));
    sU = sU.*sUdet;
    sV = sV.*sVdet;
    sS = sS.*sVdet.*sUdet;
    sS = cat(2,sS(1,1,:,:),sS(2,2,:,:),sS(3,3,:,:));

    qU = reshape(mat2quat(reshape(sU,[3 3 m*nf]))',[4 nf m]);
    qVT = reshape(mat2quat(reshape(permute(sV,[2 1 3 4]),[3 3 m*nf]))',[4 nf m]);
    for i = 1:m
      for f = 1:nf
        assert(norm(T(:,:,f,i) - [quat2mat(qU(:,f,i)')*diag(sS(:,:,f,i))*quat2mat(qVT(:,f,i)') tT(:,:,f,i)],inf)<1e-12);
      end
    end
    % gltf uses ijkw
    qU = qU([2 3 4 1],:,:);
    qVT = qVT([2 3 4 1],:,:);

    tT = permute(tT,[3 1 4 2]);
    sS = permute(sS,[3 2 4 1]);
    qU = permute(qU,[2 1 3]);
    qVT = permute(qVT,[2 1 3]);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function Conn = stat_threshConn(Conn,Stat,ThreshMethod,ThreshValue,remNonStatMethods)
% threshold connectivity values according to Stats structure
% Author: Tim Mullen, 2011

if nargin<3
    error('You must supply at least 3 arguments');
end

if nargin<4 && strcmpi(ThreshMethod,'pval')
    error('You must supply a p-value threshold');
end

if nargin<5
    remNonStatMethods = false;
end

statconnmethods = hlp_getConnMethodNames(Stat);
connmethods     = hlp_getConnMethodNames(Conn);

methodsToSkip   = ~ismember_bc(connmethods,statconnmethods);

if remNonStatMethods
    % remove all methods that we don't have stats for
    Conn = rmfield(Conn,connmethods(methodsToSkip));
end

for m=1:length(connmethods)
    if methodsToSkip(m)
        % skip thresholding for methods
        % that don't have associated Stats
        continue;
    end

    switch ThreshMethod
        case 'pval'
            % zero out any connections with p-value greater than
            % threshold
            Conn.(connmethods{m})(Stat.(connmethods{m}).(ThreshMethod) > ThreshValue) = 0;
        case 'thresh'
            % zero out any connections with amplitude less than
            % stat threshold

            Conn.(connmethods{m})( Conn.(connmethods{m}) < (Stat.(connmethods{m}).(ThreshMethod)) ) = 0;

        otherwise
            error('ThreshMethod must be ''pval'' or ''thresh''');
    end

end

repo: REST


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

[m, p] = max(sigmoid(X * all_theta'), [], 2);

% =========================================================================

end

author: zsiciarz

repo: ml-coursera


text:
function plotTree(tree, varNames)
%PLOTTREE   Plot a syntax tree (as stored in a TREEVAR object).
%   Calling sequence:
%      TREEVAR.PLOTTREE(TREE, VARNAMES)
%   where the input is:
%      TREE:     A MATLAB struct, describing the syntax tree of a mathematical
%                expression.
%      VARNAMES: An optional cellstring, whose elements are the name of the
%                variables that appear in a problem.
%
%   Usually, this method is called from within the TREEVAR plot() method.
%
%   Example:
%      % First, define a TREEVAR and carry out some operations:
%      u = treeVar();
%      v = cos(u);
%      w = sin(diff(u));
%      t = v + w;
%      % The following are equivalent:
%      plot(t)
%      treeVar.plotTree(t.tree)
%
% See also TREEVAR.PLOT.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Number of variables involved in the problem
numVars = length(tree.ID);

% Default variable names if none are passed:
if ( nargin < 2 )
    varNames = cell(1, numVars + 1);
    varNames{1} = 't';
    if ( numVars == 1 )
        % Simpler output for scalar case
        varNames{2} = 'u';
    else
        for varCounter = 1:numVars
            varNames{varCounter + 1} = sprintf('u%i', varCounter);
        end
    end
end

% Starting values for plotting
startx = .5;
deltax = .25;

% Layout the tree. The returned tree will store information about where
% it's supposed to be plotted.
tree = layoutNodes(tree, startx, deltax, tree.height + 1, tree.height + 1);

% Create a figure, and plot the first node
plot(tree.x, tree.y, '.','markersize', 25);
hold on

% Maximum differential orders that appear in the tree
maxDiffOrder = tree.diffOrder;

% Start the recursive calls to the actual plotting method:
plotTreePlot(tree, maxDiffOrder, varNames);

% Some final massaging of the figure:
xlim([0, 1])
ylim([0, 1])
hold off
set(gca, 'xtick', [])
set(gca, 'ytick', [])
title('Function evaluation tree')
end

function tree = layoutNodes(tree, treex, deltax, currHeight, maxHeight)
%LAYOUTNODES   Return coordinates for where to plot nodes of syntax trees.

if ( ~isstruct(tree) )
    % In case of scalars or CHEBFUNS appearing in the expression, the
    % corresponding nodes will not be syntax trees. Need to treat this case
    % separately.

    % Store the value of the current "tree":
    treeVal = tree;

    % Create a dummy tree, needed for storing the information required when the
    % tree is plotted:
    tree = [];
    tree.numArgs = 0;
    if ( isnumeric(treeVal) )
        tree.method = num2str(treeVal, '%2.2f');
    else
        tree.method = 'chebfun';
    end

    % Store coordinates where the node should be plotted:
    tree.y = currHeight/(maxHeight+1);
    tree.x = treex;

    % Scalars and CHEBFUNS have zero diffOrders.
    tree.diffOrder = 0;
    return
end

% Store coordinates where the current node should be plotted:
tree.y = currHeight/(maxHeight+1);
tree.x = treex;

% Lay out nodes recursively
switch tree.numArgs
    case 0
        % Do nothing
    case 1
        % Current tree has one child:
        tree.center = layoutNodes(tree.center, treex, deltax, ...
            currHeight - 1, maxHeight);
    case 2
        % Current tree has two childs:
        tree.left  = layoutNodes(tree.left, treex - deltax, ...
            deltax/2, currHeight - 1, maxHeight);
        tree.right  = layoutNodes(tree.right, treex + deltax, ...
            deltax/2, currHeight - 1, maxHeight);
end

end

function plotTreePlot(tree, maxDiffOrder, varNames)
%PLOTREEPLOT   The actual plotting of the syntax tree.

% Nice modern Matlab colours for lines and nodes:
lineColor = [0 0.447 0.741];
depVarColor = [0.85 0.325 0.098];
indVarColor = [0.929 0.694 0.125];

% Plot, depending on the number of arguments that the method has.
switch tree.numArgs
    case 0
        % Plot the constructor leaves in different colours.
        if ( strcmp(tree.method, 'constr') )
            if ( ~any(tree.ID) )
                % Independent variable
                col = indVarColor;
            else
                col = depVarColor;
            end
            plot(tree.x, tree.y, '.','markersize', 25, 'color', col);
        end

    case 1
        % Plot syntax tree for univariate methods:
        plot(tree.center.x, tree.center.y, '.','markersize', 25, ...
            'color', lineColor);
        plot([tree.x tree.center.x], [tree.y tree.center.y], 'color',lineColor)
        plotTreePlot(tree.center, maxDiffOrder, varNames)

    case 2
        % Plot syntax tree for bivariate methods. Start with left sub-tree.
        plot(tree.left.x, tree.left.y, '.','markersize', 25, ...
            'color', lineColor);
        plot([tree.x tree.left.x], [tree.y tree.left.y], 'color',lineColor)
        plotTreePlot(tree.left, maxDiffOrder, varNames)

        % Plot right sub-tree.
        plot(tree.right.x, tree.right.y, '.','markersize', 25, ...
            'color', lineColor);
        plot([tree.x tree.right.x], [tree.y tree.right.y], 'color',lineColor)
        plotTreePlot(tree.right, maxDiffOrder, varNames)
end

% Add the method name to the plot:
if ( strcmp(tree.method, 'constr') )
    % If we're at a constructor leaf, change the text slightly:
    varID = find(tree.ID == 1);
    if ( isempty(varID) )
        % Independent variable
        varString = varNames{1};
    else
        varString = varNames{varID + 1};
    end

    text(tree.x + 0.02, tree.y - 0.01, varString, 'Interpreter', 'none', ...
        'fontsize',14)
else
    text(tree.x + 0.02, tree.y - 0.01, tree.method, 'Interpreter', 'none', ...
        'fontsize',14)
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [all_surf_handles, pcl, ncl] = surface(obj, varargin)
% Render image data on brain surfaces; options for cutaways and canonical surfaces
%
% *Usage:*
%   [all_surf_handles, pcl, ncl] = surface(obj)
%   [all_surf_handles, pcl, ncl] = surface(r, ['cutaways', any optional inputs to surface_cutaway])
%
% This function uses region.surface to create surface figures.
% See help region.surface for options.
%
% :Examples:
% ::
%
%    % Create an initial surface plot from an fmri_data object:
%    han = surface(regionmasks{2});
%
%    % Now add a second region in green:
%    cluster_surf(region(regionmasks{2}), {[0 1 0]}, han, 5);
%
%    % Use optional arguments taken by surface_cutaway:
%    poscm = colormap_tor([1 .3 0], [1 1 0]); % orange to yellow
%    [all_surf_handles, pcl, ncl] = surface(t, 'cutaway', 'ycut_mm', -30, 'pos_colormap', poscm, 'existingfig');
%    [all_surf_handles2, pcl, ncl] = surface(t, 'foursurfaces', 'pos_colormap', poscm, 'neg_colormap', negcm);
%    [all_surf_handles2, pcl, ncl] = surface(t, 'foursurfaces', 'existingfig', 'color_upperboundpercentile', 95, 'color_lowerboundpercentile', 5, 'neg_colormap', colormap_tor([0 0 1], [.3 0 .5]));
%
%    % Use mediation_brain_surface_figs and re-make colors
%    all_surf_handles = mediation_brain_surface_figs([]);
%    surface(t2, 'cutaway', 'surface_handles', all_surf_handles, 'color_upperboundpercentile', 95, 'color_lowerboundpercentile', 5, 'neg_colormap', colormap_tor([0 0 1], [.2 0 .5]));
%

% :Programmers' notes:
% 9-8-2020 fixed weird bug that broke this somehow. Maybe someone deleted a line of code and committed?

if size(obj.dat, 2) > 1
    obj = mean(obj);
end

if isa(obj, 'atlas')

    r = atlas2region(obj, 'noverbose');

else

    r = region(obj, 'noverbose');

end

[all_surf_handles, pcl, ncl] = surface(r, varargin{:});

end

% Unique color option
% % Get unique color indices and expand to 256 colors for render_on_surface
% % expand
%
% obj.dat = double(obj.dat);
%
% len = max(obj.dat);
% rgbcell = scn_standard_colors(len);
%
% newvec = round(linspace(1, len, 256))';
%
% newcolors = rgbcell(newvec);
%
% newcm = cat(1, newcolors{:});
%
% render_on_surface(obj, hh, 'colormap', newcm);
%
% % works, but need to turn off edge effects!

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [C,PI,PC] = dual(V,F)
  % DUAL Construct the dual polygonal mesh of a given triangle mesh
  %
  % [C,PI,PC] = dual(V,F)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of indices into rows of V
  % Outputs:
  %   C  #C by dim list dual vertex positions
  %   PI  #PI stream of polygon indices into rows of C
  %   PC  #V+1 list of cumulative sum of dual face valences
  %
  % See also: polygons_to_triangles
  %
  assert(~any(on_boundary(F),'all'));
  [~,C] = circumradius(V,F);

  % vertex indices
  I = F(:);
  % only keep unique
  [I,J] = unique(I);
  % index in faces
  IF = mod(J-1,size(F,1))+1;
  % order in face
  IC = floor((J-1)/size(F,1))+1;

  [Fp,Fi] = triangle_triangle_adjacency(F);

  p1 = @(I) mod(I,3)+1;

  %NF = Fp(sub2ind(size(Fp),IF,p1(IC)));
  %NC = p1(Fi(sub2ind(size(Fi),IF,p1(IC))));
  %clf;
  %hold on;
  %tsurf(F,V,'FaceColor','w',falpha(0.8,1));
  %qvr(V(I,:),C(IF,:)-V(I,:),0,'LineWidth',2);
  %qvr(C(IF,:),C(NF,:)-C(IF,:),0,'LineWidth',2);
  %qvr(C(NF,:),V(F(sub2ind(size(F),NF,NC)),:)-C(NF,:),0,'LineWidth',2);
  %hold off;
  %axis equal;
  %view(3);

  % starting face
  IF0 = IF;
  % ledger of vertex-face pairs in order of observation
  L = [];
  while true
    L = [L;I IF];
    N = sub2ind(size(Fp),IF,p1(IC));
    NF = Fp(N);
    NC = p1(Fi(N));
    %clf;
    %hold on;
    %tsurf(F,V,'FaceColor','w',falpha(0.8,1));
    %qvr(V(I,:),C(IF,:)-V(I,:),0,'LineWidth',2);
    %qvr(C(IF,:),C(NF,:)-C(IF,:),0,'LineWidth',2);
    %qvr(C(NF,:),V(F(sub2ind(size(F),NF,NC)),:)-C(NF,:),0,'LineWidth',2);
    %hold off;
    %axis equal;
    %view(3);
    %pause
    IF = NF;
    IC = NC;

    keep = find(IF ~= IF0);
    IF = IF(keep);
    IC = IC(keep);
    I = I(keep);
    IF0 = IF0(keep);
    if isempty(keep)
      break;
    end
  end

  % stable sort by vertex
  [~,S] = sort(L(:,1));
  L = L(S,:);

  PC = cumsum([0;accumarray(L(:,1),1)]);
  PI = L(:,2);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function rinds = APPspInds2RegionInds(map, sinds)
% rinds = APPspInds2RegionInds(map, sinds)
% Gets the sp in each region
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Permission granted to non-commercial enterprises for
% modification/redistribution under GNU GPL.
% Current Version: 1.0  09/30/2005

nr = max(map);
rinds = cell(nr, 1);
for r = 1:nr
    count = 0;
    rs = find(map==r);
    for k = 1:length(rs)
        count = count + length(sinds{rs(k)});
    end
    rinds{r} = zeros(count, 1);

    count = 0;
    for k = 1:length(rs)
        rinds{r}(count+1:count+length(sinds{rs(k)})) = sinds{rs(k)};
        count = count + length(sinds{rs(k)});
    end
end

author: zouchuhang

repo: LayoutNet


text:
function jdate = julian (month, day, year)

% Julian date

% Input

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

% Output

%  jdate = Julian date

% special notes

%  (1) calendar year must include all digits

%  (2) will report October 5, 1582 to October 14, 1582
%      as invalid calendar dates and stop

% Orbital Mechanics with Matlab

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

y = year;
m = month;
b = 0;
c = 0;

if (m <= 2)
   y = y - 1;
   m = m + 12;
end

if (y < 0)
   c = -.75;
end

% check for valid calendar date

if (year < 1582)
   % null
elseif (year > 1582)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
elseif (month < 10)
   % null
elseif (month > 10)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
elseif (day <= 4)
   % null
elseif (day > 14)
   a = fix(y / 100);
   b = 2 - a + floor(a / 4);
else
   clc; home;

   fprintf('\n\n  this is an invalid calendar date!!\n');

   keycheck;

   return;
end

jd = fix(365.25 * y + c) + fix(30.6001 * (m + 1));

jdate = jd + day + b + 1720994.5;

repo: mcbench-benchmarks


text:
%% Transient diffusion equation
%% PDE and boundary conditions
% The transient diffusion equation reads
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc)
% , $D$ is the diffusion coefficient, and $\alpha$ is a constant.
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc

%% Define the domain and create a mesh structure
L = 50;  % domain length
Nx = 20; % number of cells
m = createMesh2D(Nx,Nx, L,L);
%% Create the boundary condition structure
BC = createBC(m); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=0; % left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
BC.top.a(:) = 0; BC.top.b(:)=1; BC.top.c(:)=0; % top boundary
BC.bottom.a(:) = 0; BC.bottom.b(:)=1; BC.bottom.c(:)=0; % bottom boundary
%% define the transfer coeffs
D_val = 1;
D = createCellVariable(m, D_val);
alfa = createCellVariable(m, 1);
%% define initial values
c_init = 1;
c_old = createCellVariable(m, c_init,BC); % initial values
c = c_old; % assign the old value of the cells to the current values
%% loop
dt = 1; % time step
final_t = 100;
for t=dt:dt:final_t
    [M_trans, RHS_trans] = transientTerm(c_old, dt, alfa);
    Dave = harmonicMean(D);
    Mdiff = diffusionTerm(Dave);
    [Mbc, RHSbc] = boundaryCondition(BC);
    M = M_trans-Mdiff+Mbc;
    RHS = RHS_trans+RHSbc;
    c = solvePDE(m,M, RHS);
    c_old = c;
    figure(1);visualizeCells(c);drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function uiwaitvec(h)
while ~isempty(h)
  uiwait(h(1));
  h = h(ishandle(h));
end

author: kristinbranson

repo: JAABA


%  Script file: doy.m
%
%  Purpose:
%    This program calculates the day of year corresponding
%    to a specified date.  It illustrates the use switch
%    and for constructs.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    01/27/07    S. J. Chapman        Original code
%
% Define variables:
%   day          -- Day (dd)
%   day_of_year  -- Day of year
%   ii           -- Loop index
%   leap_day     -- Extra day for leap year
%   month        -- Month (mm)
%   year         -- Year (yyyy)

% Get day, month, and year to convert
disp('This program calculates the day of year given the ');
disp('specified date.');
month = input('Enter specified month (1-12): ');
day   = input('Enter specified day(1-31):    ');
year  = input('Enter specified year(yyyy):   ');

% Check for leap year, and add extra day if necessary
if mod(year,400) == 0
   leap_day = 1;          % Years divisible by 400 are leap years
elseif mod(year,100) == 0
   leap_day = 0;          % Other centuries are not leap years
elseif mod(year,4) == 0
   leap_day = 1;          % Otherwise every 4th year is a leap year
else
   leap_day = 0;          % Other years are not leap years
end

% Calculate day of year by adding current day to the
% days in previous months.
day_of_year = day;
for ii = 1:month-1

   % Add days in months from January to last month
   switch (ii)
   case {1,3,5,7,8,10,12},
      day_of_year = day_of_year + 31;
   case {4,6,9,11},
      day_of_year = day_of_year + 30;
   case 2,
      day_of_year = day_of_year + 28 + leap_day;
   end

end

% Tell user
fprintf('The date %2d/%2d/%4d is day of year %d.\n', ...
         month, day, year, day_of_year);


text:
function [Y] = Rotation_ForceMoment_FromTireFixed(X,A_K_F,A_Rfl_F,A_Rfr_F,A_down_up)
%% Input parameters:
% A_K_F     [---]   Rotation Matrix from RearAxle Fixed Coordinate System to VehicleFixed Coordinate System
% A_Rfl_F   [---]   Rotation Matrix from RearAxle Fixed Coordinate System to WheelFixed Coordinate System FrontLeft
% A_Rfr_F   [---]   Rotation Matrix from RearAxle Fixed Coordinate System to WheelFixed Coordinate System FrontRight
% A_down_up [...]   Orientation Change Matrix between z-up Orientation and z-down Orientation
% X         [---]   Original Input Vectors [x1 x2 x3 x4] [3x4]
% Output parameters:
% Y         [---]   Transformed Output Vectors [y1 y2 y3 y4] [3x4]

%% Output of transformed vectors
% Quarter Vehicle Model along Vehicle Fixed z-Axis: Rotation from TireFixed Coordinate System into VehicleFixed Coordinate System
% Y=[(A_K_F*transpose(A_Rfl_F)*A_down_up*X(:,1)) (A_K_F*transpose(A_Rfr_F)*A_down_up*X(:,2)) (A_K_F*A_down_up*X(:,3)) (A_K_F*A_down_up*X(:,4))];

% Quarter Vehicle Model along Rear Axle Fixed z-Axis: Rotation from TireFixed Coordinate System into RearAxleFixed Coordinate System
Y=[(transpose(A_Rfl_F)*A_down_up*X(:,1)) (transpose(A_Rfr_F)*A_down_up*X(:,2)) (A_down_up*X(:,3)) (A_down_up*X(:,4))];

% Quarter Vehicle Model according MATLAB-Documentation
% Y=[(transpose(A_Rfl_F)*A_down_up*X(:,1)) (transpose(A_Rfr_F)*A_down_up*X(:,2)) (A_down_up*X(:,3)) (A_down_up*X(:,4))];

end

author: TUMFTM


text:
classdef DisplayFactory < handle

    methods (Access = public, Static)

        function display = create(chartType,title)
            switch chartType
                case 'plot'
                    display = Display_Plot(title);
                case 'log'
                    display = Display_SemilogY(title);
                case 'bar'
                    display = Display_Bar(title);
                otherwise
                    error('Invalid Chart Type.')
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function value = cn_geg_monomial_integral ( n, alpha, expon )

%*****************************************************************************80
%
%% CN_GEG_MONOMIAL_INTEGRAL: integral of monomial with Gegenbauer weight on CN.
%
%  Discussion:
%
%    CN_GEG is the cube [-1,+1]^N with the Gegenbauer weight function
%
%      w(alpha;x) = product ( 1 <= i <= n ) (1-x(i)^2)^alpha.
%
%    with -1.0 < alpha.
%
%    value = integral ( CN )
%      product ( 1 <= i <= n ) x(I)^expon(i) (1-x(i)^2)^alpha dx(i)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real ALPHA, the exponent of (1-X).
%    -1.0 < ALPHA.
%
%    Input, integer EXPON(N), the exponents.
%
%    Output, real VALUE, the value of the integral.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CN_GEG_MONOMIAL_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  ALPHA <= -1.0\n' );
    error ( 'CN_GEG_MONOMIAL_INTEGRAL - Fatal error!' );
  end

  value = 1.0;
  for i = 1 : n
    value2 = c1_geg_monomial_integral ( alpha, expon(i) );
    value = value * value2;
  end

  return

end

author: johannesgerer

repo: jburkardt-m


text:
close all; clear all
% pde = HodgeLaplacianEdata1;
pde = HodgeLaplacianFdata1;
[node,elem] = squaremesh([0,1,0,1],1/32);
% bdFlag = setboundary(node,elem,'Neumann');
% Pure Neumann boundary condition doesn't work.
% bdFlag = setboundary(node,elem,'Dirichlet');
bdFlag = setboundary(node,elem,'Dirichlet','x==0','Neumann','~(x==0)');

err = zeros(4,1); N = zeros(4,1);
option = [];
for i = 1:4
%     [sigma,u] = HodgeLaplacianE(node,elem,pde,bdFlag,option);
    [sigma,u,AD] = HodgeLaplacianF(node,elem,pde,bdFlag,option);
    err(i) = getL2error(node,elem,pde.exactsigma,sigma);
    N(i) = size(u,1);
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
end
showrate(N,err,2);

repo: Aero-matlab


text:
clear all; close all;

% The data
x=[1 2 3 4 5 6 7 8 9 10]
y=[0.2 0.5 0.3 3.5 1.0 1.5 1.8 2.0 2.3 2.2]

p1=fminsearch('fit1',[1 1],[],x,y);
p2=fminsearch('fit2',[1 1],[],x,y);
p3=fminsearch('fit3',[1 1],[],x,y);

xf=0:0.1:11
y1=polyval(p1,xf); y2=polyval(p2,xf); y3=polyval(p3,xf);

subplot(2,1,2)
plot(xf,y1,'k'), hold on
plot(xf,y2,'k--','Linewidth',[2])
plot(xf,y3,'k','Linewidth',[2])
plot(x,y,'ro','Linewidth',[2]), hold on

legend('E_\infty','E_1','E_2','location','NorthWest')
set(gca,'Fontsize',[15],'Ylim',[0 4],'Ytick',[0 1 2 3 4],'Xlim',[0 11],'Xtick',[0 2 4 6 8 10])

x=[1 2 3 4 5 6 7 8 9 10]
y=[0.2 0.5 0.3 0.7 1.0 1.5 1.8 2.0 2.3 2.2]

p1=fminsearch('fit1',[1 1],[],x,y);
p2=fminsearch('fit2',[1 1],[],x,y);
p3=fminsearch('fit3',[1 1],[],x,y);

xf=0:0.1:11
y1=polyval(p1,xf);
y2=polyval(p2,xf);
y3=polyval(p3,xf);

subplot(2,1,1)
plot(xf,y1,'k'), hold on
plot(xf,y2,'k--','Linewidth',[2])
plot(xf,y3,'k','Linewidth',[2])

plot(x,y,'ro','Linewidth',[2]), hold on

legend('E_\infty','E_1','E_2','location','NorthWest')
set(gca,'Fontsize',[15],'Ylim',[0 4],'Ytick',[0 1 2 3 4],'Xlim',[0 11],'Xtick',[0 2 4 6 8 10])

author: dynamicslab