text:
%%********************************************************************
%% infeaspt: generate an initial point for sdp.m
%%
%%  [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac);
%%
%%  options = 1  if want X0,Z0 to be scaled identity matrices
%%          = 2  if want X0,Z0 to be scalefac*(identity matrices).
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%********************************************************************

function [X0,y0,Z0] = infeaspt(blk,At,C,b,options,scalefac)
%%
if (nargin < 5); options = 1; end;
if (options == 1); scalefac = []; end;
if (options == 2) && (nargin < 6); scalefac = 1000; end;
if (scalefac <= 0); error('scalefac must a positive number'); end;
%%
if ~iscell(At); At = {At}; end;
if ~iscell(C);  C = {C}; end;
m = length(b);
if all(size(At) == [size(blk,1) m]);
    convertyes = zeros(size(blk,1),1);
    for p = 1:size(blk,1)
        if strcmp(blk{p,1},'s') && all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
        end
    end
    if any(convertyes)
        At = svec(blk,At,ones(size(blk,1),1));
    end
end;
%%
%%[blk,At,C,b] = validate(blk,At,C,b);
%%
X0 = cell(size(C)); Z0 = cell(size(C));
m = length(b);
for p = 1:size(blk,1);
    pblk = blk(p,:);
    blktmp = pblk{2};
    n = length(C{p});
    y0 = zeros(m,1);
    b2 = 1 + abs(b');
    if (options == 1);
        if strcmp(pblk{1},'s');
            normAni = [];
            X0{p} = sparse(n,n); Z0{p} = sparse(n,n);
            ss = [0, cumsum(blktmp)];
            tt = [0, cumsum(blktmp.*(blktmp+1)/2)];
            for i = 1:length(pblk{2})
                if ~isempty(At{p,1})
                    pos = tt(i)+1 : tt(i+1);
                    Ai = At{p,1}(pos,:);
                    normAni = 1+sqrt(sum(Ai.*Ai));
                end
                if (length(At(p,:)) >= 2) %% for low rank constraints
                    dd = At{p,3};
                    qq = [0, cumsum(pblk{3})]; normtmp = ones(1,length(pblk{3}));
                    idxD = [0; find(diff(dd(:,1))); size(dd,1)];
                    for k = 1:length(pblk{3})
                        idx = qq(k)+1 : qq(k+1);
                        idx2 = idxD(k)+1: idxD(k+1);
                        Ak = At{p,2}(:,idx);
                        ii = dd(idx2,2)-qq(k); %% undo cumulative indexing
                        jj = dd(idx2,3)-qq(k);
                        len = pblk{3}(k);
                        Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                        tmp = Ak'*Ak*Dk;
                        normtmp(1,k) = 1+sqrt(sum(sum(tmp.*tmp')));
                    end
                    normAni = [normAni, normtmp]; %#ok
                end
                pos = ss(i)+1 : ss(i+1);  ni = length(pos);
                tmp = C{p}(pos,pos);
                normCni = 1+sqrt(sum(sum(tmp.*tmp)));
                const  = 10; %%--- old: const = 1;
                constX = max([const,sqrt(ni),ni*(b2./normAni)]);
                constZ = max([const,sqrt(ni),normAni,normCni]);
                X0{p}(pos,pos) = constX*spdiags(1+1e-10*randmat(ni,1,0,'u'),0,ni,ni);
                Z0{p}(pos,pos) = constZ*spdiags(1+1e-10*randmat(ni,1,0,'u'),0,ni,ni);
            end
        elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            idenqX = zeros(sum(blktmp),1);
            idenqZ = zeros(sum(blktmp),1);
            idenqX(s(1:len)) = max([1,b2./normA])*sqrt(blktmp') ;
            idenqZ(s(1:len)) = max([sqrt(blktmp); max([normA,normC])*ones(1,len)])';
            idenqX(s(1:len)) = idenqX(s(1:len)).*(1+1e-10*randmat(len,1,0,'u'));
            idenqZ(s(1:len)) = idenqZ(s(1:len)).*(1+1e-10*randmat(len,1,0,'u'));
            X0{p} = idenqX;
            Z0{p} = idenqZ;
        elseif strcmp(pblk{1},'l');
            normC = 1+norm(C{p});
            normA = 1+sqrt(sum(At{p,1}.*At{p,1}));
            const = 10; %%--- old: const =1;
            constX = max([const,sqrt(n),sqrt(n)*b2./normA]);
            constZ = max([const,sqrt(n),normA,normC]);
            X0{p} = constX*(1+1e-10*randmat(n,1,0,'u'));
            Z0{p} = constZ*(1+1e-10*randmat(n,1,0,'u'));
        elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
        else
            error(' blk: some fields not specified correctly');
        end;
    elseif (options == 2);
        if strcmp(pblk{1},'s');
            n = sum(blktmp);
            X0{p} = scalefac*spdiags(1+1e-10*randmat(n,1,0,'u'),0,n,n);
            Z0{p} = scalefac*spdiags(1+1e-10*randmat(n,1,0,'u'),0,n,n);
        elseif strcmp(pblk{1},'q');
            s = 1+[0, cumsum(blktmp)];
            len = length(blktmp);
            idenq = zeros(sum(blktmp),1);
            idenq(s(1:len)) = 1+1e-10*randmat(len,1,0,'u');
            X0{p} = scalefac*idenq;
            Z0{p} = scalefac*idenq;
        elseif strcmp(pblk{1},'l');
            X0{p} = scalefac*(1+1e-10*randmat(n,1,0,'u'));
            Z0{p} = scalefac*(1+1e-10*randmat(n,1,0,'u'));
        elseif strcmp(pblk{1},'u');
            X0{p} = sparse(n,1);
            Z0{p} = sparse(n,1);
        else
            error(' blk: some fields not specified correctly');
        end
    end
end
%%********************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
  function  Quaternion = madgwickAHRS( Gyroscope, Accelerometer, Magnetometer, quaternion_l,Beta,SamplePeriod)
            q = quaternion_l; % short name local variable for readability
            % Normalise accelerometer measurement
            if(norm(Accelerometer) == 0), return; end	% handle NaN
            Accelerometer = Accelerometer / norm(Accelerometer);	% normalise magnitude

            % Normalise magnetometer measurement
            if(norm(Magnetometer) == 0), return; end	% handle NaN
            Magnetometer = Magnetometer / norm(Magnetometer);	% normalise magnitude

            % Reference direction of Earth's magnetic feild
            h = quaternProd(q, quaternProd([0 Magnetometer], quaternConj(q)));
            b = [0 norm([h(2) h(3)]) 0 h(4)];

            % Gradient decent algorithm corrective step
            F = [-2*(q(2)*q(4) - q(1)*q(3)) - Accelerometer(1)
                -2*(q(1)*q(2) + q(3)*q(4)) - Accelerometer(2)
                -2*(0.5 - q(2)^2 - q(3)^2) - Accelerometer(3)
                ((2*b(2)*(0.5 - q(3)^2 - q(4)^2) + 2*b(4)*(q(2)*q(4) - q(1)*q(3))) - Magnetometer(1))
                ((2*b(2)*(q(2)*q(3) - q(1)*q(4)) + 2*b(4)*(q(1)*q(2) + q(3)*q(4))) - Magnetometer(2))
                ((2*b(2)*(q(1)*q(3) + q(2)*q(4)) + 2*b(4)*(0.5 - q(2)^2 - q(3)^2)) - Magnetometer(3))];
            J = [2*q(3),                 	-2*q(4),                    2*q(1),                         -2*q(2)
                -2*q(2),                 	-2*q(1),                    	-2*q(4),                         -2*q(3)
                0,                         4*q(2),                    4*q(3),                         0
               -2*b(4)*q(3),               2*b(4)*q(4),               -4*b(2)*q(3)-2*b(4)*q(1),       -4*b(2)*q(4)+2*b(4)*q(2)
                -2*b(2)*q(4)+2*b(4)*q(2),	2*b(2)*q(3)+2*b(4)*q(1),	2*b(2)*q(2)+2*b(4)*q(4),       -2*b(2)*q(1)+2*b(4)*q(3)
                2*b(2)*q(3),                2*b(2)*q(4)-4*b(4)*q(2),	2*b(2)*q(1)-4*b(4)*q(3),        2*b(2)*q(2)];
            step = (J'*F);
            step = step / norm(step);	% normalise step magnitude

            % Compute rate of change of quaternion
            qDot = 0.5 * quaternProd(q, [0 Gyroscope]) - Beta * step';

            % Integrate to yield quaternion
            q = q + qDot * SamplePeriod;
            Quaternion = q / norm(q); % normalise quaternion
        end

repo: GPS-INS-Integrated-Navigation


text:
function model = createToyModelForGeoFBA()
%createToyModelForMDFBA creates a toy model for
%MDFBA
% The model created looks as follows:
%
%
%   <-> A -> B ---> C --> E <->
%        \                ^
%         \              /
%           --> F --> G
%
%
% A normal FBA maximizing the production of E should not yield any flux
% through A -> D
% MDFBA should show this flux.

model = createModel();
%Reactions in {Rxn Name, MetaboliteNames, Stoichiometric coefficient} format
Reactions = {'R1',{'A','B'},[-1 1];...
             'R2',{'B','C'},[-1  1];...
             'R3',{'C','E'},[-1  1];...
             'R4',{'A','F'},[-1 1];...
             'R6',{'F','G'},[-1 1];...
             'R7',{'G','E'},[-1 1]};
ExchangedMets = {'A','E'};
%Add Reactions
for i = 1:size(Reactions,1)
    %All reactions are irreversible
    model = addReaction(model,Reactions{i,1},Reactions{i,2},Reactions{i,3},0,0,1000);
end

%Add Exchangers
model = addExchangeRxn(model,ExchangedMets,-1000*ones(numel(ExchangedMets),1),1000*ones(numel(ExchangedMets),1));
model = changeObjective(model,'EX_E',1);

author: opencobra

repo: cobratoolbox


text:
function [ output_data ] = normalize_complex_mem( input_data, ...
    DeltaKCOAPoly, dim1_coords_m, dim2_coords_m, ...
    weight_fun, oversample_rate, ...
    fftsign, dim )
%NORMALIZE_COMPLEX_MEM Normalize complex SAR data to bring to a consistent
%point in the processing chain
%
%    output_data = normalize_complex_mem(input_data, ...
%       DeltaKCOAPoly, dim1_coords_m, dim2_coords_m, ... % Deskew parameters
%       weight_fun, oversample_rate, ... % Deweighting paramters
%       fftsign, dim)
%
%       Parameter name    Description
%
%       input_data        Array of complex values to deskew
%       DeltaKCOAPoly     Polynomial that describes center of frequency
%                            support of data in the processing dimension.
%                            Described in the SICD design and exploitation
%                            document.  DeltaKCOAPoly of 0 means no deskew.
%       dim1_coords_m     Coordinate of each "row" in dimension 1.  Used
%                            only for deskew.
%       dim2_coords_m     Coordinate of each "column" in dimension 2.  Used
%                            only for deskew.
%       weight_fun        Description of weighting applied.  Either a
%                            function handle to a function that takes a
%                            single argument (number of elements) and
%                            produces the weighting to apply, or a vector
%                            that is the weighting function sampled.  An
%                            empty weight_fun means to apply no
%                            deweighting.
%       oversample_rate   Amount of sampling beyond the ImpRespBW in the
%                            processing dimension.  (Nyquist sampling = 1).
%                            Used only for deweighting.
%       fftsign           FFT sign in the processing dimension (-1 or +1).
%       dim               Dimension over which to perform normalization
%       output_data       INPUT_DATA with normalization applied
%
% Author: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

output_data = input_data;
if any(DeltaKCOAPoly(:)~=0) % Do we need to center frequency support?
    output_data = deskewmem(input_data, DeltaKCOAPoly, dim1_coords_m, dim2_coords_m, dim, fftsign);
end
if fftsign==1 % -1 is the norm
    output_data = conj(output_data);
    % This transform will maintain the magnitude of the data and make the
    % IFFT behave like an FFT (and vice versa) in magnitude only:
    % abs(x)==abs(conj(x)) and
    % abs(fft(x))==abs(ifft(conj(x))), since fft(x)==conj(ifft(conj(x)))
end
if ~isempty(weight_fun) % Empty array means uniform weighting
    output_data = deweightmem(output_data, weight_fun, oversample_rate, dim);
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function fbcEnabled = isFbcEnabled()
% Checks whether the version of libSBML has been built with
% the FBC package extension enabled

% Filename    : isFbcEnabled.m
% Description : check fbc status
% Author(s)   : SBML Team <sbml-team@googlegroups.com>
% Organization: EMBL-EBI, Caltech
% Created     : 2011-02-08
%
% This file is part of libSBML.  Please visit http://sbml.org for more
% information about SBML, and the latest version of libSBML.
%
% Copyright (C) 2013-2018 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%     3. University of Heidelberg, Heidelberg, Germany
%
% Copyright (C) 2009-2011 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%
% Copyright (C) 2006-2008 by the California Institute of Technology,
%     Pasadena, CA, USA
%
% Copyright (C) 2002-2005 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. Japan Science and Technology Agency, Japan
%
% This library is free software; you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation.  A copy of the license
% agreement is provided in the file named "LICENSE.txt" included with
% this software distribution and also available online as
% http://sbml.org/software/libsbml/license.html

% assume not enabled
fbcEnabled = 0;

if (isoctave() == '0')
  filename = fullfile(tempdir, 'fbc.xml');
else
  filename = fullfile(pwd, 'fbc.xml');
end;
writeTempFile(filename);

try
  [m, e] = TranslateSBML(filename, 1, 0);

  if (length(e) == 0 && isfield(m, 'fbc_version') == 1 )
    fbcEnabled = 1;
  end;

  delete(filename);

catch

  delete(filename);

  return
end;

function writeTempFile(filename)

fout = fopen(filename, 'w');

fprintf(fout, '<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n');
fprintf(fout, '<sbml xmlns=\"http://www.sbml.org/sbml/level3/version1/core\" ');
fprintf(fout, 'xmlns:fbc=\"http://www.sbml.org/sbml/level3/version1/fbc/version1\" ');
fprintf(fout, 'level=\"3\" version=\"1\" fbc:required=\"false\">\n');
fprintf(fout, '  <model/>\n</sbml>\n');

fclose(fout);


text:
%KCENTRES Finds K center objects from a distance matrix
%
%  [LAB,J,DM] = KCENTRES(D,K,N)
%
% INPUT
%   D    Distance matrix between, e.g. M objects (may be a dataset)
%   K    Number of center objects to be found (optional; default: 1)
%   N    Number of trials starting from a random initialization
%        (optional; default: 1)
%
% OUTPUT
%   LAB  Integer labels: each object is assigned to its nearest center
%   J    Indices of the center objects
%   DM   A list of distances corresponding to J: for each center in J
%        the maximum distance of the objects assigned to this center.
%
% DESCRIPTION
% Finds K center objects from a symmetric distance matrix D. The center
% objects are chosen from all M objects such that the maximum of the
% distances over all objects to the nearest center is minimized. For K > 1,
% the results depend on a random initialization. The procedure is repeated
% N times and the best result is returned.
%
% If N = 0, initialisation is not random, but done by a systematic
% selection based on a greedy approach.
%
% SEE ALSO (PRTools Guide)
% HCLUST, PRKMEANS, EMCLUST, MODESEEK

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: kcentres.m,v 1.5 2008/07/03 09:08:43 duin Exp $

function [labels,Jopt,dm] = kcentres(d,k,n)

	if (nargin < 3) | isempty(n)
      n = 1;
      prwarning(4,'Number of trials not supplied, assuming one.');
   end

	if (nargin < 2) | isempty(k)
      k = 1;
      prwarning(4,'Number of centers not supplied, assuming one.');
   end

	if(isdataset(d))
		d = +d;
		prwarning(4,'Distance matrix is convert to double.');
	end

	[m,m2] = size(d);
	if ( ~issym(d,1e-12) | any(diag(d) > 1e-14) )
		error('Distance matrix should be symmetric and have zero diagonal')
	end

% checking for a zero diagonal
	t = eye(m) == 1;
	if(~all(d(t)==0))
		error('Distance matrix should have a zero diagonal.')
	end

	if (k == 1)
		dmax = max(d);
		[dm,Jopt] = min(dmax);
		labels = repmat(1,m,1);
		return;
	end

	if k > m
		error('Number of centres should not exceed number of objects')
	end

	% We are here only if K (> 1) centers are to be found.
	% Loop over number of trials.
	dmax = max(max(d));
	dopt = inf;
	s = sprintf('k-centres, %i attempts: ',n);
	prwaitbar(n,s,n>1);
	if n == 0
		nrep = 1;
	else
		nrep = n;
	end
	for tri = 1:nrep
		prwaitbar(n,tri,[s int2str(tri)]);
		if n == 0
			M = kcentresort(d,k);           % systematic initialisation
		else
			M = randperm(m); M = M(1:k);		% Random initializations
		end
 		J = zeros(1,k); 					    	  % Center objects to be found.

		% Iterate until J == M. See below.
		while 1,
			[dm,I] = min(d(M,:));

			% Find K centers.
			for i = 1:k
				JJ = find(I==i);
				if (isempty(JJ))
					%JJ can be empty if two or more objects are in the same position of
					% feature space in dataset
					J(i) = 0;
				else
						% Find objects closest to the object M(i)
					[dummy,j,dm] = kcentres(d(JJ,JJ),1,1);
					J(i) = JJ(j);
				end
			end

			J(find(J==0)) = [];
			k = length(J);
			if k == 0
				error('kcentres fails as some objects are identical: add some noise')
			end
			if (length(M) == k) & (all(M == J))
				% K centers are found.
				[dmin,labs] = min(d(J,:));
				dmin = max(dmin);
				break;
			end
			M = J;
		end

		% Store the optimal centers in JOPT.

		if (dmin <= dopt)
			dopt = dmin;
			labels = labs';
			Jopt = J;
		end

	end
	prwaitbar(0)

	% Determine the best centers over the N trials.
	dm = zeros(1,k);
	for i=1:k
		L = find(labels==i);
		dm(i) = max(d(Jopt(i),L));
	end

return;

%KCENTRESORT Sort objects given by dissimilarity matrix
%
%               N = KCENTRESORT(D,P,CRIT)
%
% INPUT
%   D     Square dissimilarity matrix, zeros on diagonal
%   P     Number of prototypes to be selected
%   CRIT  'dist' or 'centre'
%
% OUTPUT
%   N     Indices of selected prototypes
%
% DESCRIPTION
% Sort objects given by square dissim matrix D using a greedy approach
% such that the maximum NN distance from all objects (prototypes)
% to the first K: max(min(D(:,N(1:K),[],2)) is minimized.
%
% This routines tries to sample the objects such that they are evenly
% spaced judged from their dissimilarities. This may be used as
% initialisation in KCENTRES. It works reasonably, but not very good.
%
% SEE ALSO (PRTools Guide)
% KCENTRES

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function N = kcentresort(d,p,crit);

d = +d;
m = size(d,1);
if nargin < 3, crit = 'dist'; end
if nargin < 2 | isempty(p), p = m; end
L = [1:m];
N = zeros(1,p);
[dd,n] = min(max(d,[],2));  % this is the first (central) prototype
e = d(:,n);                 % store here the distances to the nearest prototype (dNNP)
f = min(d,repmat(e,1,m));   % replace distances that are larger than dNNP by dNNP
N(1) = n;                   % ranking of selected prototypes
L(n) = [];                  % candidate prototypes (all not yet selected objects)

for k=2:p                   % extend prototype set
        if strcmp(crit,'centre')
        [dd,n] = min(max(f(L,L),[],1)); % This selects the next prototype out of candidates in L
        e = min([d(:,L(n)) e],[],2);    % update dNNP
        f = min(d,repmat(e,1,m));       % update replacement of distances that are larger
                                              % than dNNP by dNNP
        elseif strcmp(crit,'dist')
                [dd,n] = max(mean(d(N(find(N > 0)),L)));
        else
                error('Illegal crit')
        end
  N(k) = L(n);                    % update list of selected prototypes
  L(n) = [];                      % update list of candidate prototypes
end

author: marianux

repo: ecg-kit

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


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

%%
clear; close all; clc;

%% Syntax
% |[F,V,logicValid,indFix2]=patchEdgeCollapse(F,V,E,logicKeep,meanOption);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function obj = addParamVariable(obj, bounds)
    % Adds parameters as the NLP decision variables to the problem
    %
    % Parameters:
    %  bounds: a structed data stores the boundary information of the
    %  NLP variables @type struct

    % get basic information of the variables
    params = obj.Plant.Params;

    param_names = fieldnames(params);

    for j=1:length(param_names)

        p_name = param_names{j};

        var = struct();
        var.Name = p_name;
        siz = size(params.(p_name));
        var.Dimension = prod(siz); %#ok<PSIZE>
        if isfield(bounds,p_name)
            if isfield(bounds.(p_name),'lb')
                var.lb = bounds.(p_name).lb;
            end
            if isfield(bounds.(p_name),'ub')
                var.ub = bounds.(p_name).ub;
            end
            if isfield(bounds.(p_name),'x0')
                var.x0 = bounds.(p_name).x0;
            end
        end
        % determines the nodes at which the variables to be defined.
        if obj.Options.DistributeParameters
            % time variables are defined at all nodes if distributes the weightes
            obj = addVariable(obj, p_name, 'all', var);

            % add an equality constraint between the time variable at
            % neighboring nodes to make sure they are same
            p   = flatten(params.(p_name));
            siz = size(p);
            pn  = SymVariable([p_name 'n'],siz);
            p_cont = SymFunction([p_name 'Cont_' obj.Plant.Name],transpose(p-pn),{p,pn});
            % create an array of constraints structure
            p_cstr = struct();
            p_cstr(obj.NumNode-1) = struct();
            [p_cstr.Name] = deal(p_cont.Name);
            [p_cstr.Dimension] = deal(var.Dimension);
            [p_cstr.lb] = deal(0);
            [p_cstr.ub] = deal(0);
            [p_cstr.Type] = deal('Linear');
            [p_cstr.SymFun] = deal(p_cont);
            for i=1:obj.NumNode-1
                p_cstr(i).DepVariables = [obj.OptVarTable.(p_name)(i);obj.OptVarTable.(p_name)(i+1)];
            end

            % add to the NLP constraints table
            obj = addConstraint(obj,[p_name 'Cont'],'except-last',p_cstr);

        else
            % otherwise only define at the first node
            obj = addVariable(obj, p_name, 'first', var);
        end

    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [EV, EVal] = ncuts(A, n_ev)
% Computes the n_ev smallest (non-zero) eigenvectors and eigenvalues of the
% of the Laplacian of A

D = sparse(1:size(A,1), 1:size(A,1), full(sum(A, 1)), size(A,1), size(A,2));

opts.issym = 0;
opts.isreal = 1;
opts.disp = 0;
nvec = n_ev+1;

[EV, EVal] = eigs((D - A) + (10^-10) * speye(size(D)), D, nvec, 'sm',opts);

[junk, sortidx] = sort(diag(EVal), 'descend');
EV = EV(:,sortidx(end-1:-1:1));
v = diag(EVal);
EVal = v(sortidx(end-1:-1:1));

EV = bsxfun(@rdivide, EV, sqrt(sum(EV.^2,1))); % makes the eigenvectors unit norm

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
%% Copyright (C) 2014, 2016, 2018-2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun findsymbols (@var{x})
%% Return a list (cell array) of the symbols in an expression.
%%
%% The list is sorted alphabetically.  For details, @pxref{@@sym/symvar}.
%%
%% If two variables have the same symbol but different assumptions,
%% they will both appear in the output.  It is not well-defined
%% in what order they appear.
%%
%% @var{x} could be a sym, sym array, cell array, or struct.
%%
%% @example
%% @group
%% syms x y z
%% C = @{x, 2*x*y, [1 x; sin(z) pi]@};
%% S = findsymbols (C)
%%   @result{} S = @{ ... @}
%% S@{:@}
%%   @result{} ans = (sym) x
%%   @result{} ans = (sym) y
%%   @result{} ans = (sym) z
%% @end group
%% @end example
%%
%% Note ℯ, ⅈ, π, etc are not considered as symbols.
%%
%% Note only returns symbols actually appearing in the RHS of a
%% @code{symfun}.
%%
%% @seealso{symvar, @@sym/symvar, @@sym/findsym}
%% @end defun

function L = findsymbols(obj, dosort)

  if (nargin == 1)
    dosort = true;
  elseif (nargin ~= 2)
    print_usage ();
  end

  if isa(obj, 'sym')
    cmd = { 'if isinstance(_ins[0], MatrixBase):'
            '    # having problems with Eq in matrices in SymPy >= 1.7.1'
            '    s = set().union(*(i.free_symbols for i in _ins[0] if i is not None))'
            'else:'
            '    s = _ins[0].free_symbols'
            'l = list(s)'
            'l = sorted(l, key=str)'
            'return l,' };
    L = pycall_sympy__ (cmd, obj);

  elseif iscell(obj)
    %fprintf('Recursing into a cell array of numel=%d\n', numel(obj))
    L = {};
    for i=1:numel(obj)
      temp = findsymbols(obj{i}, false);
      if ~isempty(temp)
        L = {L{:} temp{:}};
      end
    end

  elseif isstruct(obj)
    %fprintf('Recursing into a struct array of numel=%d\n', numel(obj))
    L = {};
    fields = fieldnames(obj);
    for i=1:numel(obj)
      for j=1:length(fields)
        thisobj = getfield(obj, {i}, fields{j});
        temp = findsymbols(thisobj, false);
        if ~isempty(temp)
          L = {L{:} temp{:}};
        end
      end
    end

  else
    L = {};
  end

  % sort and make unique using internal representation
  if dosort
    Ls = {};
    for i=1:length(L)
      Ls{i} = sympy (L{i});
    end
    [tilde, I] = unique(Ls);
    L = L(I);
  end
end

%!test
%! syms x b y n a arlo
%! z = a*x + b*pi*sin (n) + exp (y) + exp (sym (1)) + arlo;
%! s = findsymbols (z);
%! assert (isequal ([s{:}], [a,arlo,b,n,x,y]))
%!test
%! syms x
%! s = findsymbols (x);
%! assert (isequal (s{1}, x))
%!test
%! syms z x y a
%! s = findsymbols ([x y; 1 a]);
%! assert (isequal ([s{:}], [a x y]))
%!assert (isempty (findsymbols (sym (1))))
%!assert (isempty (findsymbols (sym ([1 2]))))
%!assert (isempty (findsymbols (sym (nan))))
%!assert (isempty (findsymbols (sym (inf))))
%!assert (isempty (findsymbols (exp (sym (2)))))

%!test
%! % empty sym for findsymbols, findsym, and symvar
%! assert (isempty (findsymbols (sym([]))))
%! assert (isempty (findsym (sym([]))))
%! assert (isempty (symvar (sym([]))))

%!test
%! % diff. assumptions make diff. symbols
%! x1 = sym('x');
%! x2 = sym('x', 'positive');
%! f = x1*x2;
%! assert (length (findsymbols (f)) == 2)

%!test
%! % symfun or sym
%! syms x f(y)
%! a = f*x;
%! b = f(y)*x;
%! assert (isequal (findsymbols(a), {x y}))
%! assert (isequal (findsymbols(b), {x y}))

%!test
%! % findsymbols on symfun does not find the argnames (unless they
%! % are on the RHS of course, this matches SMT 2014a).
%! syms a x y
%! f(x, y) = a;  % const symfun
%! assert (isequal (findsymbols(f), {a}))
%! syms a x y
%! f(x, y) = a*y;
%! assert (isequal (findsymbols(f), {a y}))

%!test
%! % sorts lexigraphically, same as symvar *with single input*
%! % (note symvar does something different with 2 inputs).
%! syms A B a b x y X Y
%! f = A*a*B*b*y*X*Y*x;
%! assert (isequal (findsymbols(f), {A B X Y a b x y}))
%! assert (isequal (symvar(f), [A B X Y a b x y]))

%!test
%! % symbols in matpow
%! syms x y
%! syms n
%! A = [sin(x) 2; y 1];
%! B = A^n;
%! L = findsymbols(B);
%! assert (isequal (L, {n x y}))

%!test
%! % array of eq
%! syms x y
%! assert (isequal (findsymbols (2 == [2 x y]), {x y}))

%!test
%! % array of ineq
%! syms x y
%! A = [x < 1  2*x < y  x >= 2  3 <= x  x != y];
%! assert (isequal (findsymbols (A), {x y}))

repo: octsympy


text:
% Read data from files generated by HList
% Author: Xiao Xiong
% Date Created: 25 Jan 2005
% Date Lasted modefied: 25 Jan 2005

function x = readHList(file_name, vector_size);

    % read in the log Mel vectors first
    FILE = fopen( file_name );
    j=0;
    fsearch(':', FILE);     % skip the first ':'
    while fsearch(':', FILE) == 0   % read N_spectral number immediately after ':'
        j = j+1;
        temp = textscan(FILE, '%n', vector_size);
        x(:,j) = temp{1};        % store the vecors as column vector of logmel
    end
    num_of_vector_in_current_file = j;

author: singaxiong

repo: SignalGraph


text:
%Install_lightspeed
% Compiles mex files for the lightspeed library.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

% thanks to Kevin Murphy for suggesting this routine.
% thanks to Ruben Martinez-Cantin for UNDERSCORE_LAPACK_CALL

fprintf('Compiling lightspeed 2.6 mex files...\n');
fprintf('Change directory to lightspeed for this to work.\n');

% Matlab version
v = sscanf(version,'%d.%d.%*s (R%d) %*s');
% v(3) is the R number
% could also use v(3)>=13
atleast65 = (v(1)>6 || (v(1)==6 && v(2)>=5));
atleast73 = (v(1)>7 || (v(1)==7 && v(2)>=3));
atleast75 = (v(1)>7 || (v(1)==7 && v(2)>=5));
atleast76 = (v(1)>7 || (v(1)==7 && v(2)>=6));
atleast78 = (v(1)>7 || (v(1)==7 && v(2)>=8)); % R2009a
if atleast73
	% largeArrayDims and mwSize were added in version 7.3 (R2006b)
	% http://www.mathworks.com/help/techdoc/rn/bqt6wtq.html
	flags = ' -largeArrayDims ';
else
	flags = ' -DmwSize=int -DmwIndex=int ';
end

% copy matlab's original repmat.m as xrepmat.m
if exist('xrepmat.m') ~= 2
  w = fullfile(matlabroot,'toolbox','matlab','elmat','repmat.m');
  cmd = ['"' w '" xrepmat.m'];
  if ispc
    system(['copy ' cmd]);
  else
    system(['cp -rp ' cmd]);
  end
end

% these are done first to initialize mex
eval(['mex' flags '-c flops.c']);
eval(['mex' flags 'sameobject.c']);
eval(['mex' flags 'int_hist.c']);
eval(['mex' flags '-c mexutil.c']);
eval(['mex' flags '-c util.c']);

libdir = '';
if ispc
	[compiler,options] = mexcompiler;
	libdir = options.LIBLOC;
	engmatopts = [compiler 'engmatopts.bat'];
elseif ismac
	options = struct;
	% this installer is set up for 64-bit MacOSX 10.6 with gcc-4.0
	% if you are using something else, run 'mex -v -c flops.c'
	% and use the output to change these strings
	options.COMPILER = 'gcc-4.0';
	options.COMPFLAGS = '-fno-common -no-cpp-precomp -arch x86_64 -isysroot /Developer/SDKs/MacOSX10.6.sdk -mmacosx-version-min=10.6  -fexceptions';
	options.OPTIMFLAGS = '-O -DNDEBUG';
else
	options = struct;
	options.COMPILER = 'cc';
	options.COMPFLAGS = '-fPIC';
	options.OPTIMFLAGS = '-O';
end

% Routines that use LAPACK
lapacklib = '';
blaslib = '';
lapackflags = flags;
if atleast78
	lapackflags = [lapackflags ' -DBLAS64'];
end
if ispc
  if strncmp(compiler,'MSVC',4)
    if atleast65
      % version >= 6.5
      lapacklib = fullfile(libdir,'libmwlapack.lib');
    end
  else
    lapacklib = fullfile(libdir,'libmwlapack.lib');
  end
  if atleast75
    blaslib = fullfile(libdir,'libmwblas.lib');
  end
  %%% Paste the location of libmwlapack.lib %%%
  %lapacklib = '';
  if ~exist(lapacklib,'file')
    lapacklib = 'dtrsm.c';
    fprintf('libmwlapack.lib was not found.  To get additional optimizations, paste its location into install_lightspeed.m\n');
  else
    fprintf('Using the lapack library at %s\n',lapacklib);
  end
else
  % in version 7.5, non-PC systems do not need to specify lapacklib,
  % but they must use an underscore when calling lapack routines
	% http://www.mathworks.com/help/techdoc/matlab_external/br_2m24-1.html
	lapackflags = [lapackflags ' -DUNDERSCORE_LAPACK_CALL'];
  if atleast76
    lapacklib = '-lmwlapack';
    blaslib = '-lmwblas';
  end
end
disp(['mex' lapackflags ' solve_triu.c "' lapacklib '" "' blaslib '"'])
eval(['mex' lapackflags ' solve_triu.c "' lapacklib '" "' blaslib '"']);
eval(['mex' lapackflags ' solve_tril.c "' lapacklib '" "' blaslib '"']);

if ispc
  % Windows
  %if exist('util.obj','file')
  eval(['mex' flags 'addflops.c flops.obj'])
	if atleast78
		eval(['mex' flags 'gammaln.c util.obj -outdir @double'])
	else
		eval(['mex' flags 'gammaln.c util.obj'])
	end
  eval(['mex' flags 'digamma.c util.obj'])
  eval(['mex' flags 'trigamma.c util.obj'])
  eval(['mex' flags 'tetragamma.c util.obj'])
	eval(['mex' flags 'setnonzeros.c'])
  if strncmp(compiler,'MSVC',4)
		clear random.dll randomseed randbinom randgamma sample_hist
    disp(['install_random.bat "' options.VSINSTALLDIR '" ' options.vcvarsopts]);
    system(['install_random.bat "' options.VSINSTALLDIR '" ' options.vcvarsopts]);
    eval(['mex' flags 'randomseed.c util.obj random.lib'])
    eval(['mex' flags 'randbinom.c mexutil.obj util.obj random.lib'])
    eval(['mex' flags 'randgamma.c mexutil.obj util.obj random.lib'])
    eval(['mex' flags 'sample_hist.c util.obj random.lib'])
  else
    fprintf('mexcompiler is not MSVC. The randomseed() function will have no effect.');
    eval(['mex' flags 'randomseed.c util.obj random.c'])
    eval(['mex' flags 'randbinom.c mexutil.obj util.obj random.c'])
    eval(['mex' flags 'randgamma.c mexutil.obj util.obj random.c'])
    eval(['mex' flags 'sample_hist.c util.obj random.c'])
  end
  eval(['mex' flags 'repmat.c mexutil.obj'])
  try
    % standalone programs
    % compilation instructions are described at:
    % http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_external/ch1_im15.html#27765
    if atleast78
			disp('lightspeed''s matfile utility is not supported for this version of Matlab');
		elseif atleast65
      % -V5 is required for Matlab >=6.5
      mex('-f',engmatopts,'matfile.c','-V5');
    else
      mex('-f',engmatopts,'matfile.c');
    end
    % uncomment the line below if you want to build test_flops.exe
    % This program lets you check the flop counts on your processor.
    % mex('-f',engmatopts,'tests/test_flops.c');
  catch
    disp('Could not install the standalone programs.');
    disp(lasterr)
  end
else
  % UNIX
  eval(['mex' flags 'addflops.c flops.o'])
	if atleast78
		eval(['mex' flags 'gammaln.c util.o -lm -outdir @double'])
	else
		eval(['mex' flags 'gammaln.c util.o -lm'])
	end
  eval(['mex' flags 'digamma.c util.o -lm'])
  eval(['mex' flags 'trigamma.c util.o -lm'])
  eval(['mex' flags 'tetragamma.c util.o -lm'])
	eval(['mex' flags 'setnonzeros.c'])
  if ismac
    % thanks to Nicholas Butko for these mac-specific lines
		clear librandom.dylib randomseed randbinom randgamma sample_hist
		cmd = [options.COMPILER ' ' options.COMPFLAGS ' ' options.OPTIMFLAGS ' -c random.c; ' options.COMPILER ' ' options.COMPFLAGS ' -dynamiclib -Wl,-install_name,`pwd`/librandom.dylib -o librandom.dylib random.o'];
		disp(cmd);
		system(cmd)
    eval(['mex' flags 'randomseed.c util.o librandom.dylib -lm'])
    eval(['mex' flags 'randbinom.c mexutil.o util.o librandom.dylib -lm'])
    eval(['mex' flags 'randgamma.c mexutil.o util.o librandom.dylib -lm'])
    eval(['mex' flags 'sample_hist.c util.o librandom.dylib -lm'])
  else
    % this command only works on linux
		clear librandom.so randomseed randbinom randgamma sample_hist
    cmd = [options.COMPILER ' ' options.COMPFLAGS ' ' options.OPTIMFLAGS ' -c random.c; ' options.COMPILER ' ' options.COMPFLAGS ' -shared -Wl,-E -Wl,-soname,`pwd`/librandom.so -o librandom.so random.o'];
		disp(cmd);
		system(cmd)
    eval(['mex' flags 'randomseed.c util.o librandom.so -lm'])
    eval(['mex' flags 'randbinom.c mexutil.o util.o librandom.so -lm'])
    eval(['mex' flags 'randgamma.c mexutil.o util.o librandom.so -lm'])
    eval(['mex' flags 'sample_hist.c util.o librandom.so -lm'])
  end
  eval(['mex' flags 'repmat.c mexutil.o'])
  try
    % standalone programs
    if atleast78
			disp('lightspeed''s matfile utility is not supported for this version of Matlab');
		elseif atleast65
      % -V5 is required only for Matlab >=6.5
      mex -f matopts.sh matfile.c -V5
    else
      mex -f matopts.sh matfile.c
    end
    % uncomment the line below if you want to build test_flops.exe
    % This program lets you check the flop counts on your processor.
    % mex -f matopts.sh tests/test_flops.c
  catch
    disp('Could not install the standalone programs.');
    disp(lasterr);
    fprintf('Note: if matlab cannot find matopts.sh, your installation of matlab is faulty.\nIf you get this error, don''t worry, lightspeed should still work.');
  end
end

addpath(genpath(pwd))
fprintf('Done.\n');
fprintf('Type "test_lightspeed" to verify the installation.\n');

author: andrewssobral

repo: lrslibrary


text:
%% Example of plotting a 2D Gaussian Mixture Model
clear all;
%% 1) Generate a random GMM

rot_mat = @(theta) [cos(theta), -sin(theta);sin(theta) cos(theta)];

Priors       = [0.3,0.3,0.3];
Priors       = Priors./sum(Priors);

Mu(:,1)      = [-2 0]';
Mu(:,2)      = [2 1]';
Mu(:,3)      = [0 -2]';
Sigma        = [];
Sigma(:,:,1) = rot_mat(pi/5)  * [3,0;0,0.1];
Sigma(:,:,2) = rot_mat(-pi/5) * [3,0;0,0.1];
Sigma(:,:,3) = rot_mat(pi/10) * [3,0;0,0.1];

%% Plot the Contours of the GMM

close all;
figure;
hold on; grid on;
plot_gmm_contour(gca,Priors,Mu,Sigma,[0 0 1]);
box on;

author: epfl-lasa


text:
function grad = B_realImag2complex(future_layers, curr_layer)

future_grad = GetFutureGrad(future_layers, curr_layer);

real_grad = real(future_grad);
imag_grad = -imag(future_grad);
grad = [real_grad; imag_grad];
end


text:
function [p1,p2,dp1,dp2] = autoGen_cartPoleKinematics(x,q,dx,dq,l,empty)
%AUTOGEN_CARTPOLEKINEMATICS
%    [P1,P2,DP1,DP2] = AUTOGEN_CARTPOLEKINEMATICS(X,Q,DX,DQ,L,EMPTY)

%    This function was generated by the Symbolic Math Toolbox version 6.3.
%    30-Nov-2015 20:13:07

p1 = [x;empty];
if nargout > 1
    t2 = cos(q);
    t3 = sin(q);
    p2 = [x+l.*t3;-l.*t2];
end
if nargout > 2
    dp1 = [dx;empty];
end
if nargout > 3
    dp2 = [dx+dq.*l.*t2;dq.*l.*t3];
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [dat, A]= proc_commonAverageReference(dat, refChans, rerefChans)
%PROC_COMMONAVERAGEREFERENCE - rereferencing to a common reference
%
%Synopsis:
% dat= proc_commonAverageReference(dat, <refChans, rerefChans>)
%
% rereference signals to common average reference. you should only
% used scalp electrodes as reference, not e.g. EMG channels.
%
%Arguments:
%      dat        - data structure of continuous or epoched data
%      refChans   - channels used as average reference,
%                   see util_chanind for format,
%                   default util_scalpChannels(dat)
%      rerefChans - those channels are rereferenced, default refChans
%
%Returns:
%      dat        - updated data structure
%
% SEE util_scalpChannels, util_chanind

dat = misc_history(dat);

misc_checkType(dat, 'STRUCT(x clab)');

if ~exist('refChans','var') || isempty(refChans)
  refChans= dat.clab(util_scalpChannels(dat));
end
if ~exist('rerefChans','var') || isempty(rerefChans),
  rerefChans= refChans;
end

misc_checkType(refChans, 'CELL{CHAR}|CHAR');
misc_checkType(rerefChans, 'CELL{CHAR}|CHAR');

rc= util_chanind(dat, refChans);
rrc= util_chanind(dat, rerefChans);
car= mean(dat.x(:,rc,:), 2);

nChans= size(dat.x,2);
A= eye(nChans, nChans);
A(rc,rrc)= A(rc,rrc) - 1/length(rc);
dat= proc_linearDerivation(dat, A, 'CLab','copy', 'Appendix',' car');

author: bbci


text:
function bessel_i1_values_test ( )

%*****************************************************************************80
%
%% BESSEL_I1_VALUES_TEST demonstrates the use of BESSEL_I1_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_I1_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_I1_VALUES stores values of \n' );
  fprintf ( 1, '  the Bessel I1 function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            I1(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_i1_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mtrSetupGastrocParams(subjID)

% Setting up parameters files for a subject based on ROIs assumes
% 1. You are in the subject directory
% 2. A met_params.txt file is already in the conTrack directory

localSubjDir = ['/biac2/wandell2/data/conTrack/muscle/' subjID];
dlhSubjDir = ['/radlab_share/home/tony/images/' subjID];
lappySubjDir = ['c:/cygwin/home/sherbond/images/' subjID];
machineListL = {'deep-lumen1','deep-lumen2','deep-lumen4'};
% machineListR = {'deep-lumen5','deep-lumen6','deep-lumen8'};
mtrTwoRoiSamplerSISR(localSubjDir, 'tendon_plate.mat', 'end_voi.mat', 'met_params.txt', 'GAS_met_params.txt', 'paths_GAS.dat', 'gasRoisMask.nii.gz');
mtrTwoRoiSamplerSISR(localSubjDir, 'tendon_plate.mat', 'end_voi.mat', 'met_params.txt', 'GAS_met_params_lap.txt', 'paths_GAS.dat', 'gasRoisMask.nii.gz', 'temp.sh', lappySubjDir, machineListL);

delete(fullfile('conTrack','temp.sh'));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function V = volt(k, d, varargin)
%VOLT      Volterra integral operator.
%   This function is deprecated. Use OPERATORBLOCK.VOLT.

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

V = linop( operatorBlock.volt(k, double(d), varargin{:}) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
report_this_filefun(mfilename('fullpath'));

figure

[xx,yy]=meshgrid(vlon,vlat);
surfl(yy,xx,tmap/1000),shading interp;

li = light('Position',[ 0 0 100],'Style','infinite');
material shiny
lighting gouraud
% axis([ min(vlat) max(vlat) min(vlon) max(vlon) ]);

set(gca,'FontSize',12,'FontWeight','bold',...
    'LineWidth',1.5,...
    'Box','on','SortMethod','childorder','TickDir','out')
axis ij
view([ -90 90])
colormap(gray)

set(gca,'Ylim',[min(vlon) max(vlon)]);
set(gca,'Xlim',[min(vlat) max(vlat)]);

hold on

%pl = plot3(a.Latitude,a.Longitude,a.Longitude*0+6000,'or');

author: CelsoReyes


text:
addpath('/home/bob/matlab/stats/');
[f,sc] = findSubjects('','_dt6',{'es041113','tk040817'});
N = length(f);
outDir = '/teal/scr1/dti/assymetry';

h=figure(99);
pause(0.5);

for(ii=1:N)
  disp(['Computing assymetry maps for ' sc{ii} '...']);
  roiPath = fullfile(fileparts(f{ii}), 'ROIs');
  dt = load(f{ii});
  % *** TO DO: we could find the optimal mid-sagittal split by
  % coregistering the two halves to minimize error.
  cc(ii) = dtiReadRoi(fullfile(roiPath,'CC_FA'));
  ccMean = mean(cc(ii).coords);
  midSag = inv(dt.xformToAcPc)*[ccMean 1]';
  midSag = round(midSag(1));
  [eVec,eVal] = dtiSplitTensor(dt.dt6);
  pdd = squeeze(eVec(:,:,:,[1 2 3],1));
  clear eVec;
  fa = dtiComputeFA(eVal);
  clear eVal;
  b0 = mrAnatHistogramClip(double(dt.b0),0.4,0.995);

  sz = size(fa);
  nLR = 40;
  [X,Y,Z] = ndgrid([1:midSag-1], [1:sz(2)], [1:sz(3)]);
  leftInds = sub2ind(sz, X(:), Y(:), Z(:));
  [X,Y,Z] = ndgrid([midSag+1:sz(1)], [1:sz(2)], [1:sz(3)]);
  rightInds = sub2ind(sz, X(:), Y(:), Z(:));
  n = length(leftInds)/3;
  for(jj=1:3)
    tmp = pdd(:,:,:,jj);
    lpdd(:,:,:,jj) = reshape(tmp(leftInds), nLR, sz(2), sz(3));
    rpdd(:,:,:,jj) = reshape(tmp(rightInds), nLR, sz(2), sz(3));
  end
  lfa = reshape(fa(leftInds), nLR, sz(2), sz(3));
  rfa = reshape(fa(rightInds), nLR, sz(2), sz(3));
  lb0 = reshape(b0(leftInds), nLR, sz(2), sz(3));
  rb0 = reshape(b0(rightInds), nLR, sz(2), sz(3));

  % Mirror-flip the left-right axis for the righ hemi
  rfa = flipdim(rfa,1);
  rb0 = flipdim(rb0,1);
  rpdd = flipdim(rpdd,1);
  % For Vectors, we might want to flip the vector direction about the midline
  %rpdd(:,:,:,1) = -rpdd(:,:,:,1);
  %rpdd(:,:,:,2) = -rpdd(:,:,:,2);

  mask = zeros(sz);
  delta = lb0>0.1 & lfa>0.05 & rb0>0.1 & rfa>0.05;
  delta = dtiCleanImageMask(delta, 7);
  mask(leftInds) = delta;
  mask(rightInds) = flipdim(delta,1);

  % compute angle difference between left/right PDD
  delta = acos(dot(lpdd, rpdd, 4));
  delta = reshape(delta, nLR, sz(2), sz(3));
  pdd = zeros(sz);
  pdd(leftInds) = delta;
  pdd(rightInds) = flipdim(delta,1);
  pdd(~mask) = 0;
  %figure(h);subplot(3,1,1);
  %imagesc(makeMontage(pdd,[20:60]));colormap(hot);axis image;colorbar;title(sc{ii});

  delta = acos(dot(abs(lpdd), abs(rpdd), 4));
  delta = reshape(delta, nLR, sz(2), sz(3));
  abspdd = zeros(sz);
  abspdd(leftInds) = delta;
  abspdd(rightInds) = flipdim(delta,1);
  abspdd(~mask) = 0;

  fa(:) = 0;
  delta = lfa-rfa;
  fa(leftInds) = delta;
  fa(rightInds) = flipdim(delta,1);

  b0(:) = 0;
  delta = lb0-rb0;
  b0(leftInds) = delta;
  b0(rightInds) = flipdim(delta,1);

  % Spatially normalize the symetry maps
  sn = dt.t1NormParams(2).sn;
  tMm = sqrt(sum(sn.VG(1).mat(1:3,1:3).^2));
  tOrig  = sn.VG.mat\[0 0 0 1]';
  tOrig  = tOrig(1:3)';
  bb = [-tMm .* (tOrig-1) ; tMm.*(sn.VG(1).dim(1:3)-tOrig)];
  og  = -tMm .* tOrig;
  M1  = [tMm(1) 0 0 og(1) ; 0 tMm(2) 0 og(2) ; 0 0 tMm(3) og(3) ; 0 0 0 1];
  outMmPerVox = [2 2 2];
  of  = -outMmPerVox.*(round(-bb(1,:)./outMmPerVox)+1);
  M2  = [outMmPerVox(1) 0 0 of(1) ; 0 outMmPerVox(2) 0 of(2) ; 0 0 outMmPerVox(3) of(3) ; 0 0 0 1];
  d.inMat = inv(sn.VG(1).mat*inv(M1)*M2);
  dField = mrAnatSnToDeformation(sn, dt.mmPerVox, bb);
  d.deformX = dField(:,:,:,1);
  d.deformY = dField(:,:,:,2);
  d.deformZ = dField(:,:,:,3);
  clear dField;
  d.outMat = inv(dt.xformToAcPc);
  [npdd,xf] = mrAnatResliceSpm(pdd, d, bb, dt.mmPerVox, [1 1 1 0 0 0], 0);
  npdd(isnan(npdd)) = 0;
  [nabspdd,xf] = mrAnatResliceSpm(abspdd, d, bb, dt.mmPerVox, [1 1 1 0 0 0], 0);
  nabspdd(isnan(nabspdd)) = 0;
  [nfa,xf] = mrAnatResliceSpm(fa, d, bb, dt.mmPerVox, [1 1 1 0 0 0], 0);
  nfa(isnan(nfa)) = 0;
  [nb0,xf] = mrAnatResliceSpm(b0, d, bb, dt.mmPerVox, [1 1 1 0 0 0], 0);
  nb0(isnan(nb0)) = 0;
  [nmask,xf] = mrAnatResliceSpm(mask, d, bb, dt.mmPerVox, [1 1 1 0 0 0], 0);
  nmask(isnan(nmask)) = 0;
  figure(h);subplot(2,2,1);
  imagesc(makeMontage(npdd,[20:60]));colormap(hot);axis image;colorbar;
  figure(h);subplot(2,2,2);hist(npdd(logical(nmask(:)>=0.5)),100);
  title(sprintf('%0.2f', mean(npdd(logical(nmask(:)>=0.5)))));
  figure(h);subplot(2,2,3);
  imagesc(makeMontage(nabspdd,[20:60]));colormap(hot);axis image;colorbar;
  figure(h);subplot(2,2,4);hist(nabspdd(logical(nmask(:)>=0.5)),100);
  title(sprintf('%0.2f', mean(nabspdd(logical(nmask(:)>=0.5)))));
  refresh(h);;
  save(fullfile(outDir,[sc{ii} '_assymMaps']), 'pdd', 'abspdd', 'fa', 'b0', 'mask', ...
       'npdd', 'nabspdd', 'nfa', 'nb0', 'nmask');
  % Save as analyze?
end

am = load(fullfile(outDir,[sc{1} '_assymMaps']));
sz = size(am.nfa);
pdd = zeros([sz N]);
apdd = zeros([sz N]);
fa = zeros([sz N]);
b0 = zeros([sz N]);
mnMask = ones(sz);
mask = zeros([sz N]);
for(ii=1:N)
  disp(['Loading assymetry maps for ' sc{ii} '...']);
  am = load(fullfile(outDir,[sc{ii} '_assymMaps']));
  % wrap to 0-pi/2, since the PDD is sign-invariant.
  wrapThese = am.npdd>pi/2;
  am.npdd(wrapThese) = am.npdd(wrapThese)-pi/2;
  pdd(:,:,:,ii) = am.npdd;
  apdd(:,:,:,ii) = am.nabspdd;
  b0(:,:,:,ii) = am.nb0;
  fa(:,:,:,ii) = am.nfa;
  %pdd(:,:,:,ii) = dtiSmooth3(am.npdd,3);
  %apdd(:,:,:,ii) = dtiSmooth3(am.nabspdd,3);
  %b0(:,:,:,ii) = dtiSmooth3(am.nb0,3);
  %fa(:,:,:,ii) = dtiSmooth3(am.nfa,3);
  mask(:,:,:,ii) = am.nmask;
  mnMask = mnMask & am.nmask>0.5;
end
mnB0 = mean(b0,4);
sdB0 = std(b0,0,4)+0.00000001;
zB0 = mnB0./sdB0;
mnFa = mean(fa,4);
sdFa = std(fa,0,4)+0.00000001;
zFa = mnFa./sdFa;
mnPdd = mean(pdd,4);
sdPdd = std(pdd,0,4)+0.00000001;
zPdd = mnPdd./sdPdd;
mnAPdd = mean(apdd,4);
sdAPdd = std(apdd,0,4)+0.00000001;
zPdd = mnPdd./sdPdd;
figure(50);imagesc(makeMontage(mnPdd,[20:60]));axis image;colorbar;title('mean PDD');
figure(51);imagesc(makeMontage(sdPdd,[20:60]));axis image;colorbar;title('stdev PDD');
figure(52);imagesc(makeMontage(zPdd,[20:60]));axis image;colorbar;title(['zPDD']);
figure(53);imagesc(makeMontage(zFa,[20:60]));axis image;colorbar;title('zFA');
figure(54);imagesc(makeMontage(zB0,[20:60]));axis image;colorbar;title(['zB0']);

[behData,colNames] = dtiGetBehavioralData(sc);
behIndex = 6;

figure;hist(pdd(mask(:)>=0.5),100);
title(sprintf('Delta PDD (mean=%0.3f)', mean(pdd(mask(:)>=0.5))));

totalAssym = zeros(1,N);
for(ii=1:N)
  tmp = pdd(mask(:,:,:,ii)>=0.5);
  totalAssym(ii) = mean(tmp);
end

% Compute correlations

%symZ = (apdd-repmat(mnAPdd, [1 1 1 N]))./repmat(sdAPdd,[1 1 1 N]);
symZ = (pdd-repmat(mnPdd, [1 1 1 N]))./repmat(sdPdd,[1 1 1 N]);
%symZ = (fa-repmat(mnFa, [1 1 1 N])) ./ repmat(sdFa, [1 1 1 N]);
%symZ = (b0-repmat(mnB0, [1 1 1 N])) ./ repmat(sdB0, [1 1 1 N]);
beh = behData(:,behIndex);
mnBeh = mean(beh);
sdBeh = std(beh);
behZ = (beh-mnBeh) ./ sdBeh;
r = zeros(size(mnPdd));
for(ii=1:N)
    r = r + symZ(:,:,:,ii).*behZ(ii);
end
r = r./(N-1);
r(~mnMask) = 0;

% compute Fischer's z': Z = 0.5*(log(1+r)-log(1-r));
% std err of Fishcer's z is 1/sqrt(n-3)
Z = 0.5*(log((1+r)./(1-r)));
df = N-3;
p = erfc((abs(Z)*sqrt(df))/sqrt(2));
%p = normpdf(Z,0,1/sqrt(df))
pn = -log10(p);
pn(pn>10) = 10;

%figure;hist(Z(mnMask(:)),100);
%title(sprintf('Fischer''s z (mean=%0.3f)', mean(Z(mnMask(:)))));

q = 0.2;
p_sorted = sort(p(mnMask(:)));
numP = length(p_sorted);
I = [1:numP]'./numP*q;
cVID = 1;
pID = p_sorted(max(find(p_sorted<=I/cVID)));
%cVN = sum(1./(1:numP));
%pN = p_sorted(max(find(p_sorted<=I/cVN)));
%H = fdrEmpDistr(Z(:), .01);
%H0 = fdrTheoNull(H, 't', df);
%[thr, fdrCurve] = fdrThresh({'FDR',1}, H, H0, 0.1);
%pThresh = spm_P_FDR(Z(:),[df],'Z',1,p_sorted);
pn(~mnMask) = 0;
figure;imagesc(makeMontage(pn,[20:60]));axis image;colorbar;

% There is an assymetry/PA correlation at image coord [15 41 31]
% (acpc coord [-51 -39 1]) which is on the posterior inferior
% temporal gyrus. Fiber tracking shows this region as clearly part
% of the arcuate, connecting MTG with frontal (language?) regions.

return;

% MTG location:
c = [15 41 31]; % acpc = [-51 -39 1]
% lateral occipital location:
%c = [27 20 36];
% ILF:
c = [60 47 28]; % acpc = [38 -26 -5];
for(behIndex=1:length(colNames))
  assymIndex = squeeze(b0(c(1), c(2), c(3), :));
  gv = ~isnan(behData(:,behIndex)) & ~isnan(assymIndex);
  [s.p,s.r,s.df] = statTest(assymIndex(gv), behData(gv,behIndex), 'r');
  if(s.p<0.0001) sig='***';
  elseif(s.p<0.001) sig='**';
  elseif(s.p<0.01) sig='*';
  else sig=''; end
  msg = sprintf('%s %s (r=%0.2f, p=%0.5f, df=%d)', sig, colNames{behIndex}, s.r, s.p, s.df);
  disp(msg);
end
behIndex = 6;
figure; scatter(assymIndex, behData(:,behIndex));
title(msg);

%txform = eye(4); txform(4,1:3) = [-81 -121 -61]';
txform = dt.xformToAcPc;
upsamp = 1;
imPerRow = 7;
firstIm = 20;
sp=ginput(1)./upsamp;y=mod(sp(1),sz(2));x=mod(sp(2),sz(1));z=floor(sp(2)/sz(1))*imPerRow+floor(sp(1)/sz(2))+firstIm;c=round([x y z]);acpc=mrAnatXformCoords(txform,c);
fprintf(['motage coords = [%0.1f %0.1f]; Image coords = [%0.1f %0.1f ' ...
'%0.1f]; AcPc coords = [%0.1f %0.1f %0.1f]\n'], sp, x, y, z, acpc);

%Threshold probs
p_norm = -log10(p);
p_norm = p_norm./max(p_norm(:));
p_norm = round(p_norm*255+1);


text:
function F = conj( F )
%CONJ Complex conjugate of a DISKFUNV.
%   CONJ(F) returns the complex conjugate of F. For a complex F, CONJ(F) =
%   REAL(F) - i*IMAG(F).
%
%   Since only real DISKFUNV objects are supported, the result of this
%   function is always the original function F.
%
%   See also DISKFUNV/IMAG, DISKFUNV/REAL.

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

% Empty check:
if ( isempty( F ) )
    return
end

% Take conj part of each component:
F.components = cellfun( @conj, F.components, 'UniformOutput', false );

end


text:
function blkStruct = slblocks
%SLBLOCKS Define the Simulink library block representation.
%
%    Define the Simulink library block representation for the Testing

%    CP 10-11-05
%    Copyright 2005-2006, Fraunhofer FOKUS
%    $Revision: 1.0.0.0 $  $Date: 2005/11/10 20:11:50 $
%    $Revision: 1.0.0.0 $  $Date: 2006/07/10 20:11:50 $

blkStruct.Name    = 'MIL Test Blockset';
% The function that is called when
% the user double-clicks on this icon.
blkStruct.OpenFcn = 'MIL_Test';
blkStruct.MaskInitialization = '';

% Define the library list for the Simulink Library browser.
% Return the name of the library model and the name for it
Browser(1).Library = 'MIL_Test';
Browser(1).Name    = 'MIL Test Blockset';
%Browser(1).IsFlat  = 1;% Is this library "flat" (i.e. no subsystems)?

blkStruct.Browser = Browser;

% End of slblocks.m

author: Sable

repo: mcbench-benchmarks


text:
function alist_generate

load H
% H = dvbs2ldpc(1/2);
% H = [0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0; ...
%  0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0; ...
%  0 1 0 0 1 0 1 0 0 1 0 0 0 0 0 0; ...
%  0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0; ...
%  0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 1; ...
%  1 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0; ...
%  0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0; ...
%  0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1; ...
%  1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1; ...
%  0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0; ...
%  0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0; ...
%  1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0];

[N, M] = size(H);

fid = fopen('NR_ldpc.alist', 'w');
% fid = fopen('dvbs2_ldpc_12.alist', 'w');
fprintf(fid, '%d %d\n', N, M);

H = full(H);

row_non_zero_entry_number = sum(H, 2).';
column_non_zero_entry_number = sum(H, 1);

fprintf(fid, '%d %d\n', max(row_non_zero_entry_number), max(column_non_zero_entry_number));

for index = 1:N
  fprintf(fid, '%d ', row_non_zero_entry_number(index));
end

fprintf(fid, '\n');

for index = 1:M
  fprintf(fid, '%d ', column_non_zero_entry_number(index));
end

fprintf(fid, '\n');

for row_index = 1:N
  non_zero_entry_positions = find(H(row_index, :) ~= 0);

  for index = 1:max(row_non_zero_entry_number) % length(non_zero_entry_positions)
    if index <= length(non_zero_entry_positions)
      fprintf(fid, '%d ', non_zero_entry_positions(index));
    else
      fprintf(fid, '0 ');
    end
  end

  fprintf(fid, '\n');
end

for column_index = 1:M
  non_zero_entry_positions = find(H(:, column_index) ~= 0);

  for index = 1:max(column_non_zero_entry_number)
    if index <= length(non_zero_entry_positions)
      fprintf(fid, '%d ', non_zero_entry_positions(index));
    else
      fprintf(fid, '0 ');
    end
  end

  fprintf(fid, '\n');
end

fclose(fid);

end

author: xiaoshaoning


text:
function geometry_test012 ( )

%*****************************************************************************80
%
%% TEST012 tests CIRCLE_LUNE_AREA_2D, CIRCLE_SECTOR_AREA_2D, CIRCLE_TRIANGLE_AREA_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  n_test = 12;

  center(1:2) = [ 0.0, 0.0 ];
  r = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST012\n' );
  fprintf ( 1, '  CIRCLE_LUNE_AREA_2D computes the area of a\n' );
  fprintf ( 1, '    circular lune, defined by joining the endpoints\n' );
  fprintf ( 1, '    of a circular arc.\n' );
  fprintf ( 1, '  CIRCLE_SECTOR_AREA_2D computes the area of a\n' );
  fprintf ( 1, '    circular sector, defined by joining the endpoints\n' );
  fprintf ( 1, '    of a circular arc to the center.\n' );
  fprintf ( 1, '  CIRCLE_TRIANGLE_AREA_2D computes the signed area of a\n' );
  fprintf ( 1, '    triangle, defined by joining the endpoints\n' );
  fprintf ( 1, '    of a circular arc and the center.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '      R            Theta1      Theta2        Sector       Triangle     Lune\n' );
  fprintf ( 1, '\n' );

  for i = 0 : n_test

    theta1 = 0.0;
    theta2 = i * 2.0 * pi / n_test;

    area1 = circle_sector_area_2d ( r, center, theta1, theta2 );

    area2 = circle_triangle_area_2d ( r, center, theta1, theta2 );

    area3 = circle_lune_area_2d ( r, center, theta1, theta2 );

    fprintf ( 1, '  %10f  %10f  %10f  %10f  %10f  %10f\n', ...
      r, theta1, theta2, area1, area2, area3 );

  end

  return
end


text:
function [rI, Par]   =  MCWNNM_ADMM1_Denoising( nI, I, Par )
rI           = nI;   % Estimated Image
[h, w, ch]  = size(rI);
Par.h = h;
Par.w = w;
Par.ch = ch;
Par = SearchNeighborIndex( Par );
% noisy image to patch
NoiPat =	Image2Patch( nI, Par );
Par.TolN = size(NoiPat, 2);
Sigma_arrCh = zeros(Par.ch, Par.TolN);
for iter = 1 : Par.Iter
    Par.iter = iter;
    % iterative regularization
    rI =	rI + Par.delta * (nI - rI);
    % image to patch
    CurPat =	Image2Patch( rI, Par );
    % estimate local noise variance
    for c = 1:Par.ch
        if(iter == 1)
            TempSigma_arrCh = sqrt(max(0, repmat(Par.nSig0(c)^2, 1, size(CurPat, 2)) - mean((NoiPat((c-1)*Par.ps2+1:c*Par.ps2, :) - CurPat((c-1)*Par.ps2+1:c*Par.ps2, :)).^2)));
            %             TempSigma_arrCh = sqrt(abs(repmat(Par.nSig0(c)^2, 1, size(CurPat, 2)) - mean((NoiPat((c-1)*Par.ps2+1:c*Par.ps2, :) - CurPat((c-1)*Par.ps2+1:c*Par.ps2, :)).^2)));
        else
            TempSigma_arrCh = Par.lambda*sqrt(max(0, repmat(Par.nSig0(c)^2, 1, size(CurPat, 2)) - mean((NoiPat((c-1)*Par.ps2+1:c*Par.ps2, :) - CurPat((c-1)*Par.ps2+1:c*Par.ps2, :)).^2)));
            %             TempSigma_arrCh = Par.lambda*sqrt(abs(repmat(Par.nSig0(c)^2, 1, size(CurPat, 2)) - mean((NoiPat((c-1)*Par.ps2+1:c*Par.ps2, :) - CurPat((c-1)*Par.ps2+1:c*Par.ps2, :)).^2)));
        end
        Sigma_arrCh((c-1)*Par.ps2+1:c*Par.ps2, :) = repmat(TempSigma_arrCh, [Par.ps2, 1]);
    end
    if (mod(iter-1, Par.Innerloop) == 0)
        Par.nlsp = Par.nlsp - 10;  % Lower Noise level, less NL patches
        NL_mat  =  Block_Matching(CurPat, Par);% Caculate Non-local similar patches for each
    end
    % Denoising by MCWNNM
    [Y_hat, W_hat]  =  MCWNNM_ADMM1_Estimation( NL_mat, Sigma_arrCh, CurPat, Par );   % Estimate all the patches
    rI = PGs2Image(Y_hat, W_hat, Par);
    PSNR  = csnr( I, rI, 0, 0 );
    SSIM      =  cal_ssim( I, rI, 0, 0 );
    fprintf( 'Iter = %2.3f, PSNR = %2.2f, SSIM = %2.2f \n', iter, PSNR, SSIM );
    Par.PSNR(iter, Par.image)  =   PSNR;
    Par.SSIM(iter, Par.image)      =  SSIM;
end
return;

author: csjunxu


text:
function mesh2off(OFFfilename, faces,vertices)
% output a mesh into an OFF file
% input:
% OFFfilename is the file name for the off file
% faces is Nx3 matrix for vertex index
% vertices is KxD matrix. D is the number of polygon size. =3 if triangles

file = fopen(OFFfilename,'w');
fprintf (file, 'OFF %d %d 0\n', size(vertices,1), size(faces,1));
for v=1:size(vertices,1)
    fprintf(file, '%f %f %f\n', vertices(v,1), vertices(v,2), vertices(v,3));
end

faces = faces -1; % matlab starts from 1, office starts from 0

for f=1:size(faces,1)
    fprintf(file, '%d', size(faces,2));
    for i=1:size(faces,2)
        fprintf(file, ' %f', faces(f,i));
    end
    fprintf(file, '\n');
end

fclose(file);

author: jianxiongxiao

repo: ProfXkit


text:
 function y = Gtomo_nufft_filter(omega, ob, do_phase)
%function y = Gtomo_nufft_filter(omega, ob, do_phase)
% build the sinogram-spectrum-sized matrix
% that is .* multiplied after 2D FT, before iFFT or NUiFFT
% in
%	omega	[M 2]	frequency sample locations (radians)
% out
%	y	[M 2]	filter
% Copyright 2001-1, Jeff Fessler, The University of Michigan
% Extend to fan-beam geometry, 2003-11, Yingying Zhang

%
% effect of image-domain shift
%
y = exp(1i * (omega * ob.nxy_shift(:)));

%
% effect of image basis function (extension to non-rect by S. Matej)
% basis: b(x/Dx,y/Dy) <-FT-> Dx*Dy * B(Dx*u,Dy*v)
%
y = y .* (ob.dx).^2;

if isempty(ob.basis.type) | streq(ob.basis.type, 'pixel')
	% sinc_2 due to square pixels
	if ob.chat, printf('pixel basis, dx=%g', ob.dx), end
	y = y .* sinc(omega(:,1)/(2*pi)) .* sinc(omega(:,2)/(2*pi));

elseif streq(ob.basis.type, 'no')
	if ob.chat, printf('no image basis modeled'), end

elseif streq(ob.basis.type, 'KB')
	if ob.chat, printf('KB basis: J=%g, alpha=%g, m=%g, n=%g', ...
		ob.basis.diam, ob.basis.shape, ob.basis.m, ob.basis.dim); end
	norm = kaiser_bessel_ft(0, ...
		ob.basis.diam, ob.basis.shape, ob.basis.m, ob.basis.dim);
	nd = sqrt(omega(:,1)/(2*pi) .* omega(:,1)/(2*pi) + ...
		omega(:,2)/(2*pi) .* omega(:,2)/(2*pi));
	y = y .* kaiser_bessel_ft(nd, ...
		ob.basis.diam, ob.basis.shape, ob.basis.m, ob.basis.dim) / norm;

elseif streq(ob.basis.type, 'Gauss')
	error('Gaussian basis not yet implemented')

else
	error(sprintf('basis function %s not implemented', ob.basis.type))
end

% form into sinogram shape
K = ob.Krho;
y = reshape(y, K, ob.na);

%
% detector radial response resolution-loss effect in frequency domain
% extension to parallel-beam non-rect case by S. Matej
%

kk = [-K/2 : K/2-1]';
if isempty(ob.beam.type) | streq(ob.beam.type, 'rect')

	% Trick: in fan-beam, the detector blur is not shift-invariant.
	% We approximate by the blur at the center of object.
	% The effective strip_width / ds at the center is the same.

	strip_ray = ob.strip_width / ob.ds;
	blur = sinc(strip_ray * kk / K);
	if ob.chat
		printf('strip integrals, relative beam width=%g', strip_ray)
	end

elseif streq(ob.beam.type, 'line')
	blur = ones(size(kk));
	if ob.chat, printf('line integrals modeled'), end

elseif streq(ob.beam.type, 'KB')
	if ob.chat, printf('KB beam shape: J=%g, alpha=%g, m=%g', ...
		ob.beam.diam, ob.beam.shape, ob.beam.m), end
	rel_bdiam = ob.beam.diam;	% KB_diameter relative to ob.ds
	norm = kaiser_bessel_ft(0, rel_bdiam, ob.beam.shape, ob.beam.m, 1);
	blur = kaiser_bessel_ft(...
		kk/K, rel_bdiam, ob.beam.shape, ob.beam.m, 1) / norm;

elseif streq(ob.beam.type, 'Gauss')
	g_sigma = ob.beam.shape / sqrt(8*log(2));
	blur = exp(-(2*pi*kk/K).^2 .* g_sigma^2/2);

else
	error(sprintf('beam shape %s not implemented', ob.beam.type))

end

y = y .* repmat(blur, [1 ob.na]);	% include blur effect

%
% phase "shift" to effect a half-pixel shift in each row
% corresponding to "offset=0" in tomographic projection.
% build in the post-fft shift too while at it.
%
%if streq(ob.geometry, 'par') % do we need this in fan-beam ???
if do_phase % do we need this in fan-beam ???
	phase = exp(1i*2*pi*(K+ob.is.shift0)/2 * kk / K);
	y = y .* repmat(phase, [1 ob.na]);
	y = y ./ ob.ds;	% see JF tech. report
end

% trick: for parallel case, build in the phase shift due to offset_s
% added 2005-8-24 since it had been omitted previously
if ~ob.is.fan
	phase = exp(-2i*pi * ob.offset_s * kk / K);
	y = y .* repmat(phase, [1 ob.na]);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [Points,CorCoeffVec,F,CF] = extractNcorrResults(handles_ncorr,IMref)
%% function for Extracting the results from Ncorr and calculate correlated image points, correlation coefficients, faces and face colors, in step 2
%[Points,CorCoeffVec,F,CF] = extractNcorrResults(handles_ncorr,IM);
%
% INPUT:
% * handles_ncorr - the Ncorr handle, after running the analysis (up until "Displacement", no need for "Strains").
% * IM - the reference image, for obtaining the grayscale values for the
% faces
%
% OUTPUT:
% * Points - a 2nX1 cell containing 2D image points of the correlated
% points between the 1st image (reference) and all the other images. Points
% which could not be correlated get NaN instead of a value, but the number
% of points is equal for all images, to preserve correspondance.
% * CorCoeffVec - a 2nX1 cell containing vectors the same length as Points,
% with the correlation coefficient for each point (see Ncorr documentation
% for understanding the meaning of the correlation coefficient).
% * F - triangular faces defined on the reference image (regular grid)
% based on Delaunay
% * CF - grayscale values for each face in F, based on the reference image

%%
DICresults = handles_ncorr.data_dic;

% Extract results
Disp=DICresults.displacements;
DispInfo=DICresults.dispinfo;
nCur=size(Disp,2);

Factor=DispInfo.spacing+1;
ROI_DIC=cell(nCur,1);
CorCoeff=cell(nCur,1); CorCoeffVec=cell(nCur,1); Points=cell(nCur,1);
Uvec=cell(nCur,1); Vvec=cell(nCur,1);

for ii=1:nCur

    ROI_DIC{ii}=Disp(ii).roi_dic.mask;

    CorCoeff{ii}=Disp(ii).plot_corrcoef_dic;
    Uref=Disp(ii).plot_u_ref_formatted;
    Vref=Disp(ii).plot_v_ref_formatted;

    [YrefROIVec,XrefROIVec] = find(ROI_DIC{1});
    [YcurROIVec,XcurROIVec] = find(ROI_DIC{ii});

    PtempRef=[XrefROIVec,YrefROIVec];
    PtempRef=(PtempRef-1)*Factor+1; % switch from sapcing to pixels
    Pref=PtempRef;

    CorCoeffVec{ii}=CorCoeff{ii}(ROI_DIC{1});
    CorCoeffVec{ii}(CorCoeffVec{ii}==0)=NaN;

    % displacements from ref to cur
    UrefROIVec=Uref(ROI_DIC{1});
    UrefROIVec(UrefROIVec==0)=NaN;
    VrefROIVec=Vref(ROI_DIC{1});
    VrefROIVec(VrefROIVec==0)=NaN;

    Uvec{ii}=UrefROIVec;
    Vvec{ii}=UrefROIVec;
    % current points
    Points{ii}=[Pref(:,1)+UrefROIVec,Pref(:,2)+VrefROIVec];

    % save face colors for further 3D analysis
    if ii==1
        % pixel colors
        IMrefSmall=IMref(1:Factor:end,1:Factor:end);
        IMrefSmallMasked=IMrefSmall;
        IMrefSmallMasked(~ROI_DIC{1})=[];
        ColorRef=IMrefSmallMasked(:);
    end

end

% Create triangulation
DT = delaunayTriangulation(Pref);
F = DT.ConnectivityList;
V = Pref;

% remove irregular triangles
EdgeLengths = patchEdgeLengths(F,V);
EdgeLengths = [EdgeLengths(1:3:length(EdgeLengths)) EdgeLengths(2:3:length(EdgeLengths)) EdgeLengths(3:3:length(EdgeLengths))];
EdgeLengthsMax=max(EdgeLengths,[],2);

F(EdgeLengthsMax>1.1*sqrt(2)*Factor,:)=[];

% flip direction to have normals pointing out
F=F(:,[1 3 2]);

% face colors (average node colors)
CF=mean(ColorRef(F),2);

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function [parameters, ll, ht, intercept, VCV, scores, diagnostics] = scalar_vt_vech(data,dataAsym,p,o,q,composite,startingvals,options)
% Estimation of symmetric and asymmetric scalar multivariate vech ARCH models using variance
% targeting to reduce the number of parameters needing to be estimated simultaneously
%
% USAGE:
%  [PARAMETERS,LL,HT,INTERCEPT,VCV,SCORES,DIAGNOSTICS] = scalar_vt_vech(DATA,DATAASYM,P,O,Q,COMPOSITE,STARTINGVALS,OPTIONS)
%
% INPUTS:
%   DATA         - A T by K matrix of zero mean residuals -OR-
%                    K by K by T array of covariance estimators (e.g. realized covariance)
%   DATAASYM     - [OPTIONAL] K by K by T array of asymmetric covariance estimators
%   P            - Positive, scalar integer representing the number of lags of the innovation process
%   O            - Non-negative scalar integer representing the number of asymmetric lags to include
%   Q            - Non-negative scalar integer representing the number of lags of conditional covariance
%   COMPOSUTE    - [OPTIONAL] String, one of:
%                    'None' - Standard K-dimensional QMLE (Default)
%                    'Diagonal' - Use only pairwise likelihoode for i,i+1, i=1,2,...,K-1
%                    'Full' - Use all pairwise likelihoods
%   STARTINGVALS - [OPTIONAL] (p+o+q) x 1 vector of starting values
%   OPTIONS      - [OPTIONAL] Options to use in the optimization (fminunc)
%
% OUTPUTS:
%   PARAMETERS   - A p+o+q column vector of parameters.  The intercept is reported in DIAGNOSTICS
%   LL           - The log likelihood at the optimum
%   HT           - A [K K T] dimension matrix of conditional covariances
%   INTERCEPT    - K by K matrix containing the intercept computed from the unconditional variance
%                    and parameters
%   VCV          - A numParams^2 square matrix of robust parameter covariances (A^(-1)*B*A^(-1)*t^(-1))
%   SCORES       - A T by numParams matrix of individual scores
%
% COMMENTS:
%    The conditional variance, H(t), of a scalar variance-targeting vech is modeled
%    as follows:
%
%      H(t) = (1-alpha(1)+...+alpha(p)-beta(1)-...-beta(q))*C + ...
%              alpha(1)*r_{t-1}'*r_{t-1} + ... + alpha(p)*r_{t-p}'*r_{t-p}+...
%              gamma(1)*n_{t-1}'*n_{t-1} + ... + gamma(o)*n_{t-p}'*n_{t-p}+...
%              beta(1)*H(t-1) +...+ beta(q)*H(t-q)
%
%    where n_{t} = r_{t} .* (r_{t}<0)
%
% EXAMPLES:
%     Estimation of a Scalar VECH(1,0,1)
%       [simulatedData, Ht, pseudoRC] = scalar_vt_vech_simulate(1000, [.02 .04 .95], .01*eye(2), 1, 1, 1, 72);
%       parameters = scalar_vt_vech(simulatedData,[],1,0,1)
%     Estimation of an asymmetric Scalar VECH(1,1,1)
%       parameters = scalar_vt_vech(simulatedData,[],1,1,1)
%     Estimation of an asymmetric Scalar VECH(1,1,1) using realized-type data
%       asymPseudoRC = zeros(size(pseudoRC));
%       for i = 1:1000
%         asymPseudoRC(:,:,i) = pseudoRC(:,:,i).*(double(simulatedData(i,:)<0)'*double(simulatedData(i,:)<0));
%       end
%       parameters = scalar_vt_vech(pseudoRC,asymPseudoRC,1,1,1);

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 10/28/2009

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Argument Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 3
        o=0;
        q=0;
        composite = 'None';
        startingvals=[];
        options = [];
    case 4
        q=0;
        composite = 'None';
        startingvals=[];
        options = [];
    case 5
        composite = 'None';
        startingvals=[];
        options = [];
    case 6
        startingvals=[];
        options = [];
    case 7
        options = [];
    case 8
        % Nothing
    otherwise
        error('Between 3 and 8 arguments required.')
end

%data should be TxK, T>K
if ndims(data)==2
    [t,k]=size(data);
    if ~isempty(dataAsym)
        error('If DATA is a T by K matrix, DATAASYM must be empty.');
    end
    temp = zeros(k,k,t);
    dataAsym = zeros(k,k,t);
    for i=1:t
        temp(:,:,i) = data(i,:)'*data(i,:);
        dataAsym(:,:,i) = (data(i,:).*(data(i,:)<0))'*(data(i,:).*(data(i,:)<0));
    end
    data = temp;
elseif ndims(data)==3
    [k,m,t] = size(data);
    if m~=k
        error('DATA must be K by K by T is a 3D array.');
    end
    if ~isempty(dataAsym)
        if ndims(dataAsym)~=3
            error('DATAASYM must be a 3D array with the same dimensions as DATA');
        end
        [k2,m2,t2]=size(dataAsym);
        if any([k m t]~=[k2 m2 t2])
            error('DATAASYM must be a 3D array with the same dimensions as DATA');
        end
    end
end

if min(t,k)<2 || t<k
    error('DATA must be a T by K matrix or a K by K by T 3D array, T>K>1');
end

%p, o, q much be non-negative scalars
if length(p)>1 || any(p<1) || floor(p)~=p
    error('P must be a positive scalar');
end
if isempty(o)
    o = 0;
end
if length(o)>1 || any(o<0) || floor(o)~=o
    error('O must be a non-negative scalar');
end
if o>0 && isempty(dataAsym)
    error('DATAASYM must be non-empty if O>0.')
end

if isempty(q)
    q = 0;
end
if length(q)>1 || any(q<0) || floor(q)~=q
    error('Q must be a non-negative scalar');
end

if isempty(composite)
    composite = 'None';
end
switch lower(composite)
    case {'none'}
        useComposite = 0;
    case {'diagonal'}
        useComposite = 1;
    case {'full'}
        useComposite = 2;
end

% Startingvals must have p+o+q parameters sum(alpha)+kappa*sum(gamma)+sum(beta)<1
C = mean(data,3);
if o>0
    Casym = mean(dataAsym,3);
    kappa = 1 / (max(eig(C^(-0.5)*Casym*C^(-0.5))) + eps);
else
    Casym = zeros(k);
    kappa = 2;
end
if  ~isempty(startingvals)
    if size(startingvals,2)>size(startingvals,1)
        startingvals = startingvals';
    end
    if length(startingvals)<(p+o+q)
        error('STARTINGVALS should be a P+O+Q by 1 vector');
    end
    %Only validate if provided
    A=startingvals(1:p);
    G=startingvals(p+1:p+o);
    B=startingvals(p+o+1:p+o+q);
    if (sum(A)+sum(G)/kappa+sum(B))>=.999998
        error('Weighted sum of STATINGVALUES must be less than 1. See Comments.');
    end
    if any(A<0) || any(B<0) || any(G<0)
        error('STARTINGVALS must all be nonnegative.');
    end
end

%Make sure options is a valid option structure
if isempty(options)
    options=optimset('fminunc');
    options.Display='iter';
    options.Diagnostics='on';
    options.LargeScale='off';
end
try
    optimset(options);
catch ME
    error('OPTIONS is not a valid options structure');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Argument Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the backCast
backCast = zeros(k);
backCastAsym = zeros(k);
tau = max(ceil(sqrt(t)),k);
weights = .06 * .94.^(0:tau);
weights = weights / sum(weights);
for i=1:tau
    backCast = backCast + weights(i) * data(:,:,i);
    if o>0
        backCastAsym = backCastAsym + weights(i) * dataAsym(:,:,i);
    end
end

%Augment the data with backcasts
[startingvals,lls,output_parameters] = scalar_vt_vech_starting_values(startingvals,data,dataAsym,p,o,q,C,Casym,kappa,useComposite,backCast,backCastAsym); %#ok<NASGU>

% finally to transform the parameters to the unrestricted equivalents
startingvals=scalar_vt_vech_transform(startingvals,p,o,q,kappa);

[parameters,ll,exitflag,output]=fminunc('scalar_vt_vech_likelihood',startingvals,options,data,dataAsym,p,o,q,C,Casym,kappa,backCast,backCastAsym,false,useComposite,true);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimation Robustification
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if exitflag<=0 &&  ll<lls(1)
    %Did not converge, but function improved
    options.MaxFunEvals=4*100*(p+q);
    options.MaxIter=2*100*(p+q);
    parameters=fminunc('scalar_vt_vech_likelihood',parameters,options,data,dataAsym,p,o,q,C,Casym,kappa,backCast,backCastAsym,false,useComposite,true);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimation Robustification
%%%%%%%%%%%%%%%%%%%%%%%%%%%

parameters=scalar_vt_vech_itransform(parameters,p,o,q,kappa);
if nargout>1
    [ll,~,ht]=scalar_vt_vech_likelihood(parameters,data,dataAsym,p,o,q,C,Casym,kappa,backCast,backCastAsym,false,useComposite,false);
    ll=-ll;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the VCV
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargout>4
    k2 = k*(k+1)/2;
    CVech=vech(C);
    if o==0
        momentCount = k2;
        CasymVech = [];
    else
        momentCount = 2*k2;
        CasymVech=vech(Casym);
    end

    A = zeros(momentCount + p + o + q);
    A(1:momentCount, 1:momentCount) = - t * eye(momentCount);
    jointParameters=[CVech;CasymVech;parameters];

    % A will be k2 + p + q square
    scores = zeros(t,momentCount+p+q);
    scoreCount = 1;
    for i=1:k
        for j=i:k
            scores(:,scoreCount) = squeeze(data(j,i,:));
            if o>0
                scores(:,k2+scoreCount) = squeeze(dataAsym(j,i,:));
            end
            scoreCount = scoreCount  + 1;
        end
    end
    [~,gt] = gradient_2sided(@scalar_vt_vech_likelihood,parameters,data,dataAsym,p,o,q,C,Casym,kappa,backCast,backCastAsym,false,useComposite,false);
    scores(:,momentCount+1:momentCount+p+o+q) = gt;
    A(momentCount+1:momentCount+p+o+q,:) = hessian_2sided_nrows(@scalar_vt_vech_likelihood,jointParameters,p+o+q,data,dataAsym,p,o,q,C,Casym,kappa,backCast,backCastAsym,true,useComposite,false);

    A = A/t;
    B = covnw(scores);
    Ainv = inv(A);
    VCV = Ainv*B*Ainv'/t; %#ok<MINV>
    VCV = VCV(momentCount+1:momentCount+p+o+q,momentCount+1:momentCount+p+o+q);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the VCV
%%%%%%%%%%%%%%%%%%%%%%%%%%%
diagnostics.EXITFLAG=exitflag;
diagnostics.ITERATIONS=output.iterations;
diagnostics.FUNCCOUNT=output.funcCount;
diagnostics.MESSAGE=output.message;
diagnostics.kappa=kappa;
alpha = sum(parameters(1:p));
gamma = sum(parameters(p+1:p+o));
beta = sum(parameters(p+o+1:p+o+q));
diagnostics.intercept = C*(1-alpha-beta);
diagnostics.composite = composite;
if o>0
    diagnostics.intercept  = diagnostics.intercept  - gamma*Casym;
end
intercept = diagnostics.intercept;

author: bashtage

repo: mfe-toolbox


text:
function [ fx, fy ] = f11_f1 ( n, x, y )

%*****************************************************************************80
%
%% F11_F1 returns first derivatives of function f11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real FX(N,1), FY(N,1), the derivative values.
%
  fx(1:n,1) = y(1:n,1) + 1.0;
  fy(1:n,1) = x(1:n,1);

  return
end


text:
function [iout] = center_crop(im,a,b)

[w,h,~] = size(im);
c1 = w-a-(w-a)/2;
c2 = h-b-(h-b)/2;
iout = im(c1+1:c1+a,c2+1:c2+b,:);

end

author: cszn

repo: SRMD

save_path: github-repos/MATLAB/cszn-SRMD


text:
function B = multitransp(A, unused) %#ok<INUSD>
% Transpose the matrix slices of an N-D array (no complex conjugate)
%
% function B = multitransp(A)
%
% If A is a 3-D array, then B is a 3-D array such that
%
%     B(:, :, i) = A(:, :, i).'
%
% for each i. If A is an N-D array, then B is an N-D array with the slices
% A(:, :, i, j, k, ...) transposed.
%
% This function is just a wrapper for pagetranspose, with a fallback call
% to multitransp_legacy in case pagetranspose is not available.
% If pagetranspose is available, it is better to call it directly.
% Note that pagemtimes also allows to compute products with transposes
% without explicitly transposing arrays.
%
% See also: multiprod multihconj multiscale multiskew multiskewh multitrace

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Aug. 12, 2021.
% Contributors: Xiaowen Jiang
% Change log:
%
%   Aug. 12, 2021 (NB):
%       Matlab R2020b introduced a built-in function pagetranspose which
%       does essentially everything we ever needed to do with multitransp
%       in Manopt. Accordingly, multitransp became a wrapper for
%       pagetranspose, and the old code for multitransp remains available
%       as multitransp_legacy.

    assert(nargin == 1, ...
           'The new multitransp only takes one input. Check multitransp_legacy.');

    if exist('pagetranspose', 'file') % Added to Matlab R2020b
        B = pagetranspose(A);
    else
    %   warning('manopt:multi', ...
    %          ['Matlab R2020b introduced pagetranspose.\n' ...
    %           'Calling the old code multitransp_legacy instead.\n' ...
    %           'To disable this warning: warning(''off'', ''manopt:multi'')']);
        B = multitransp_legacy(A);
    end

end

author: NicolasBoumal

repo: manopt


text:
function [dynpcm2] = Pa2dynpcm2(Pa)
% Convert units of pressure from pascals to dynes per centimeter squared.
% Chad A Greene 2012
dynpcm2 = Pa*10;


text:
close all; clear all; clc;

%  image = imread('Resources/images/3096.jpg');
 image = imread('Resources/images/28075.jpg');
%   image = imread('Resources/images/113016.jpg');

 %image = imread('Resources/images/3096.jpg');
%     image = rgb2gray(image);

% neighborhoodType = 4; % 4-point connectivity
neighborhoodType = 4; % 8-point connectivty

tic;
[ segmentedImage, binaryImage, regionMatrix ] = RegionGrowingSegmentation(image, neighborhoodType);
display('Segmentation time:');
toc;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function value = r4_binom ( n, m )

%*****************************************************************************80
%
%% R4_BINOM evaluates the binomial coefficient using R4 arithmetic.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, integer N, M, the arguments.
%
%    Output, real VALUE, the binomial coefficient.
%
  persistent sq2pil
  persistent bilnmx
  persistent fintmx

  sq2pil = 0.91893853320467274;

  if ( isempty ( bilnmx ) )
    bilnmx = log ( r4_mach ( 2 ) );
    fintmx = 0.9 / r4_mach ( 3 );
  end

  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_BINOM - Fatal error!\n' );
    fprintf ( 1, '  N < 0.\n' );
    error ( 'R4_BINOM - Fatal error!' )
  end

  if ( m < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_BINOM - Fatal error!\n' );
    fprintf ( 1, '  M < 0.\n' );
    error ( 'R4_BINOM - Fatal error!' )
  end

  if ( n < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_BINOM - Fatal error!\n' );
    fprintf ( 1, '  N < M.\n' );
    error ( 'R4_BINOM - Fatal error!' )
  end

  k = min ( m, n - m );

  if ( k <= 20 && k * log ( max ( n, 1 ) ) <= bilnmx )

    value = 1.0;

    for i = 1 : k
      value = value * ( n - i + 1 ) / i;
    end

  else

    if ( k < 9 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R4_BINOM - Fatal error!\n' );
      fprintf ( 1, '  Result overflows.\n' );
      fprintf ( 1, '  N or M is too big.\n' );
      error ( 'R4_BINOM - Fatal error!' )
    end

    xn = n + 1;
    xk = k + 1;
    xnk = n - k + 1;

    corr = r4_lgmc ( xn ) - r4_lgmc ( xk ) - r4_lgmc ( xnk );

    value = xk * log ( xnk / xk ) ...
      - xn * r4_lnrel ( - ( xk - 1.0 ) / xn ) ...
      - 0.5 * log ( xn * xnk / xk ) + 1.0 - sq2pil + corr;

    if ( bilnmx < value )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R4_BINOM - Fatal error!\n' );
      fprintf ( 1, '  Result overflows.\n' );
      fprintf ( 1, '  N or M is too big.\n' );
      error ( 'R4_BINOM - Fatal error!' )
    end

    value = exp ( value );

  end

  if ( value < fintmx )
    value = r4_aint ( value + 0.5 );
  end

  return
end


text:
function varargout = process_import_surfaces( varargin )
% PROCESS_IMPORT_SURFACES: Import three surfaces in the database.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Import surfaces';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import anatomy'};
    sProcess.Index       = 3;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ImportAnatomy#Import_the_anatomy';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    sProcess.isSeparator = 1;
    % File selection options
    SelectOptions = {...
        '', ...                                % Filename
        '', ...                                % FileFormat
        'open', ...                            % Dialog type: {open,save}
        'Import surface file...', ...          % Window title
        'ImportAnat', ...                      % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                          % Selection mode: {single,multiple}
        'files', ...                           % Selection mode: {files,dirs,files_and_dirs}
        bst_get('FileFilters', 'surface'), ... % Get all the available file formats
        'SurfaceIn'};                          % DefaultFormats: {ChannelIn,DataIn,DipolesIn,EventsIn,MriIn,NoiseCovIn,ResultsIn,SspIn,SurfaceIn,TimefreqIn
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = 'NewSubject';
    % Option: Head surface
    sProcess.options.headfile.Comment = 'Head surface:';
    sProcess.options.headfile.Type    = 'filename';
    sProcess.options.headfile.Value   = SelectOptions;
    % Option: Hemisphere 1
    sProcess.options.cortexfile1.Comment = 'Hemisphere 1:';
    sProcess.options.cortexfile1.Type    = 'filename';
    sProcess.options.cortexfile1.Value   = SelectOptions;
    % Option: Hemisphere 2
    sProcess.options.cortexfile2.Comment = 'Hemisphere 2:';
    sProcess.options.cortexfile2.Type    = 'filename';
    sProcess.options.cortexfile2.Value   = SelectOptions;
    % Option: Inner skull
    sProcess.options.innerfile.Comment = 'Inner skull:';
    sProcess.options.innerfile.Type    = 'filename';
    sProcess.options.innerfile.Value   = SelectOptions;
    % Option: Outer skull
    sProcess.options.outerfile.Comment = 'Outer skull:';
    sProcess.options.outerfile.Type    = 'filename';
    sProcess.options.outerfile.Value   = SelectOptions;
    % Option: Number of vertices
    sProcess.options.nverthead.Comment = 'Number of vertices (head): ';
    sProcess.options.nverthead.Type    = 'value';
    sProcess.options.nverthead.Value   = {7000, '', 0};
    % Option: Number of vertices
    sProcess.options.nvertcortex.Comment = 'Number of vertices (cortex): ';
    sProcess.options.nvertcortex.Type    = 'value';
    sProcess.options.nvertcortex.Value   = {15000, '', 0};
    % Option: Number of vertices
    sProcess.options.nvertskull.Comment = 'Number of vertices (skull): ';
    sProcess.options.nvertskull.Type    = 'value';
    sProcess.options.nvertskull.Value   = {7000, '', 0};
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return
    end
    % Get filenames to import
    HeadFile      = sProcess.options.headfile.Value{1};
    HeadFormat    = sProcess.options.headfile.Value{2};
    CortexFile1   = sProcess.options.cortexfile1.Value{1};
    CortexFormat1 = sProcess.options.cortexfile1.Value{2};
    CortexFile2   = sProcess.options.cortexfile2.Value{1};
    CortexFormat2 = sProcess.options.cortexfile2.Value{2};
    InnerFile     = sProcess.options.innerfile.Value{1};
    InnerFormat   = sProcess.options.innerfile.Value{2};
    OuterFile     = sProcess.options.outerfile.Value{1};
    OuterFormat   = sProcess.options.outerfile.Value{2};
    if isempty([HeadFile, CortexFile1, InnerFile, OuterFile])
        bst_report('Error', sProcess, [], 'Not enough files selected.');
        return
    end
    % Number of vertices
    nVertHead   = sProcess.options.nverthead.Value{1};
    nVertCortex = sProcess.options.nvertcortex.Value{1};
    nVertSkull  = sProcess.options.nvertskull.Value{1};
    if isempty(nVertCortex) || (nVertCortex < 50) || isempty(nVertHead) || (nVertHead < 50) || isempty(nVertSkull) || (nVertSkull < 50)
        bst_report('Error', sProcess, [], 'Invalid number of vertices.');
        return
    end
    % Divide by number of cortex files
    if ~isempty(CortexFile2)
        nVertCortex = nVertCortex / 2;
    end

    % ===== GET/CREATE SUBJECT =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    % Create subject is it does not exist yet
%     if isempty(sSubject)
%         [sSubject, iSubject] = db_add_subject(SubjectName);
%     end
    if isempty(iSubject)
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" does not exist.']);
        return
    end
    % Check if a MRI is available for the subject
    if isempty(sSubject.Anatomy)
        bst_report('Error', sProcess, [], ['No MRI available for subject "' SubjectName '".']);
        return
    end
    % The subject can't be using the default anatomy
    if (iSubject ~= 0) && sSubject.UseDefaultAnat
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" is using the default anatomy (read-only).']);
        return
    end

    % ===== HEAD FILE =====
    if ~isempty(HeadFile)
        % Import file
        [iScalp, OldHeadFile] = import_surfaces(iSubject, HeadFile, HeadFormat, 0);
        OldHeadFile = OldHeadFile{1};
        % Set the file type
        OldHeadFile = db_surface_type(OldHeadFile, 'Scalp');
        % Downsample
        NewHeadFile = tess_downsize(OldHeadFile, nVertHead, 'reducepatch');
        % Delete intial file
        if ~file_compare(OldHeadFile, NewHeadFile)
            file_delete(file_fullpath(OldHeadFile), 1);
            NewHeadFile = file_fullpath(NewHeadFile);
        end
        % Update Comment field
        HeadMat.Comment = 'Head';
        bst_save(file_fullpath(NewHeadFile), HeadMat, 'v7', 1);
    end

    % ===== CORTEX FILE =====
    % First hemisphere
    if ~isempty(CortexFile1)
        % Import file
        [iCortex1, OldCortexFile1] = import_surfaces(iSubject, CortexFile1, CortexFormat1, 0);
        OldCortexFile1 = OldCortexFile1{1};
        % Downsample
        NewCortexFile1 = tess_downsize(OldCortexFile1, nVertCortex, 'reducepatch');
        % Delete intial file
        if ~file_compare(NewCortexFile1, OldCortexFile1)
            file_delete(file_fullpath(OldCortexFile1), 1);
            NewCortexFile1 = file_fullpath(NewCortexFile1);
        end
    end
    % Second hemisphere
    if ~isempty(CortexFile2)
        % Import file
        [iCortex2, OldCortexFile2] = import_surfaces(iSubject, CortexFile2, CortexFormat2, 0);
        OldCortexFile2 = OldCortexFile2{1};
        % Downsample
        NewCortexFile2 = tess_downsize(OldCortexFile2, nVertCortex, 'reducepatch');
        % Delete intial file
        if ~file_compare(NewCortexFile2, OldCortexFile2)
            file_delete(file_fullpath(OldCortexFile2), 1);
            NewCortexFile2 = file_fullpath(NewCortexFile2);
        end
    end
    % Merge hemispheres
    if ~isempty(CortexFile1) && ~isempty(CortexFile2)
        CortexFile = tess_concatenate({NewCortexFile1, NewCortexFile2}, 'cortex');
        % Delete separate hemispheres
        file_delete({file_fullpath(NewCortexFile1), file_fullpath(NewCortexFile2)}, 1);
    % Only one hemisphere
    elseif ~isempty(CortexFile1)
        CortexFile = NewCortexFile1;
        % Update Comment field for Head file
        CortexMat.Comment = 'Cortex';
        bst_save(CortexFile, CortexMat, 'v7', 1);
    end
    % Set file type
    if ~isempty(CortexFile1) || ~isempty(CortexFile2)
        CortexFile = db_surface_type(CortexFile, 'Cortex');
    end

    % ===== INNER SKULL =====
    if ~isempty(InnerFile)
        % Import file
        [iInner, OldInnerFile] = import_surfaces(iSubject, InnerFile, InnerFormat, 0);
        OldInnerFile = OldInnerFile{1};
        % Set the file type
        OldInnerFile = db_surface_type(OldInnerFile, 'InnerSkull');
        % Downsample
        NewInnerFile = tess_downsize(OldInnerFile, nVertSkull, 'reducepatch');
        % Update Comment field
        InnerMat.Comment = 'Inner skull';
        bst_save(file_fullpath(NewInnerFile), InnerMat, 'v7', 1);
        % Delete intial file
        if ~file_compare(OldInnerFile, NewInnerFile)
            file_delete(file_fullpath(OldInnerFile), 1);
        end
    end

    % ===== OUTER SKULL =====
    if ~isempty(OuterFile)
        % Import file
        [iOuter, OldOuterFile] = import_surfaces(iSubject, OuterFile, OuterFormat, 0);
        OldOuterFile = OldOuterFile{1};
        % Set the file type
        OldOuterFile = db_surface_type(OldOuterFile, 'OuterSkull');
        % Downsample
        NewOuterFile = tess_downsize(OldOuterFile, nVertSkull, 'reducepatch');
        % Update Comment field
        OuterMat.Comment = 'Outer skull';
        bst_save(file_fullpath(NewOuterFile), OuterMat, 'v7', 1);
        % Delete intial file
        if ~file_compare(OldOuterFile, NewOuterFile)
            file_delete(file_fullpath(OldOuterFile), 1);
        end
    end

    % Reload subject folder (necessary, because we deleted manually some files,
    % without removing their references in the database)
    db_reload_subjects(iSubject);

    OutputFiles = {'import'};
end

author: brainstorm-tools


text:
function P = tapas_uniqc_spm_imatrix(M)
% Return the parameters for creating an affine transformation matrix
% copy of spm_imatrix from Version 6906 (SPM12) 20-Oct-16
% FORMAT P = spm_imatrix(M)
% M   - Affine transformation matrix
% P   - Parameters (see spm_matrix for definitions)
%__________________________________________________________________________
%
% See also: spm_matrix.m
%__________________________________________________________________________
% Copyright (C) 1996-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner & Stefan Kiebel

%-Translations and Zooms
%--------------------------------------------------------------------------
R         = M(1:3,1:3);
C         = chol(R'*R);
P         = [M(1:3,4)' 0 0 0  diag(C)'  0 0 0];
if det(R)<0, P(7)=-P(7); end % Fix for -ve determinants

%-Shears
%--------------------------------------------------------------------------
C         = diag(diag(C))\C;
P(10:12)  = C([4 7 8]);
R0        = spm_matrix([0 0 0  0 0 0 P(7:12)]);
R0        = R0(1:3,1:3);
R1        = R/R0;

%-This just leaves rotations in matrix R1
%--------------------------------------------------------------------------
%[          c5*c6,           c5*s6, s5]
%[-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
%[-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]

% There may be slight rounding errors making x>1 or x<-1.
rang      = @(x) min(max(x, -1), 1);

P(5)      = asin(rang(R1(1,3)));
if (abs(P(5))-pi/2)^2 < 1e-9
    P(4)  = 0;
    P(6)  = atan2(-rang(R1(2,1)), rang(-R1(3,1)/R1(1,3)));
else
    c     = cos(P(5));
    P(4)  = atan2(rang(R1(2,3)/c), rang(R1(3,3)/c));
    P(6)  = atan2(rang(R1(1,2)/c), rang(R1(1,1)/c));
end

author: translationalneuromodeling

repo: tapas