text:
function [res,qualMeasOut]=MLEM(proj,geo,angles,niter,varargin)
%MLEM solves the tomographic problem by using Maximum Likelihood Expectation
% Maximisation algorithm.
%
%   MLEM(PROJ,GEO,ALPHA,NITER,opt) solves the reconstruction problem
%   using the projection data PROJ taken over ALPHA angles, corresponding
%   to the geometry described in GEO, using NITER iterations.
%
% 'verbose': Get feedback or not. Default: 1
%
% 'init':    Describes different initialization techniques.
%             •  'none'     : Initializes the image to ones (default)
%             •  'FDK'      : Initializes image to FDK reconstruction
%
% 'QualMeas':    Asks the algorithm for a set of quality measurement
%                parameters. Input should contain a cell array of desired
%                quality measurement names. Example: {'CC','RMSE','MSSIM'}
%                These will be computed in each iteration.
%
% 'groundTruth':  An image as ground truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
[verbose,res,QualMeasOpts,gpuids,gt]=parse_inputs(proj,geo,angles,varargin);
measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    res_prev=gt;
    clear gt
end
if nargout<2 && measurequality
    warning("Image metrics requested but none caught as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),niter);

res = max(res,0);

% Back-projection weight, V
V = Atb(ones(size(proj),'single'),geo,angles,'matched','gpuids',gpuids);
V(V<=0.) = inf;

for ii=1:niter
    if measurequality && ~strcmp(QualMeasOpts,'error_norm')
        res_prev = res; % only store if necessary
    end
    if (ii==1);tic;end

    den = Ax(res,geo,angles,'gpuids',gpuids);
    den(den<=0.)=inf;

    imgupdate = Atb(proj./den, geo,angles,'matched','gpuids',gpuids)./V;
    res = max(res.*imgupdate,0.);

    if measurequality
        qualMeasOut(:,ii)=Measure_Quality(res_prev,res,QualMeasOpts);
    end

    if (ii==1)&&(verbose==1)
        expected_time=(toc)*niter;
        disp('MLEM');
        disp(['Expected duration   :    ',secs2hms(expected_time)]);
        disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
        disp('');
    end

end
end

%% Parse inputs
function [verbose,f0,QualMeasOpts,gpuids,gt]=parse_inputs(proj,geo,angles,argin)
opts = {'verbose','init','qualmeas','gpuids','groundtruth'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:MLEM:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('TIGRE:MLEM:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option is not default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('TIGRE:MLEM:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end

    switch opt
        % % % % % % % Verbose
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
        % Initial image
        %  =========================================================================
        case 'init'
            if default || strcmp(val,'none')
                f0=ones(geo.nVoxel','single');
            else
                if strcmp(val,'FDK')
                    f0=FDK(proj, geo, angles);
                else
                    error('TIGRE:MLEM:InvalidInput','Invalid init')
                end
            end
        % Image Quality Measure
        %  =========================================================================
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:MLEM:InvalidInput','Invalid quality measurement parameters');
                end
            end
        % GPUIDS
        % =========================================================================
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        otherwise
            error('TIGRE:MLEM:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in MLEM()']);
    end
end
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% demo with operators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Please adjust pp_code/postprocessing.m with your custom functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all

filename = fullfile('..','images','a0280-IMG_0854.JPG');

device = 'gpu';

load(fullfile('models','model_sRGB-XYZ-sRGB.mat'));
           %|res1|I-res1|XYZ|RGB|res2|RGB+res2
pp_method = 'none|none|denoise+deblur|none|none';
%post-processing methods: denoise, deblur, dehaze, editdetails,
% exposure-fusion, transfer-colors, chrom-adapt, super-res
%The order is: localLayer|sRGB-localLayer|CIE XYZ|sRGB|localLayer

show = 1;

save_output = 1;

opt = [];  %extra options, if needed

if save_output == 1
    out_dir = fullfile('..','results');
    if exist(out_dir, 'dir') == 0
        mkdir(out_dir);
    end
end

image = im2double(imread(filename));

if length(size(image)) ~= 3
    error('cannot deal with grayscale images');
end

fprintf('processing image %s...\n', filename);

[~,name,~] = fileparts(filename);

tasks = strsplit(pp_method,'|');

output_XYZ = applyLocalMapping(nets.local_sRGB, image, ...
    'to-xyz', device, tasks{1}, opt);

output_XYZ = applyGlobalMapping(nets.global_sRGB, output_XYZ,device, ...
    tasks{2}, opt);

if strcmpi(tasks{3},'none') == 0
    output_XYZ = postprocessing(output_XYZ,tasks{3},opt);
end

output_sRGB = applyGlobalMapping(nets.global_XYZ, output_XYZ,device, ...
    tasks{4}, opt);

output_sRGB = applyLocalMapping(nets.local_XYZ, output_sRGB, ...
    'to-srgb', device, tasks{5}, opt);

output_sRGB(output_sRGB>1) = 1;
output_sRGB(output_sRGB<0) = 0;

if save_output == 1
    imwrite(output_sRGB,fullfile(out_dir, [name '_result.png']));
end

if show == 1
    subplot(1,2,1);imshow(image);  title('input');
    subplot(1,2,2);imshow(output_sRGB);  title('result');
    linkaxes
end

author: mahmoudnafifi


text:
function stats = chip_histogram_features( varargin )
% ------------
% Description:
% ------------
%  This function is to obtain state of the art histogram based features
% such as:
%   Mean
%   Variance
%   Skewness
%   Kurtosis
%   Energy
%   Entropy
% ---------
% History:
% ---------
% Creation: beta         Date: 09/11/2007
%----------
% Example:
%----------
% Stats = chip_histogram_features( I,'NumLevels',9,'G',[] )
%
% -----------
% Author:
% -----------
%    (C)Xunkai Wei <xunkai.wei@gmail.com>
%    Beijing Aeronautical Technology Research Center
%    Beijing %9203-12,10076
%

% Parameter checking
[I, NL, GL] = ParseInputs(varargin{:});
% Scale I so that it contains integers between 1 and NL.
if GL(2) == GL(1)
    SI = ones(size(I));
else
    slope = (NL-1) / (GL(2) - GL(1));
    intercept = 1 - (slope*(GL(1)));
    SI = round(imlincomb(slope,I,intercept,'double'));
end
% Clip values if user had a value that is outside of the range, e.g., double
% image = [0 .5 2;0 1 1]; 2 is outside of [0,1]. The order of the following
% lines matters in the event that NL = 0.
SI(SI > NL) = NL;
SI(SI < 1) = 1;
%--------------------------------------------------------------------------
% 1. Calculate histogram for all scaled gray level from 1 to NL
%--------------------------------------------------------------------------
% Get image size
s = size(SI);
% Generate gray level vector
Gray_vector = 1:NL;
% intialize parameters
Histogram = zeros(1,NL);
% Using inline function numel, make it easy
for i =1:NL
    Histogram(i) = numel(find(SI==i));
end
%--------------------------------------------------------------------------
% 2. Now calculate its histogram statistics
%--------------------------------------------------------------------------
% Calculate obtains the approximate probability density of occurrence of the intensity
% levels
Prob                = Histogram./(s(1)*s(2));
% 2.1 Mean
Mean                = sum(Prob.*Gray_vector);
% 2.2 Variance
Variance            = sum(Prob.*(Gray_vector-Mean).^2);
% 2.3 Skewness
Skewness            = calculateSkewness(Gray_vector,Prob,Mean,Variance);
% 2.4 Kurtosis
Kurtosis            = calculateKurtosis(Gray_vector,Prob,Mean,Variance);
% 2.5 Energy
Energy              = sum(Prob.*Prob);
% 2.6 Entropy
Entropy             = -sum(Prob.*log(Prob));
%-------------------------------------------------------------------------
% 3. Insert all features and return
%--------------------------------------------------------------------------
stats =[Mean Variance Skewness Kurtosis  Energy  Entropy];
% End of funtion
%--------------------------------------------------------------------------
% Utility functions
%--------------------------------------------------------------------------
function Skewness = calculateSkewness(Gray_vector,Prob,Mean,Variance)
% Calculate Skewness
term1    = Prob.*(Gray_vector-Mean).^3;
term2    = sqrt(Variance);
Skewness = term2^(-3)*sum(term1);

function Kurtosis = calculateKurtosis(Gray_vector,Prob,Mean,Variance)
% Calculate Kurtosis
term1    = Prob.*(Gray_vector-Mean).^4;
term2    = sqrt(Variance);
Kurtosis = term2^(-4)*sum(term1);

function [I, nl, gl] = ParseInputs(varargin)
% parsing parameter checking
% Inputs must be max seven item
iptchecknargin(1,5,nargin,mfilename);
%
% Check I
I = varargin{1};
iptcheckinput(I,{'logical','numeric'},{'2d','real','nonsparse'}, ...
    mfilename,'I',1);
% ------------------------
% Assign Defaults
% -------------------------
%
if islogical(I)
    nl = 2;
else
    nl = 8;
end
gl = getrangefromclass(I);

% Parse Input Arguments
if nargin ~= 1

    paramStrings = {'NumLevels','GrayLimits'};

    for k = 2:2:nargin

        param = lower(varargin{k});
        inputStr = iptcheckstrs(param, paramStrings, mfilename, 'PARAM', k);
        idx = k + 1;  %Advance index to the VALUE portion of the input.
        if idx > nargin
            eid = sprintf('Images:%s:missingParameterValue', mfilename);
            msg = sprintf('Parameter ''%s'' must be followed by a value.', inputStr);
            error(eid,'%s', msg);
        end

        switch (inputStr)
            case 'NumLevels'
                nl = varargin{idx};
                iptcheckinput(nl,{'logical','numeric'},...
                    {'real','integer','nonnegative','nonempty','nonsparse'},...
                    mfilename, 'NL', idx);
                if numel(nl) > 1
                    eid = sprintf('Images:%s:invalidNumLevels',mfilename);
                    msg = 'NL cannot contain more than one element.';
                    error(eid,'%s',msg);
                elseif islogical(I) && nl ~= 2
                    eid = sprintf('Images:%s:invalidNumLevelsForBinary',mfilename);
                    msg = 'NL must be two for a binary image.';
                    error(eid,'%s',msg);
                end
                nl = double(nl);

            case 'GrayLimits'

                gl = varargin{idx};
                iptcheckinput(gl,{'logical','numeric'},{'vector','real'},...
                    mfilename, 'GL', idx);
                if isempty(gl)
                    gl = [min(I(:)) max(I(:))];
                elseif numel(gl) ~= 2
                    eid = sprintf('Images:%s:invalidGrayLimitsSize',mfilename);
                    msg = 'GL must be a two-element vector.';
                    error(eid,'%s',msg);
                end
                gl = double(gl);
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function d = delta(a)

d = a.delta;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [f,g]=ycfun2(x);
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
   x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
   2*x(1)+x(1)^2+2*x(2)+x(3)-40];
g=x(1)^2+x(3)-2;

author: Eurus-Holmes


text:
function my_grad(V,F,u,k)
%MY_GRAD Compute the numerical gradient operator for triangle meshes
%
% my_grad(V,F);
%
% Inputs:
%  V,F  the input mesh for which to compute the gradient operator
% Outputs:
%  G  the sparse gradient matrix
%

G = ...

end

author: odedstein

repo: sgi-introduction-course


text:
%ROBOT	robot object constructor
%
%	ROBOT
%	ROBOT(robot)		create a copy of an existing ROBOT object
%	ROBOT(LINK, ...)	create from a cell array of LINK objects
%	ROBOT(DH, ...)		create from legacy DYN matrix
%	ROBOT(DYN, ...)		create from legacy DYN matrix
%
%	optional trailing arguments are:
%		Name			robot type or name
%		Manufacturer		who built it
%		Comment			general comment
%
%  If the legacy matrix forms are used the default name is the workspace
% variable that held the data.

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function r = SerialLink(L, varargin)
	r.name = 'noname';
	r.manufacturer = '';
	r.comment = '';
	r.links = [];
	r.n = 0;
	r.mdh = 0;
	r.gravity = [0; 0; 9.81];
	r.base = eye(4,4);
	r.tool = eye(4,4);
	r.handle = [];  % graphics handles
	r.q = [];   % current joint angles

	r.lineopt = {'Color', 'black', 'Linewidth', 4};
	r.shadowopt = {'Color', 0.7*[1 1 1], 'Linewidth', 3};
	r.plotopt = {};
	r = class(r, 'SerialLink');

	if nargin == 0
	%zero argument constructor
		return;

	else
		if isa(L, 'SerialLink')
			if length(L) == 1
				% clone the passed robot
				for j=1:L.n
					newlinks(j) = L.links(j);
				end
				r.links = newlinks;
			else
				% compound the robots in the vector
				r = L(1);
				for k=2:length(L)
					r = r * L(k);
				end
			end
		elseif isa(L,'Link')
			r.links = L; %attach the links

		elseif isa(L, 'double')
			% legacy matrix
			dh_dyn = L;
			clear L
			for j=1:numrows(dh_dyn)
				L(j) =        Link();
				L(j).theta =  dh_dyn(j,1);
				L(j).d =      dh_dyn(j,2);
				L(j).a =      dh_dyn(j,3);
				L(j).alpha =  dh_dyn(j,4);

				if numcols(dh_dyn) > 4
					L(j).sigma = dh_dyn(j,5);
				end
			end
			r.links = L;
		else
			error('unknown type passed to SerialLink');
		end
		r.n = length(r.links);
	end
	%process the rest of the arguments in key, value pairs
	opt.name = 'robot';
	opt.comment = [];
	opt.manufacturer = [];
	opt.base = [];
	opt.tool = [];
	opt.offset = [];
	opt.qlim = [];
	opt.plotopt = [];

	[opt,out] = tb_optparse(opt, varargin);
	if ~isempty(out)
		error( sprintf('unknown option <%s>', out{1}));
	end

	% copy the properties to robot object
            p = fieldnames(r);
            for i=1:length(p)
                if isfield(opt, p{i}) && ~isempty(getfield(opt, p{i}))

					r = setfield(r, p{i}, getfield(opt, p{i}));

				end
            end

	% set the robot object mdh status flag
	mdh = [r.links.mdh];
	if all(mdh == 0)
		r.mdh = mdh(1);
	elseif all (mdh == 1)
		r.mdh = mdh(1);
	else
		error('SerialLink has mixed D&H links conventions');
	end
	r = class(r, 'SerialLink');

endfunction

author: petercorke

repo: robotics-toolbox-matlab


text:
function printOcclusionResult(im, bndinfo, lab, fn, fignum)

if ~exist('fignum', 'var') || isempty(fignum)
    fignum = 1;
end

if isfield(bndinfo, 'labels')
    objlab = bndinfo.labels;
elseif isfield(bndinfo, 'result') && isfield(bndinfo.result, 'geomProb')
    objlab = zeros(bndinfo.nseg, 1);
    [maxval, maxlab] = max(bndinfo.result.geomProb, [], 2);
    ind = find(maxlab>=2 & maxlab<=4);
    objlab(ind) = (3:2+numel(ind));
    objlab(maxlab==1) = 1;
    objlab(maxlab==5) = 2;
else
    objlab = [];
end

displayOcclusionResult(im, bndinfo, objlab, lab, {}, fignum);
set(gcf, 'PaperPositionMode', 'auto');
%print(['-f' num2str(fignum)], '-djpeg99', fn);
print(['-f' num2str(fignum)], '-depsc2', [fn(1:end-4) '.eps']);

author: Cloud-CV

repo: object-proposals


text:
function [u,v,px,py] = perform_tv_hilbert_projection(f,kernel,lam,options)

% Aujol & Chambolle projection for solving TV-K regularization
%
%   [u,v,px,py] = perform_tv_hilbert_projection(f,kernel,lam,options);
%
%   Solve the image separation problem :
%       u = argmin_u |f-u|_K^2 + lam*|u|_TV
%
%   f = u+v
%   u is the cartoon part.
%   v is the texture part.
%
%   where |.|_K is an hilber space norm defined by some operator K
%       <f,g>_K = <f,K*g>
%       |f|_K^2 = <f,f>_K
%
%   f is the input image.
%   kernel is a callback function of the kernel operator, should be like
%       f = kernel(f, options);
%
%	Set options.use_gabor==1 to setup the Gabor Kernel automatically
%	and options.use_gabor==-1 to set up the L2 identity kernel.
%
% Original code by Guy Gilboa, modified by Gabriel Peyre
%
%   Copyright (c) 2006 Gabriel Peyre

[ny,nx]=size(f);

options.null = 0;
if isfield(options, 'niter')
    niter = options.niter;
else
    niter = 30;
end
if isfield(options, 'dt')
    dt = options.dt;
else
    dt = 0.01;
end
if isfield(options, 'p0x') && ~isempty(options.p0x)
    p0x = options.p0x;
else
    p0x=zeros(ny,nx);
end
if isfield(options, 'p0y') && ~isempty(options.p0y)
    p0y = options.p0y;
else
    p0y=p0x;
end

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

px=p0x; py=p0y;
lami=1/lam;  % here the algorithm works with inverse of lambda

if use_gabor==1
    sig_k=1; freq=0.25; bias=0;
    n=11; %kernel size
    x=ones(n,1)*(-(n-1)/2:(n-1)/2); y=x';
    gs=exp(-((x/sig_k).^2+(y/sig_k).^2)/2);  % 2D gaussian
    gs=gs/sum(sum(gs)); % normalize
    r=sqrt(x.^2+y.^2); cs=cos(2*pi*r*freq); %radially symmetric
    iKx = cs.*gs; %% "Gabor filter" for inverse K
    iKx=iKx-mean(mean(iKx)); % zero mean filter
    cp=(n+1)/2; % center of kernel
    iKx(cp,cp) = iKx(cp,cp)-bias;
elseif use_gabor==-1
    iKx = 1;
end

%% perform iterations
if use_gabor~=0
    for i=1:niter,  %% do niterations
        progressbar(i,niter);
        %% compute projection
        km1div = filter2(iKx,div(px,py))-f*lam;
        %[Gx,Gy] = grad(km1div);
        Gx=gradx(km1div); Gy=grady(km1div);
        aG = sqrt(Gx.^2+Gy.^2); %abs(G)
        px = (px + dt*Gx)./(1+dt*aG);
        py = (py + dt*Gy)./(1+dt*aG);
    end % for i
    v = filter2(iKx,div(px,py))/lam;
else
    for i=1:niter,  %% do niterations
        progressbar(i,niter);
        % compute projection
        km1div = feval( kernel, div(px,py), options ) - f*lam;
        % compute gradient
        Gx=gradx(km1div); Gy=grady(km1div);
%        [Gx,Gy] = grad(km1div);
        aG = sqrt(Gx.^2+Gy.^2); %abs(G)
        px = (px + dt*Gx)./(1+dt*aG);
        py = (py + dt*Gy)./(1+dt*aG);
    end % for i
    v = lami*feval( kernel, div(px,py), options );
end

u=f-v;

%% additional functions
%%%%%%%%%%%%%%%%%%%%%%%%%%% Gradient (forward difference)
function [fx,fy] = grad(P)
error('Should not be used.');
fx = P(:,[2:end end])-P;
fy = P([2:end end],:)-P;
%%%%%%%%%%%%%%%%%%%%%%%%%%% Divergence (backward difference)
function M=div(px,py)
[m,n]=size(px);
M=zeros(m,n);
Mx=M; My=M;
Mx(2:m-1,1:n)=px(2:m-1,1:n)-px(1:m-2,1:n);
Mx(1,:)=px(1,:);
Mx(m,:)=-px(m-1,:);
My(1:m,2:n-1)=py(1:m,2:n-1)-py(1:m,1:n-2);
My(:,1)=py(:,1);
My(:,n)=-py(:,n-1);
M=Mx+My;
%%%%%%%%%%%%%%%%%%%%%%%%%%% Laplacian (central difference)
function fl = lap(P)
[gx,gy]=grad(P);
fl=div(gx,gy);
%%%%%%%%%%%%%%%%%%%%%%%%%%% Gradient on X (forward difference)
function M=gradx(I)
[m,n]=size(I);
M=zeros(m,n);
M(1:m-1,1:n)=-I(1:m-1,:)+I(2:m,:);
M(m,1:n)=zeros(1,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%% Gradient on Y (forward difference)
function M=grady(I)
[m,n]=size(I);
M=zeros(m,n);
M(1:m,1:n-1)=-I(:,1:n-1)+I(:,2:n);
M(1:m,n)=zeros(m,1);

author: gpeyre

repo: matlab-toolboxes


text:
function y=teager(x,d,m)
%TEAGER calculate teager energy waveform Y=(X,D,M)
%
%  Inputs:  x         speech signal
%           d         dimension to apply filter along [default 1st non-singleton]
%           m         Normally Y has the same length as X and the first
%                     and last output samples are extrapolated. Setting m='x'
%                     supresses this extrapolation and Y will be two
%                     samples shorter than X
%
% Outputs:  Y         output signal: y(n)=abs(x(n))^2 - x(n+1)*conj(x(n-1))
%
% Calculates the Teager energy waveform [1]. The following waveforms give
% a constant output (independent of n) where A, B, C are real constants:
%  (a) x(n) = A*sin(B*n+C)    -->   y(n) = (A*sin(B))^2
%  (b) x(n) = A*n + B         -->   y(n) = A^2
%  (c) x(n) = A*exp(j(B*n+C)) -->   y(n) = A^2*(1-exp(2jB))
%  (d) x(n) = A*exp(B*n+C)    -->   y(n) = 0
%
% Reference:
%  [1]	J. Kaiser. On a simple algorithm to calculate the energy of a signal.
%       In Proc IEEE Intl Conf Acoustics, Speech and Signal Processing,
%       pages 381384, vol.1, Apr. 1990. doi: 10.1109/ICASSP.1990.115702.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: teager.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e=size(x);
p=prod(e);
if nargin<2             % if no dimension given, find the first non-singleton
    d=find(e>1,1);
    if ~numel(d)
        d=1;
    end
end
k=e(d);                 % size of active dimension
q=p/k;                  % size of remainder
if d==1
    z=reshape(x,k,q);
else
    z=shiftdim(x,d-1);
    r=size(z);
    z=reshape(z,k,q);
end
if nargin>2 && any(m=='x')
    y=z(2:k-1,:).*conj(z(2:k-1,:))-z(3:k,:).*conj(z(1:k-2,:));
    k=k-2;              % we have lost two elements
elseif k>=4
    y=zeros(k,q);
    y(2:k-1,:)=z(2:k-1,:).*conj(z(2:k-1,:))-z(3:k,:).*conj(z(1:k-2,:));
    y(1,:)=2*y(2,:)-y(3,:);             % linearly interpolate the end points
    y(k,:)=2*y(k-1,:)-y(k-2,:);
elseif k==3
    y=repmat(x(2,:).*conj(x(2,:))-x(3,:).*conj(x(1,:)),3,1);
else
    y=zeros(k,q);
end
if d==1
    e(d)=k;
    y=reshape(y,e);
else
    r(1)=k;
    y=shiftdim(reshape(y,r),length(e)+1-d);
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [ cuboid ] = getSegHypothesisID( labelMap, minSize, vanishing, ID )
%GETBINARYMAP Summary of this function goes here
%   Detailed explanation goes here
segment = labelMap;% + 1;
numSeg = round(max(segment(:)));
[H,W] = size(segment);

invalidMask = false(H,W);
invalidMask(:,1:10) = true;   invalidMask(:,W-9:W) = true;
invalidMask(1:10,:) = true;   invalidMask(H-9:H,:) = true;

fprintf('*1');

if nargin<=3
    n = histc(segment(:),1:numSeg);
    segValid = n>minSize;
else
    segValid = false(numSeg,1);
    segValid(ID) = true;
end
counter = 0;
% maxHypot= struct([]);
% hypot = zeros(6,2,1000);
% score = zeros(1,1000);

% binaryMaps = zeros(H,W,100);
SE = strel('disk', 5, 4);
SS = strel('disk', 2, 4);

fprintf('*2');

for i = find(segValid)'
    seg = segment == i;

    binaryMap = imdilate( seg, SE);
    binaryMap = imerode( binaryMap, SE);
    binaryMap = imerode( binaryMap, SS);
    binaryMap = imdilate( binaryMap, SS);
    binaryMap = imfill(binaryMap, 'hole');

    CC = bwconncomp(binaryMap);
    valid = true(CC.NumObjects,1);
    for j = 1:CC.NumObjects
        if length(CC.PixelIdxList{j})<minSize
            valid(j) = false;
        end
        if any(invalidMask(CC.PixelIdxList{j}))
            valid(j) = false;
        end
    end
    CC.NumObjects = sum(valid);
    CC.PixelIdxList = CC.PixelIdxList(valid);

%     R = regionprops( CC, 'BoundingBox', 'Area', 'ConvexArea', 'Centroid');
    for j = 1:CC.NumObjects
        subMask = false(H,W);
        subMask(CC.PixelIdxList{j}) = true;
        counter = counter + 1;
        maxHypot(counter) = region2cuboid( subMask, vanishing );

%         figure; imshow(subMask); hold on;
%         color = [1 0 0; 0.5 0 0; 0 1 0; 0 0.5 0; 0 0 1; 0 0 0.5];
%         points = maxHypot(counter).points;
%         for k = 1:6
%             line(points([1 2 3 4 5 6 1],1), points([1 2 3 4 5 6 1],2), 'Color', color(k,:), 'LineWidth', 2);
%         end
%         for k = 1:6
%             plot(points(k,1), points(k,2), '--rs','LineWidth',2,...
%                         'MarkerEdgeColor','k',...
%                         'MarkerFaceColor',color(k,:),...
%                         'MarkerSize',5);
%         end

    end

end

fprintf('*3');

if exist('maxHypot','var')
    score = [maxHypot.score];
    cuboid = maxHypot(score>0.75);
else
    cuboid = [];
end

fprintf('*4');

%%
% numCuboid = length(score);
% % rectangle = repmat(struct('xyzBox',zeros(1000,12),'xyBox',zeros(1000,8),'count',0),3,1);
% cuboid = struct('views', zeros(numCuboid,3), 'xyBox', zeros(numCuboid,24), 'score', zeros(numCuboid,1));
% for i = 1:numCuboid
%     points = hypot(:,:,i);
%     [x2,y2] = poly2cw(points(:,1), points(:,2));
%     points = [x2 y2];
%
%     lines = points([2 3 4 5 6 1],:) - points([1 2 3 4 5 6],:);
%     direct = lines./repmat(sqrt(sum(lines.^2,2)), 1, 2);
%     midpoints = (points([2 3 4 5 6 1],:) + points([1 2 3 4 5 6],:))/2;
%     A = zeros(6,3);
%     for j = 1:3
%         vanline = [midpoints(:,1)-vanishing(j,1) midpoints(:,2)-vanishing(j,2)];
%         vandire = vanline./repmat(sqrt(sum(vanline.^2,2)), 1, 2);
%         A(:,j) = dot(direct, vandire, 2);
%     end
%     [~,direct_assign] = max( abs(A), [], 2);
%     % define type
%     type_indicate = direct_assign(1:3)==direct_assign(4:6);
%     if sum(type_indicate)==1 % two faces
%         I = find(type_indicate);
%         pointID1 = rem((I-1)+2,6)+1;
%         pointID2 = rem((pointID1-1)+3,6)+1;
%         testdir = points(pointID1,:)-points(pointID2,:);
%         testdir = testdir ./ norm(testdir);
%         refedir = vanishing(direct_assign(I),:) - (points(pointID1,:)+points(pointID2,:))/2;
%         refedir = refedir ./ norm(refedir);
%         if abs(dot(testdir, refedir))>0.95 % valid, connect line follow 3rd direction
%             pointID = rem([I I+1 I+2 I+5]-1,6)+1;
%             directIDs = direct_assign(pointID([4 1 2]));
%             directID = setdiff(1:3,directIDs);
% %             count1 = rectangle(directID).count + 1;
% %             rectangle(directID).xyBox(count1,:) = reshape(points(pointID,:)',1,[]);
% %             rectangle(directID).count = rectangle(directID).count + 1;
%             cuboid.xyBox(i,1:8) = reshape(points(pointID,:)',1,[]);
%             cuboid.views(i,1) = directID;
%
%             pointID = rem([I+2 I+3 I+4 I+5]-1,6)+1;
%             directIDs = direct_assign(pointID([1 2 3]));
%             directID = setdiff(1:3,directIDs);
% %             count2 = rectangle(directID).count + 1;
% %             rectangle(directID).xyBox(count2,:) = reshape(points(pointID,:)',1,[]);
% %             rectangle(directID).count = rectangle(directID).count + 1;
%             cuboid.xyBox(i,9:16) = reshape(points(pointID,:)',1,[]);
%             cuboid.views(i,2) = directID;
%
% %             cuboid.combo(i,1:2) = [count1 count2];
%             cuboid.score(i) = score(i);
%         end
%     elseif sum(type_indicate)==3
%         % three faces
% %         [~,lowID] = max(points(:,2));
% %         ptrID = rem( lowID+[-2 0 +2 -1 +1 +3]-1, 6) + 1;
% %         ptr = zeros(7,2);
% %         ptr(2:7,:) = points(ptrID,:);
% %         [point2D, point3D, K] = fitCuboid(ptr(:,[2 1]), 1);
% %         points(ptrID,:) = point2D([2 1],2:7)';
%
% %         pointIDs = zeros(1,3);
% %         directIDs = zeros(1,3);
%          I = find(direct_assign==3);
%          verticalID = rem([I(1)+2 I(1)+5]-1,6) + 1;
%         if norm(points(verticalID(1),:)-vanishing(3,:))<norm(points(verticalID(2),:)-vanishing(3,:))
%             pointIDs = rem([verticalID(1) verticalID(1)+2 verticalID(1)+4]-1,6)+1;
%         else
%             pointIDs = rem([verticalID(2) verticalID(2)+2 verticalID(2)+4]-1,6)+1;
%         end
%
% %         if points(verticalID(1),2)<points(verticalID(2),2)
% %             pointIDs = rem([verticalID(2) verticalID(2)+2 verticalID(2)+4]-1,6)+1;
% %         else
% %             pointIDs = rem([verticalID(1) verticalID(1)+2 verticalID(1)+4]-1,6)+1;
% %         end
%         directIDs = direct_assign(rem(pointIDs-1+1,6)+1);
% %         for j = 1:3
% %             directID = direct_assign( j+1 );
% %             directIDs(j) = directID;
% %             if norm(points(j,:)-vanishing(directID,:))<norm(points(j+3,:)-vanishing(directID,:))
% %                 pointIDs(j) = j;
% %             else
% %                 pointIDs(j) = j+3;
% %             end
% %         end
%         lines = [points(pointIDs(1),:) vanishing(directIDs(1),:); ...
%                  points(pointIDs(2),:) vanishing(directIDs(2),:); ...
%                  points(pointIDs(3),:) vanishing(directIDs(3),:)];
%         intersect = intersectPoint(lines([1 2 3],:), lines([2 3 1],:));
% %         if sum(sum((intersect([1 2 3],:)-intersect([2 3 1],:)).^2,2))<10
%
%             addPoint = sum(intersect,1)./3;
%
% %             count1 = rectangle(directIDs(1)).count + 1;
% %             rectangle(directIDs(1)).count = rectangle(directIDs(1)).count + 1;
% %             rectangle(directIDs(1)).xyBox(count1,:) = reshape([points([1 2 3],:);addPoint]',1,[]);
%             cuboid.xyBox(i,1:8) = reshape([points(rem(pointIDs(1)+[0 1 2]-1,6)+1,:);addPoint]',1,[]);
%             cuboid.views(i,1) = setdiff(1:3, direct_assign([1 2]));
% %             count2 = rectangle(directIDs(2)).count + 1;
% %             rectangle(directIDs(2)).count = rectangle(directIDs(2)).count + 1;
% %             rectangle(directIDs(2)).xyBox(count2,:) = reshape([points([3 4 5],:);addPoint]',1,[]);
%             cuboid.xyBox(i,9:16) = reshape([points(rem(pointIDs(2)+[0 1 2]-1,6)+1,:);addPoint]',1,[]);
%             cuboid.views(i,2) = setdiff(1:3, direct_assign([3 4]));
% %             count3 = rectangle(directIDs(3)).count + 1;
% %             rectangle(directIDs(3)).count = rectangle(directIDs(3)).count + 1;
% %             rectangle(directIDs(3)).xyBox(count3,:) = reshape([points([5 6 1],:);addPoint]',1,[]);
%             cuboid.xyBox(i,17:24) = reshape([points(rem(pointIDs(3)+[0 1 2]-1,6)+1,:);addPoint]',1,[]);
%             cuboid.views(i,3) = setdiff(1:3, direct_assign([5 6]));
%
% %             cuboid.views(i,:) = directIDs;
% %             cuboid.combo(i,:) = [count1 count2 count3];
%             cuboid.score(i) = score(i);
% %         end
%     end
% end
% valid = cuboid.score>0.8;
% cuboid.views = cuboid.views(valid,:);
% cuboid.score = cuboid.score(valid);
% cuboid.xyBox = cuboid.xyBox(valid,:);
% % cuboid.combo = cuboid.combo(valid,:);
% % cuboid.xyzBox = zeros(length(cuboid.score),12);
% % for i = 1:length(cuboid.score)
% %
% % end
% %

end

author: zouchuhang

repo: LayoutNet


text:
function [f,g,H] = autoHess(x,useComplex,funObj,varargin)
% Numerically compute Hessian of objective function from gradient values

p = length(x);

if useComplex % Use Complex Differentials
    mu = 1e-150;

    diff = zeros(p);
    for j = 1:p
        e_j = zeros(p,1);
        e_j(j) = 1;
        [f(j) diff(:,j)] = funObj(x + mu*i*e_j,varargin{:});
    end
    f = mean(real(f));
    g = mean(real(diff),2);
    H = imag(diff)/mu;
else % Use finite differencing
    mu = 2*sqrt(1e-12)*(1+norm(x))/norm(p);

    [f,g] = funObj(x,varargin{:});
    diff = zeros(p);
    for j = 1:p
        e_j = zeros(p,1);
        e_j(j) = 1;
        [f diff(:,j)] = funObj(x + mu*e_j,varargin{:});
    end
    H = (diff-repmat(g,[1 p]))/mu;
end

% Make sure H is symmetric
H = (H+H')/2;

if 0 % DEBUG CODE
    [fReal gReal HReal] = funObj(x,varargin{:});
    [fReal f]
    [gReal g]
    [HReal H]
    pause;
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [dlnPart, m] = gradLnDiffErfs(x1, x2, fact1, fact2),

% GRADLNDIFFERFS Compute the gradient of the log difference of two erfs.
%
%
% COPYRIGHT : Antti Honkela, 2007

% NDLUTIL

m = min(x1.^2, x2.^2);
dlnPart = 2/sqrt(pi) * (exp(-x1.^2 + m) .* fact1 ...
			- exp(-x2.^2 + m) .* fact2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [legal] = isLegal(v)
legal = sum(any(imag(v(:))))==0 & sum(isnan(v(:)))==0 & sum(isinf(v(:)))==0;

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function a = householder ( n, x )

%*****************************************************************************80
%
%% HOUSEHOLDER constructs a HOUSEHOLDER matrix.
%
%  Discussion:
%
%    A Householder matrix is also called an elementary reflector.
%
%  Formula:
%
%     A = I - ( 2 * X * X' ) / ( X' * X )
%
%  Example:
%
%    N = 5, X = ( 1, 1, 1, 0, -1 )
%
%   1/2 -1/2 -1/2  0  1/2
%  -1/2  1/2 -1/2  0  1/2
%  -1/2 -1/2  1/2  0  1/2
%    0    0    0   1   0
%   1/2  1/2  1/2  0  1/2
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is orthogonal: A' * A = A * A' = I.
%
%    inverse ( A ) = A.
%
%    det ( A ) = -1.
%
%    A is unimodular.
%
%    If Y and Z are nonzero vectors of equal length, and
%      X = ( Y - Z ) / NORM(Y-Z),
%    then
%      A * Y = Z.
%
%    A represents a reflection through the plane which
%    is perpendicular to the vector X.  In particular, A*X = -X.
%
%    LAMBDA(1) = -1;
%    LAMBDA(2:N) = +1.
%
%    If X is the vector used to define H, then X is the eigenvector
%    associated with the eigenvalue -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub, Charles Van Loan,
%    Matrix Computations, second edition,
%    Johns Hopkins University Press, Baltimore, Maryland, 1989.
%
%    Pete Stewart,
%    Introduction to Matrix Computations,
%    Academic Press, 1973,
%
%    James Wilkinson,
%    The Algebraic Eigenvalue Problem,
%    Oxford University Press, 1965.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(N), the vector that defines the
%    Householder matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    a(i,i) = 1.0;
  end

  xdot = x(1:n) * x(1:n)';

  if ( 0.0 < xdot )

    for i = 1 : n
      for j = 1 : n
        a(i,j) = a(i,j) - 2.0 * x(i) * x(j) / xdot;
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%*******************************************************
%%  Doptdesign: D-optimal experiment design.
%%
%%  max  log(det(sum_{i=1}^p lambda_i v_i v_i^T)) + const
%%
%%  s.t. lambda_i >= 0, sum_{i=1}^p lambda_i = 1.
%%
%%  Note: const = n.
%%  V:       nxp matrix with n <= p.
%%  lambda:  lambda_i is the fraction of the experiments
%%           allocated to test vector v_i.
%%  S: = V*diag(lambda)*V'.
%%
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*******************************************************

   function [blk,At,C,b,OPTIONS,lambda,bblk,AAt] = Doptdesign(V,solve);

   if (nargin == 1); solve = 0; end

   [n,p] = size(V);
   if (n > p);
     error(' size(V,1) > size(V,2)');
   end
%%
%% form blk, At, C, b
%%
   b = zeros(p,1);

   blk{1,1} = 's'; blk{1,2} = n;
   F = cell(1,p);
   for k = 1:p
      F{1,k} = -V(:,k)*V(:,k)';
   end
   At(1) = svec(blk(1,:),F,1);
   C{1,1} = sparse(n,n);

   blk{2,1} = 'l'; blk{2,2} = p;
   At{2,1} = -speye(p,p);
   C{2,1} = zeros(p,1);

   blk{3,1} = 'u'; blk{3,2} = 1;
   At{3,1} = ones(1,p);
   C{3,1} = 1;

   OPTIONS.parbarrier{1,1} = 1;
   OPTIONS.parbarrier{2,1} = 0;
   OPTIONS.parbarrier{3,1} = 0;
%%
%% form bblk, AAt to take into account of
%% low-rank constraint matrices of the form: -vk*vk'.
%%
    bblk = blk; AAt = At;
    bblk{1,1} = 's'; bblk{1,2} = n; bblk{1,3} = ones(1,p);
    AAt{1,1} = []; AAt{1,2} = V; AAt{1,3} = -ones(p,1);
%%
   if (solve)
      [obj,X,y,Z] = sqlp(blk,At,C,b,OPTIONS);
      lambda = -y;
   else
      lambda = [];
   end
%%*******************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% makenoise.m
% written by: Duncan Po
% Date: August 24, 2002
% Generate additive Gaussian noise and add the noise to an image
% This function assumes that the image is stored in uint8 format with
% values from 0 to 255
% Usage: noisepic = makenoise(imname, imformat, nvar)
% Inputs:   imname      - name of the image file
%           imformat    - format of the image file
%           nvar        - variance of the noise normalized to the image range
% Output:   noisepic    - the resulting noisy image

function noisepic = makenoise(imname, imformat, nvar)

pic = imread(imname,imformat);
figure;
imshow(pic);
pic = double(pic);

% the following step is needed for the command imnoise (normalize)
pic2 = pic/255;

% adds zero mean Gaussian noise of power (0.1*255)^2 = 650.25
noisepic = imnoise(pic2, 'Gaussian', 0, nvar);

% converts the image back to a scale of 8 bits (0 to 255)
noisepic = noisepic*255;

% calculate initial MSE of the noisy image
MSE = sqrt(sum(sum((noisepic - pic).*(noisepic-pic)/(size(pic,1)*size(pic,2))...
    /(size(pic,1)*size(pic,2)))))
noisepic = uint8(noisepic);

figure;
imshow(noisepic);


text:
% getkeyval() - get variable value from a 'key', 'val' sequence string.
%
% Usage:
%   >> val = getkeyval( keyvalstr, varname, mode, defaultval);
%
% Inputs:
%   keyvalstr  - string containing 'key', 'val' arguments
%   varname    - string for the name of the variable or index
%                of the value to retrieve (assuming arguments are
%                separated by comas).
%   mode       - if the value extracted is an integer array, the
%                'mode' variable can contain a subset of indexes to return.
%                If mode is 'present', then either 0 or 1 is returned
%                depending on wether the variable is present.
%   defaultval - default value if the varible is not found
%
% Outputs:
%   val  - a value for the variable
%
% Note: this function is helpful for finding arguments in string commands
%       stored in command history.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 29 July 2002
%
% See also: gethelpvar()

% Copyright (C) 2002 arno@salk.edu, Arnaud Delorme, CNL / Salk Institute
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function txt = getkeyval(lastcom, var, mode, default)
    % mode can be present for 0 and 1 if the variable is present
	if nargin < 1
        help getkeyval;
        return;
    end;
    if nargin < 4
		default = '';
	end;
	if isempty(lastcom)
		txt = default; return;
	end;
	if nargin < 3
		mode = '';
	end;
	if isstr(mode) & strcmp(mode, 'present')
		if ~isempty(findstr(var, lastcom))
			txt = 1; return;
		else
			txt = 0; return;
		end;
	end;
	if isnumeric(var)
		comas = findstr(lastcom, ',');
		if length(comas) >= var
			txt = lastcom(comas(var-1)+1:comas(var)-1);
            tmpval = eval( txt );
            if isempty(tmpval), txt = '';
            else txt = vararg2str( tmpval );
            end;
            return;

% 			txt = deblank(txt(end:-1:1));
% 			txt = deblank(txt(end:-1:1));
%
% 			if ~isempty(txt) & txt(end) == '}', txt = txt(2:end-1); end;
% 			if ~isempty(txt)
% 				txt = deblank(txt(end:-1:1));
% 				txt = deblank(txt(end:-1:1));
% 			end;
% 			if ~isempty(txt) & txt(end) == ']', txt = txt(2:end-1); end;
% 			if ~isempty(txt)
% 				txt = deblank(txt(end:-1:1));
% 				txt = deblank(txt(end:-1:1));
% 			end;
% 			if ~isempty(txt) & txt(end) == '''', txt = txt(2:end-1); end;
		else
			txt = default;
		end;
		%fprintf('%s:%s\n', var, txt);
		return;
	else
		comas  = findstr(lastcom, ','); % finding all comas
		comas  = [ comas  findstr(lastcom, ');') ]; % and end of lines
		varloc = findstr(lastcom, var);
		if ~isempty(varloc)
			% finding comas surrounding 'val' var in 'key', 'val' sequence
			comas = comas(find(comas >varloc(end)));
			txt = lastcom(comas(1)+1:comas(2)-1);
			txt = deblank(txt(end:-1:1));
			txt = deblank(txt(end:-1:1));
			if strcmp(mode, 'full')
				parent = findstr(lastcom, '}');
				if ~isempty(parent)
					comas = comas(find(comas >parent(1)));
					txt = lastcom(comas(1)+1:comas(2)-1);
				end;
				txt = [ '''' var ''', ' txt ];
			elseif isnumeric(mode)
				txt = str2num(txt);
				if ~isempty(mode)
					if length(txt) >= max(mode)
						if all(isnan(txt(mode))), txt = '';
						else txt = num2str(txt(mode));
						end;
					elseif length(txt) >= mode(1)
						if all(isnan(txt(mode(1)))), txt = '';
						else txt = num2str(txt(mode(1)));
						end;
					else
						txt = default;
					end;
				else
                    txt = num2str(txt);
				end;
            elseif txt(1) == ''''
                txt = txt(2:end-1); % remove quotes for text
			end;
		else
			txt = default;
		end;
	end;
	%fprintf('%s:%s\n', var, txt);

author: PatternRecognition

repo: OpenBMI


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with %regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

% Cost function
h = sigmoid(X * theta);
theta1 = [0; theta(2:size(theta)) ];

cost = (1 / m) * sum( (-y .* log(h)) - ((1 - y) .* log(1 - h)) );
reg_cost = (lambda / (2 * m)) .* (theta1' * theta);

J = cost + reg_cost;

grad = (1 / m) * X' * (h - y);
reg_grad = (lambda / m) .* theta1;

grad = grad + reg_grad;

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

grad = grad(:);

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function options = opticheckset(Names,Defaults,checkfun,varargin)
%OPTICHECKSET  Default Algorithm for dealing with solver option setting

%   Copyright (C) 2013 Jonathan Currie (I2C2)

%Transpose as required
if(size(Names,2) > 1), Names = Names'; end
if(size(Defaults,2) > 1), Defaults = Defaults'; end

%Collect Sizes and lowercase matches
m = size(Names,1);
numberargs = numel(varargin);
%Create structure with all names and default values
st = [Names,Defaults]'; options = struct(st{:});

% Check we have char or structure input. If structure, insert arguments
i = 1;
while i <= numberargs
    arg = varargin{i};
    if ischar(arg)
        break;
    end
    if ~isempty(arg)
        if ~isa(arg,'struct')
            error('An argument was supplied that wasn''t a string or struct!');
        end
        for j = 1:m
            if any(strcmp(fieldnames(arg),Names{j,:}))
                val = arg.(Names{j,:});
            else
                val = [];
            end
            if ~isempty(val)
                checkfun(Names{j,:},val);
                options.(Names{j,:}) = val;
            end
        end
    end
    i = i + 1;
end

%Check we have even number of args left
if rem(numberargs-i+1,2) ~= 0
    error('You do not have a value specified for each field!');
end

%Go through each argument pair and assign to correct field
expectval = 0; %first arg is a name
while i <= numberargs
    arg = varargin{i};
    switch(expectval)
        case 0 %field
            if ~ischar(arg)
                error('Expected fieldname (argument %d) to be a string!',i);
            end
            j = find(strcmp(arg,Names) == 1);
            if isempty(j)  % if no matches
                error('Unrecognised parameter ''%s''',arg);
            elseif(length(j) > 1)
                error('Ambiguous parameter ''%s''',arg);
            end
            expectval = 1; %next arg is a value
        case 1
            if(~isempty(arg))
                if(ischar(arg)), arg = lower(arg); end
                checkfun(Names{j,:},arg);
                options.(Names{j,:}) = arg;
            end
            expectval = 0;
    end
    i = i + 1;
end

if expectval %fallen off end somehow
    error('Missing value for final argument ''%s''',arg);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function y=rdct(x,n,a,b)
%RDCT     Discrete cosine transform of real data Y=(X,N,A,B)
% Data is truncated/padded to length N.
%
% This routine is equivalent to multiplying by the matrix
%
%   rdct(eye(n)) = diag([sqrt(2)*B/A repmat(2/A,1,n-1)]) * cos((0:n-1)'*(0.5:n)*pi/n)
%
% Default values of the scaling factors are A=sqrt(2N) and B=1 which
% results in an orthogonal matrix. Other common values are A=1 or N and/or B=1 or sqrt(2).
% If b~=1 then the columns are no longer orthogonal.
%
% see IRDCT for the inverse transform

% BUG: in line 51 we should do chopping after transform and not before

fl=size(x,1)==1;
if fl x=x(:); end
[m,k]=size(x);
if nargin<2 n=m;
end
if nargin<4 b=1;
    if nargin<3 a=sqrt(2*n);
    end
    end
if n>m x=[x; zeros(n-m,k)];
elseif n<m x(n+1:m,:)=[];
end

x=[x(1:2:n,:); x(2*fix(n/2):-2:2,:)];
z=[sqrt(2) 2*exp((-0.5i*pi/n)*(1:n-1))].';
y=real(fft(x).*z(:,ones(1,k)))/a;
y(1,:)=y(1,:)*b;
if fl y=y.'; end

author: bastamon


text:
function [Q1,Q2] = mpc(T1_meas,T1_sp,T2_meas,T2_sp)

s = 'http://byu.apmonitor.com';
c = 'mpc';

% input measurement
apm_meas(s,c,'TC1',T1_meas);
apm_meas(s,c,'TC2',T2_meas);

% input setpoint with deadband +/- DT
DT = 0.1;
apm_option(s,c,'TC1.sphi',T1_sp+DT);
apm_option(s,c,'TC1.splo',T1_sp-DT);
apm_option(s,c,'TC2.sphi',T2_sp+DT);
apm_option(s,c,'TC2.splo',T2_sp-DT);

% solve MPC
output = apm(s,c,'solve');

% test for successful solution
if (apm_tag(s,c,'apm.appstatus')==1)
    % retrieve the first Q value
    Q1 = apm_tag(s,c,'Q1.Newval');
    Q2 = apm_tag(s,c,'Q2.Newval');
else
    % display output for debugging
    disp(output)
    % not successful, set voltage to zero
    Q1 = 0;
    Q2 = 0;
end

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function planC = deleteStructure(planC, structNumV)
% function planC = deleteStructure(planC, structNumV)
%
% Function to delete structure in a batch mode.
%
% APA, 09/25/2012

indexS = planC{end};
planC = delUniformStr(structNumV, planC); %Update the uniform data.
planC{indexS.structures}(structNumV) = [];

%len = length(planC{indexS.structures});
%planC = delUniformStr(structNum, planC); %Update the uniform data.
%planC{indexS.structures}(structNum:len-1) = planC{indexS.structures}(structNum+1:len);
%planC{indexS.structures} = planC{indexS.structures}(1:len-1);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
  function ob = Gblock(base, nblock)
%|function ob = Gblock(base, nblock)
%|
%| Construct Gblock object, which is a 'meta' system object,
%| designed to work with ordered-subsets (aka block iterative) algorithms.
%| See Gblock_test.m for example usage.
%|
%| A Gblock system is created from another 'block-izable' base system object
%| 'base' by calling: Gb = Gblock(base, nblock, 1);
%|
%| Any projector object can be used if nblock=1, otherwise the base object
%| must have a "mtimes_block" method for block multiplications.
%|
%| The overloaded capabilities of a Gblock object include:
%|	Ab * x			full forward projection
%|	Ab' * y			full back projection
%|	Ab{i_block} * x		project one block (e.g., a set of views)
%|	Ab{i_block}' * x	back project one block
%|				i_block = 1,2,...,nubset
%|
%| Other properties are inherited from the base system object.
%|
%| Copyright 2002-2-18, Jeff Fessler, University of Michigan

if nargin < 2, nblock = 1; end

% create default object, as required by Mathworks
ob.base = [];
ob.nblock = nblock;
ob.i_block = 0;

if nargin == 0	% required by Mathworks
	ob = class(ob, 'Gblock');
return
end

% handle matrix as a special case
if isnumeric(base)
	error 'Gblock no longer supports numeric; call Gmatrix or Gsparse first'
%	r = input('do you wish to continue? [n|y]', 's');
%	if ~streq(r, 'y'), error '', end
%	base = Gsparse(base);
%	if nblock ~= 1, error 'call Gsparse first', end
end

% ensure that a "mtimes_block" method is available (if needed)
if nblock > 1
	try
		mtimes_block(base, 'exists');
	catch
		error 'this base object does not have a mtimes_block method'
	end
end

ob.base = base;
ob = class(ob, 'Gblock');

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren, Li Xu, Qiong Yan, Wenxiu Sun, "Shepard Convolutional Neural Networks",
% Advances in Neural Information Processing Systems (NIPS 2015)
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
addpath applications/Shepard_CNN/Shepard_super_res/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath optimization/
addpath pipeline/
addpath data/

clearvars -global config;
clearvars -global mem;
clear gen_mask_patch_cat_idx_for_super_res;
clear;
global config mem;
shepard_sr_x2_configure();
init(0);

load('data/Shepard_CNN/Shepard_super_res/x2/val_1ch/val_1');

images_t = reshape(images_t, size(images_t,1), size(images_t,2), 1, size(images_t,3));
labels_t = reshape(labels_t, size(labels_t,1), size(labels_t,2), 1, size(labels_t,3));

perm = randperm(size(images_t, 4));
images_t = images_t(:,:,:,perm);
labels_t = labels_t(:,:,:,perm);
test_samples = config.NEW_MEM(images_t(:,:,:,1:1000));
test_labels = config.NEW_MEM(labels_t(:,:,:,1:1000));

mask = config.NEW_MEM([1 0; 0 0]);
mask = repmat(mask, config.input_size(1)/2, config.input_size(2)/2, config.chs);
mask = repmat(mask, 1,1,1,config.batch_size);

count = 0;
cost_avg = 0;
epoc = 0;
points_seen = 0;
display_points = 5000;
save_points = 50000;
max_grad = 1;
fprintf('%s\n', datestr(now, 'dd-mm-yyyy HH:MM:SS FFF'));
for pass = 1:10
    for p = 1:50
        load(strcat('data/Shepard_CNN/Shepard_super_res/x2/train_1ch/patches_', num2str(p), '.mat'));

        images = reshape(images, size(images,1), size(images,2), 1, size(images,3));
        labels = reshape(labels, size(labels,1), size(labels,2), 1, size(labels,3));

        perm = randperm(20000);
        images = images(:,:,:,perm);
        labels = labels(:,:,:,perm);
        train_imgs = config.NEW_MEM(images);
        train_labels = config.NEW_MEM(labels);

        for i = 1:size(train_labels, 4) / config.batch_size
            points_seen = points_seen + config.batch_size;
            in = train_imgs(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            %in = train_labels(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            out = train_labels(:,:,:,(i-1)*config.batch_size+1:i*config.batch_size);
            out = out((size(in, 1) - config.output_size(1)) / 2 + 1:(size(in, 1) - config.output_size(1)) / 2 + config.output_size(1), ...
                      (size(in, 2) - config.output_size(2)) / 2 + 1:(size(in, 2) - config.output_size(2)) / 2 + config.output_size(2), :, :);

            % make the mask list
            mask_li = {};
            mask_li{1} = mask;
            % operate the training pipeline
            op_train_pipe_with_mask(in.*mask, mask_li, out);
            % update the weights
            config.UPDATE_WEIGHTS();

            if(cost_avg == 0)
                cost_avg = config.cost;
            else
                cost_avg = (cost_avg + config.cost) / 2;
            end

            % display point
            if(mod(points_seen, display_points) == 0)
                count = count + 1;
                fprintf('%d ', count);
            end
            % save point
            if(mod(points_seen, save_points) == 0)
                fprintf('\n%s', datestr(now, 'dd-mm-yyyy HH:MM:SS FFF'));
                epoc = epoc + 1;
                test_cost = 0;
                for t = 1:size(test_samples, 4) / config.batch_size
                    t_label = test_labels(:,:,:,(t-1)*config.batch_size+1:t*config.batch_size);
                    t_label = config.NEW_MEM(t_label((size(in, 1) - config.output_size(1)) / 2 + 1:(size(in, 1) - config.output_size(1)) / 2 + config.output_size(1), ...
                                            (size(in, 2) - config.output_size(2)) / 2 + 1:(size(in, 2) - config.output_size(2)) / 2 + config.output_size(2), :));

                    op_test_pipe_with_mask(test_samples(:,:,:,(t-1)*config.batch_size+1:t*config.batch_size).*mask, mask_li, t_label);

                    test_out = gather(mem.output);
                    test_cost = test_cost + config.cost;
                end
                test_cost = test_cost / (size(test_samples, 4) / config.batch_size);
                fprintf('\nepoc %d, training avg cost: %f, test avg cost: %f\n', epoc, cost_avg, test_cost);

                save_weights(strcat('applications/Shepard_CNN/Shepard_super_res/results/shepard_layer_x2/epoc', num2str(epoc), '.mat'));

                cost_avg = 0;
            end
        end
    end
end

author: jimmy-ren


text:
function testAxis()
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

hFig = figure;
hAxis = gca;
table('init', hAxis);
table('data', hAxis, dir);
set(hAxis, 'units', 'normalized');
set(hFig, 'resizefcn', 'myResize');


text:
function [AF,BF]=wpfbtbounds(wt,varargin)
%WPFBTBOUNDS Frame bounds of WPFBT
%   Usage: fcond=wpfbtbounds(wt,L);
%          [A,B]=wpfbtbounds(wt,L);
%          [...]=wpfbtbounds(wt);
%
%   `wpfbtbounds(wt,L)` calculates the ratio $B/A$ of the frame bounds
%   of the wavelet packet filterbank specified by *wt* for a system of length
%   *L*. The ratio is a measure of the stability of the system.
%
%   `wpfbtbounds(wt)` does the same, except *L* is chosen to be the next
%   compatible length bigger than the longest filter from the identical
%   filterbank.
%
%   `[A,B]=wpfbtbounds(...)` returns the lower and upper frame bounds
%   explicitly.
%
%   See |wfbt| for explanation of parameter *wt*.
%
%   Additionally, the function accepts the following flags:
%
%   `'intsqrt'`(default),`'intnoscale'`, `'intscale'`
%       The filters in the filterbank tree are scaled to reflect the
%       behavior of |wpfbt| and |iwpfbt| with the same flags.
%
%   `'scaling_notset'`(default),`'noscale'`,`'scale'`,`'sqrt'`
%     Support for scaling flags as described in |uwpfbt|. By default,
%     the bounds are caltulated for |wpfbt|, passing any of the non-default
%     flags results in bounds for |uwpfbt|.
%
%   See also: wpfbt, filterbankbounds

% AUTHOR: Zdenek Prusa

complainif_notenoughargs(nargin,1,'WPFBTBOUNDS');

definput.keyvals.L = [];
definput.flags.interscaling = {'intsqrt', 'intscale', 'intnoscale'};
definput.import = {'uwfbtcommon'};
definput.importdefaults = {'scaling_notset'};
[flags,~,L]=ltfatarghelper({'L'},definput,varargin);

wt = wfbtinit({'strict',wt},'nat');

if ~isempty(L) && flags.do_scaling_notset
   if L~=wfbtlength(L,wt)
       error(['%s: Specified length L is incompatible with the length of ' ...
              'the time shifts.'],upper(mfilename));
   end;
end

for ii=1:numel(wt.nodes)
   a = wt.nodes{ii}.a;
   assert(all(a==a(1)),sprintf(['%s: One of the basic wavelet ',...
                                'filterbanks is not uniform.'],...
                                upper(mfilename)));
end

% Do the equivalent filterbank using multirate identity property
[gmultid,amultid] = wpfbt2filterbank(wt,flags.interscaling,flags.scaling);

if isempty(L)
   L = wfbtlength(max(cellfun(@(gEl) numel(gEl.h),gmultid)),wt);
end

% Do the equivalent uniform filterbank
if any(amultid~=amultid(1))
   [gu,au] = nonu2ufilterbank(gmultid,amultid);
else
   [gu,au] = deal(gmultid,amultid);
end

if nargout<2
   AF = filterbankbounds(gu,au,L);
elseif nargout == 2
   [AF, BF] = filterbankbounds(gu,au,L);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%io_writepta.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% RF=io_writepta(rf,outfile);
%
% DESCRIPTION:
% Write a matlab RF pulse structure (containing N x 3 waveform array field
% with rf.waveform(:,1)= phase, rf.waveform(:,2)=amplitude, and rf.waveform(:,3)=timestep),
% to a siemens format .pta file.
%
% INPUTS:
% rf         = matlab RF pulse.
% outfile    = name of the output .pta file to be written.
%
% OUTPUTS:
% RF            = Same as input.  Not used.  The primary output of this
%                   function is a text file in Siemens .pta format.

function RF=io_writepta(rf,outfile);

B1INT=sum(rf.waveform(:,2));

l=length(rf.waveform(:,2));
RF.waveform=zeros(3,l);
RF.waveform(1,:)=rf.waveform(:,2)/max(rf.waveform(:,2));
RF.waveform(2,:)=mod(rf.waveform(:,1)*pi/180,2*pi);
RF.waveform(3,:)=[0:l-1];

%For simplicity:  Make a complex valued version of the RF pulse
RF_complex=RF.waveform(1,:).*exp(1i*RF.waveform(2,:));

%Calculate the POWER INTEGRAL (POWERINT).  This is used for calculating the
%SAR:
POWERINT=sum((real(RF_complex).^2)+(imag(RF_complex).^2));

%Calculate the MAGNITUDE/ABSOLUTE INTEGRAL (ABSINT):
ABSINT=sum(abs(RF_complex));

%AMPLITUDE INTEGRAL (AMPINT)...  This is the parameter that is used to
%calculate the transmitter voltage at runtime, and therefore it determines
%the flip angle that will be realied.  It must be calculated properly to
%ensure that the flip angle will be correct.  In addition, the calculation
%of AMPINT depends on the type of pulse that you are dealing with.  For a
%straigt plain wave (That is, the phase of each data point is either a
%constant ph, or ph+180 or ph-180)  the AMPINT is simply the sum of the
%waveform points when the maximum value is scaled to 1.  For example, for a
%three point pulse with point values of -1, 3 and -1, the AMPINT would be
%1.  This will be the method used for automatic calculation of AMPINT.  If
%the pulse is not a plane pulse, do not use this method to calculate
%AMPINT.

AIcalc=input('Would you like to perform an automatic Amplitude Inegral Calculation?(y or n)  ', 's');

if AIcalc=='y' || AIcalc=='Y'
    AMPINT=sqrt((sum(real(RF_complex)).^2)+(sum(imag(RF_complex)).^2));
elseif AIcalc=='n' || AIcalc=='N'
    % If the pulse is not a plane wave (Adiabatic pulses, off resonance pulses,
    % or dual Band Pulses), you should calculate the Amplitude integral
    % yourself and input it manually.  For methods of calculating AMPINT in this
    % case, see the document entitled
    % "RF_amplitude_integral_in_pulse_sequences.pdf', which can be found in
    % /Users/jnear/Documents/RF pulses/
    AMPINT=input('Now you must input the Amplitude Integral Manually.  Lets have it?  ');
else
    error('ERROR!  You must input a value for AMPINT!  ABORTING!');
end

%Reference gradient is only used in slice selective pulses.  This refers to
%the gradient required to excite a 10mm slice if the pulse duration is
%5.12ms  Given that the Time-bandwidth product is included in the rf
%structure, the refgrad value can be calculated here, but the option is
%provided to enter it manually, in case a different refgrad value is
%desired.

refgrad=1;
refg=input('Would you like to input a reference gradient (REFGRAD)? ','s');
if refg=='y' || refg=='Y'
    disp('INFORMATION:  REFGRAD is the gradient amplitude in [mT/m] required ');
    disp('to excite a 10mm slice with a pulse duration of 5.12 ms. ');
    man=input('Calculate refgrad automatically (y or n) ? ','s');
    if man=='y' || man=='Y'
        %refgrad [mT/m] = (1000 [mT]/[1T]) * tbw [unitless] ) / ( (0.00512 [s]) * (42577000 [Hz]/[T]) * (0.01 [m]) );
        refgrad = 1000 * rf.tbw/(0.00512 * 42577000 * 0.01);
    elseif man=='n' || man=='N'
        refgrad=input('Manually enter the reference gradient (mT/m):  ');
    else
        error('ERROR:  Response not recognized!!');
    end
end

%write to pta file for siemens
fid=fopen(outfile,'w+');
fprintf(fid,['PULSENAME: \t' outfile]);
fprintf(fid,'\nCOMMENT:  \tRF pulse generated using the FID-A toolkit (github.com/CIC-methods/FID-A).');
fprintf(fid,'\nREFGRAD:  \t%5.6f' ,refgrad);
fprintf(fid,'\nMINSLICE: \t1.00000000');
fprintf(fid,'\nMAXSLICE: \t200.00000000');
fprintf(fid,'\nAMPINT: \t%5.6f',AMPINT);
fprintf(fid,'\nPOWERINT: \t%5.6f',POWERINT);
fprintf(fid,'\nABSINT: \t%5.6f',ABSINT);
fprintf(fid,'\n\n');
n=0;
fprintf(fid,'%1.9f %1.9f ; (%1.0f)\n',RF.waveform);
fclose(fid);

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


%% ANAP-building
% %{
imfolder = 'images\ANAP-building';
im_n = 2;
imfile = cell(im_n,1);
imfile{1} = [imfolder '\' 'building0.jpg'];
imfile{2} = [imfolder '\' 'building1.jpg'];

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = [1,2];

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 1, 'persp', 200, imfolder );
%}


text:
function a=update_a(dvec,n,alpha)
% a=update_a(dvec,n, alpha)
% Computes update for a_ki (KxI) in variational approximation (Buntine & Jakulin
% DCA 2006)
% dvec: data - each column is one image (J x I)
% n: (J x K x I)
% alpha, beta: parameters of the Gamma distribution for latent h (1xK each)
% J=#pixels
% K=#components
% I=#images.

[j,k,i]=size(n);

dvec_reshaped = reshape(dvec,j,1,i); % (Jx1xI)
atmp = sum(bsxfun(@times,dvec_reshaped,n),1); % (1xKxI)
a=bsxfun(@plus, squeeze(atmp), alpha'); % (KxI)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function btv_test01 ( )

%*****************************************************************************80
%
%% BTV_TEST01 tests BURGERS_TIME_VISCOUS with the gaussian initial condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTV_TEST01\n' );
  fprintf ( 1, '  Test BURGERS_TIME_VISCOUS with the gaussian initial condition.\n' );

  nx = 81;
  nt = 200;
  t_max = 2.0;
  nu = 0.01;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial condition: gaussian\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Viscosity = %g\n', nu );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_viscous ( @ic_gaussian, nx, nt, t_max, nu, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 1 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers equation solutions over time, initial condition gaussian' )

  filename = 'btv_test01.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end


text:
%% test_all.m
%|
%| Attempt to run all test routines
%| to verify completeness of software distribution.
%|
%| I recommend running these tests to verify the completeness
%| of your installation!  If some test fails, please email me
%| the relevant output messages and tell me what version of Matlab
%| and what OS you are using.  Then go ahead and try to use the
%| toolbox because often most things work even if one test fails.
%|
%| If you do not have all the matlab toolboxes (e.g., image, signal, wavelet)
%| then some tests will fail.  But don't worry about it, I am trying to
%| make it work with just plain matlab, so most of it should work even
%| without any extra toolboxes.
%|
%| Jeff Fessler

%double6 double
%disp 'Note: forcing double precision because sparse(double)*single fails :-('

im off-quiet

prompt draw % do all plots without pausing/prompting
%printm 'Note: hit "r" at the prompt to disable subsequent prompts'

diary_file = ['~/' date '-test_all.txt'];
diary(diary_file)

% uncomment the suite(s) you want to test

test_all_mex % see if all the mex files work (this will fail on Windows)

list = {...
'test_all_util', ...
'test_all_graph', ...
'test_all_reg', ...
'test_all_nufft', ...
'test_all_systems', ...
'test_all_emission', ...
'test_all_transmission', ...
'test_all_wls', ...
'test_all_ct', ...
'test_all_mri', ...
'test_all_example'
};

run_mfile_local(list, 'pause', 0)

diary off


text:
function [roi_table, region_obj] = ttest_table_by_condition(region_obj, DATA_OBJ, varargin)
% Extract region-by-region data from DATA_OBJ fmri_data object(s) in
% regions specified by region_obj, and return a table of means and stats,
% along with extracted data.
%
% First line: One-line summary description of function
%
% :Usage:
% ::
%
%     [roi_table, region_obj] = ttest_table_by_condition(region_obj, DATA_OBJ, [optional inputs])
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2017 Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **region_obj:**
%        A region object from CANlab core tools
%
%   **DATA_OBJ:**
%        an fmri_data object or cell vector of fmri_data objects with data
%        to extract.  Does not need to be in the same space as region
%        object -- extract_data.m method converts using mm coordinates.
%        Reslicing can induce variation in results due to interpolation.
%
% :Optional Inputs:
%   **'conditionnames', 'condition_names':**
%        Followed by cell vector of names for each condition
%        No spaces or special characters.
%
%
% :Outputs:
%
%   **roi_table:**
%        A table object containing ROI coordinates, names, and indices for
%        each region, along with means, and t-stats
%
%   **region_obj:**
%        The region object you entered in, with appended data,
%        images x conditions, stored in .dat field for each region.
%
% :Examples:
% ::
%
%    % give examples of code here
%    param1 = abc();
%    param2 = xyz();
%    [out1,out2] = func_call(param1, param2)
%
% :References:
%   CITATION(s) HERE
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
% ..

% -------------------------------------------------------------------------
% DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

% Get volume from a region
volfunction = @(x) x.numVox .* prod(abs(diag(x.M(1:3, 1:3))));

if ~iscell(DATA_OBJ)
    DATA_OBJ = {DATA_OBJ};
end

nconditions = length(DATA_OBJ);
conditionnames = cell(1, nconditions);

for i = 1:length(DATA_OBJ)
    conditionnames{i} = sprintf('Cond_%d', i);
end

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case {'conditionnames', 'condition_names', 'conditions'}
                conditionnames = varargin{i+1}; varargin{i+1} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

if length(conditionnames) ~= nconditions
    error('Condition names and DATA_OBJ must be cells of same length');
end

nrois = length(region_obj);

% Initial table
% ------------------------------------------------------------------------
Index = (1:nrois)';
Name_of_Region = char(region_obj.shorttitle);
xyz_mm = cat(1, region_obj.mm_center);

for i = 1:nrois
    volume_mm3(i, 1) = volfunction(region_obj(i));
end

roi_table = table(Index, Name_of_Region, xyz_mm, volume_mm3);

% Initial table
% ------------------------------------------------------------------------

conditionnames = strrep(conditionnames, ' - ', '_vs_');
conditionnames = strrep(conditionnames, '-', '_vs_');
conditionnames = strrep(conditionnames, ' ', '_');

% extract ROI average data from each region
% Make regions x conditions matrix
% ------------------------------------------------------------------------

stars_by_condition = cell(nrois, nconditions);
[means_by_condition, t_by_condition, p_by_condition] = deal(zeros(nrois, nconditions));

for i = 1:nconditions

    rois_tmp = extract_data(region_obj, DATA_OBJ{i});

    % cat within regions to return data in region_obj object
    clear mydat
    for j = 1:nrois

        mydat = rois_tmp(j).dat(:, 1);

        % pad with longer of two
        [Anew Bnew] = padwithnan(region_obj(j).dat, mydat, 1);

        % region_obj(j).dat(:, i)
        region_obj(j).dat = [Anew Bnew];

    end

    % cat into image x roi matrix
    % for statistical tests
    dat = cat(2, rois_tmp.dat);

    % t-tests on each region across images (often subjects)
    [h, p, ci, stat] = ttest(dat);

    means_by_condition(:, i) = nanmean(dat)';

    t_by_condition(:, i) = stat.tstat';

    p_by_condition(:, i) = p';

    % Stars for each region
    for j = 1:length(p)

        if p(j) < .0015, mystr = '***';
        elseif p(j) < .015, mystr = '**';
        elseif p(j) < .055, mystr = '*';
        elseif p(j) < .105, mystr = '+';
        else mystr = ''; xadj = 0;
        end

        stars_by_condition{j, i} = mystr;

    end % loop through regions

    roi_table.(['t_' conditionnames{i}]) =  t_by_condition(:, i);

    signame = sprintf('sig%d', i);
    roi_table.(signame) = stars_by_condition(:, i);

end % loop through conditions

roi_table.means_by_condition = means_by_condition;
roi_table.p_by_condition = p_by_condition;

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function bk = shepard_basis_1d ( nd, xd, k, p, ni, xi )

%*****************************************************************************80
%
%% SHEPARD_BASIS_1D evaluates a 1D Shepard basis function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Shepard,
%    A two-dimensional interpolation function for irregularly spaced data,
%    ACM '68: Proceedings of the 1968 23rd ACM National Conference,
%    ACM, pages 517-524, 1969.
%
%  Parameters:
%
%    Input, integer ND, the number of data points.
%
%    Input, real XD(ND,1), the data points.
%
%    Input, integer K, the index of the desired basis function,
%    1 <= K <= ND.
%
%    Input, real P, the power.
%
%    Input, integer NI, the number of interpolation points.
%
%    Input, real XI(NI,1), the interpolation points.
%
%    Output, real BK(NI,1), the basis function at the interpolation points.
%
  bk = zeros ( ni, 1 );

  for i = 1 : ni

    if ( p == 0.0 )

      w(1:nd,1) = 1.0 / nd;

    else

      w = zeros ( nd, 1 );

      z = -1;
      for j = 1 : nd
        w(j) = abs ( xi(i) - xd(j) );
        if ( w(j) == 0.0 )
          z = j;
          break
        end
      end

      if ( z ~= -1 )
        w = zeros ( nd, 1 );
        w(z) = 1.0;
      else
        w(1:nd,1) = 1.0 ./ w(1:nd,1) .^ p;
        s = sum ( w );
        w(1:nd,1) = w(1:nd,1) / s;
      end

    end

    bk(i,1) = w(k,1);

  end

  return
end


text:
% CVX: Keywords for declaring variables and objectives
%   dual      - Declares one or more dual variables.
%   epigraph  - Declares an epigraph variable.
%   hypograph - Declares a hypograph variable.
%   maximize  - Specifies a concave (or affine) objective to be maximized.
%   minimize  - Specifies a convex (or affine) objective to be maximized.
%   subject   - Implements the "subject to" keyword.
%   variable  - Declares a single CVX variable with optional matrix structure.
%   variables - Declares one or more CVX variables.

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [R] = F2R(F)
% Convert temperature from degrees Fahrenheit to Rankine.
R = F+459.67;

author: Sable

repo: mcbench-benchmarks


text:
function Angles = FindAngles(vertices)

% Purpose: To find the angles of Triangle or Quadrilateral
% Synopsis: Angles = FindAngles(vertices)
% Variable description:
% INPUT: vertices - coordinates of Triangle and Quadrilateral
%                   order 3X2 if vertices are of Triangle and 4X2 if vertices
%                   are of Quadrilateral. The vertices are given input in
%                   anticlock wise direction.
% OUTPUT: Angles - Row vector of angles of Triangle/ Quadrilateral in degrees.
%                  3X1 for Trianle and 4X1 for Quadrilateral
% Example: vertices = [4 0;5 1;4 1]                 } Triangle vertices
%          angles = FindAngles(vertices);
%          angles = [45.0000   45.0000   90.0000]
%
%          vertices = [4  4;4  6;5 6; 5  4]         } Quadrilateral vertices
%          angles = FindAngles(vertices);
%          angles = [90 90 90 90]
%
% Algorithm: Calculates the angles using dot product
%            theta = arccos((a.b)/(|a||b|)
% Release : 1.0
% Release Date: 03/07/2012
%---------------------------------------------------------------------------
% Author : Siva Srinivas Kolukula
%          Senior Research Fellow
%          Structural Mechanics Laboratory
%          Indira Gandhi Center for Atomic Research
%          India
% E-mail : allwayzitzme@gmail.com
%          http://sites.google.com/site/kolukulasivasrinivas/
%---------------------------------------------------------------------------

N = length(vertices) ;
if N == 3
    disp('vertices are of Triangle')
    order = [1 2 3; 2 1 3;3 1 2] ;
elseif N == 4
    disp('Vertices are of Quadrilateral')
    order = [ 1 2 4; 2 1 3; 3 4 2;4 3 1] ;
else
    error('Wrong input')
end

    Angles = zeros(1,N) ;

for i = 1:N
    % Form the vectors along the vertices
    vec1 = [vertices(i,:); vertices(order(i,2),:)];
    vec2 = [vertices(i,:); vertices(order(i,3),:)];
    % Use dor ptodict to find the angle between vector1 and vector2
    Angles(i) = acos((diff(vec1)*diff(vec2)')/.....
        (sqrt(sum(diff(vec1,[],1).^2,2))*sqrt(sum(diff(vec2,[],1).^2,2)))) ;
end

% Convert the angles in radians to degrees
Angles = Angles*180./pi ;        %


text:
% This file is part of TREEQSM.
%
% TREEQSM 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.
%
% TREEQSM 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 TREEQSM.  If not, see <http://www.gnu.org/licenses/>.

function pmdistance = point_model_distance(P,cylinder)

% ---------------------------------------------------------------------
% POINT_MODEL_DISTANCE.M    Computes the distances of the points to the
%                               cylinder model
%
% Version 2.1.1
% Latest update     8 Oct 2021
%
% Copyright (C) 2015-2021 Pasi Raumonen
% ---------------------------------------------------------------------

% Changes from version 2.1.0 to 2.1.1, 8 Oct 2021:
% 1) Changed the determinationa NE, the number of empty edge layers, so
%     that is now limited in size, before it is given as input for
%     cubical_partition function.

% Changes from version 2.0.0 to 2.1.0, 26 Nov 2019:
% 1) Bug fix: Corrected the computation of the output at the end of the
%    code so that trees without branches are computed correctly.

% Cylinder data
Rad = cylinder.radius;
Len = cylinder.length;
Sta = cylinder.start;
Axe = cylinder.axis;
BOrd = cylinder.BranchOrder;

% Select randomly 25 % or max one million points for the distance comput.
np0 = size(P,1);
a = min(0.25*np0,1000000);
I = logical(round(0.5/(1-a/np0)*rand(np0,1)));
P = P(I,:);

% Partition the points into cubes
L = 2*median(Len);
NE = max(3,min(10,ceil(max(Len)/L)))+3;
[Partition,~,Info] = cubical_partition(P,L,NE);
Min = Info(1:3);
EL = Info(7);
NE = Info(8);

% Calculates the cube-coordinates of the starting points
CC = floor([Sta(:,1)-Min(1) Sta(:,2)-Min(2) Sta(:,3)-Min(3)]/EL)+NE+1;

% Compute the number of cubes needed for each starting point
N = ceil(Len/L);

% Correct N so that cube indexes are not too small or large
I = CC(:,1) < N+1;
N(I) = CC(I,1)-1;
I = CC(:,2) < N+1;
N(I) = CC(I,2)-1;
I = CC(:,3) < N+1;
N(I) = CC(I,3)-1;
I = CC(:,1)+N+1 > Info(4);
N(I) = Info(4)-CC(I,1)-1;
I = CC(:,2)+N+1 > Info(5);
N(I) = Info(5)-CC(I,2)-1;
I = CC(:,3)+N+1 > Info(6);
N(I) = Info(6)-CC(I,3)-1;

% Calculate the distances to the cylinders
n = size(Rad,1);
np = size(P,1);
Dist = zeros(np,2); % Distance and the closest cylinder of each points
Dist(:,1) = 2; % Large distance initially
Points = zeros(ceil(np/10),1,'int32'); % Auxiliary variable
Data = cell(n,1);
for i = 1:n
  Par = Partition(CC(i,1)-N(i):CC(i,1)+N(i),CC(i,2)-N(i):CC(i,2)+N(i),...
    CC(i,3)-N(i):CC(i,3)+N(i));
  if N(i) > 1
    S = cellfun('length',Par);
    I = S > 0;
    S = S(I);
    Par = Par(I);
    stop = cumsum(S);
    start = [0; stop]+1;
    for k = 1:length(stop)
      Points(start(k):stop(k)) = Par{k}(:);
    end
    points = Points(1:stop(k));
  else
    points = vertcat(Par{:});
  end
  [d,~,h] = distances_to_line(P(points,:),Axe(i,:),Sta(i,:));
  d = abs(d-Rad(i));
  Data{i} = [d h double(points)];
  I = d < Dist(points,1);
  J = h >= 0;
  K = h <= Len(i);
  L = d < 0.5;
  M = I&J&K&L;
  points = points(M);
  Dist(points,1) = d(M);
  Dist(points,2) = i;
end

% Calculate the distances to the cylinders for points not yet calculated
% because they are not "on side of cylinder
for i = 1:n
  if ~isempty(Data{i})
    d = Data{i}(:,1);
    h = Data{i}(:,2);
    points = Data{i}(:,3);
    I = d < Dist(points,1);
    J = h >= -0.1 & h <= 0;
    K = h <= Len(i)+0.1 & h >= Len(i);
    L = d < 0.5;
    M = I&(J|K)&L;
    points = points(M);
    Dist(points,1) = d(M);
    Dist(points,2) = i;
  end
end

% Select only the shortest 95% of distances for each cylinder
N = zeros(n,1);
O = zeros(np,1);
for i = 1:np
  if Dist(i,2) > 0
    N(Dist(i,2)) = N(Dist(i,2))+1;
    O(i) = N(Dist(i,2));
  end
end
Cyl = cell(n,1);
for i = 1:n
  Cyl{i} = zeros(N(i),1);
end
for i = 1:np
  if Dist(i,2) > 0
    Cyl{Dist(i,2)}(O(i)) = i;
  end
end
DistCyl = zeros(n,1); % Average point distance to each cylinder
for i = 1:n
  I = Cyl{i};
  m = length(I);
  if m > 19 % select the smallest 95% of distances
    d = sort(Dist(I,1));
    DistCyl(i) = mean(d(1:floor(0.95*m)));
  elseif m > 0
    DistCyl(i) = mean(Dist(I,1));
  end
end

% Define the output
pmdistance.CylDist = single(DistCyl);
pmdistance.median = median(DistCyl(:,1));
pmdistance.mean = mean(DistCyl(:,1));
pmdistance.max = max(DistCyl(:,1));
pmdistance.std = std(DistCyl(:,1));

T = BOrd == 0;
B1 = BOrd == 1;
B2 = BOrd == 2;
B = DistCyl(~T,1);
T = DistCyl(T,1);
B1 = DistCyl(B1,1);
B2 = DistCyl(B2,1);

pmdistance.TrunkMedian = median(T);
pmdistance.TrunkMean = mean(T);
pmdistance.TrunkMax = max(T);
pmdistance.TrunkStd = std(T);

if ~isempty(B)
  pmdistance.BranchMedian = median(B);
  pmdistance.BranchMean = mean(B);
  pmdistance.BranchMax = max(B);
  pmdistance.BranchStd = std(B);
else
  pmdistance.BranchMedian = 0;
  pmdistance.BranchMean = 0;
  pmdistance.BranchMax = 0;
  pmdistance.BranchStd = 0;
end

if ~isempty(B1)
  pmdistance.Branch1Median = median(B1);
  pmdistance.Branch1Mean = mean(B1);
  pmdistance.Branch1Max = max(B1);
  pmdistance.Branch1Std = std(B1);
else
  pmdistance.Branch1Median = 0;
  pmdistance.Branch1Mean = 0;
  pmdistance.Branch1Max = 0;
  pmdistance.Branch1Std = 0;
end

if ~isempty(B2)
  pmdistance.Branch2Median = median(B2);
  pmdistance.Branch2Mean = mean(B2);
  pmdistance.Branch2Max = max(B2);
  pmdistance.Branch2Std = std(B2);
else
  pmdistance.Branch2Median = 0;
  pmdistance.Branch2Mean = 0;
  pmdistance.Branch2Max = 0;
  pmdistance.Branch2Std = 0;
end

author: InverseTampere

repo: TreeQSM


text:
function ns_RESULT = ns_CloseFile(hFile);

%ns_CloseFile   Closes a neural data file
%
%   Usage:
%      ns_RESULT = ns_CloseFile(hFile)
%
%   Description:
%       Closes a previously opened file specified by the file handle hFile.
%
%   Parameters:
%       hFile	    Handle/Indentification number to an open file.
%
%   Return Values:
%       ns_RESULT   This function returns ns_OK if the file is successfully
%                   opened. Otherwise one of the following error codes is
%                   generated:
%
%                       ns_BADFILE	Invalid file handle passed to function.
%
%   Copyright (C) 2003 Neuroshare Project
%   Author: Almut Branner
%   Last modification: 2/4/2003

ns_RESULT = mexprog(14, hFile);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function fxk = Rate_function(opt, varargin)
% USAGE:
%
%    Rate_function(opt, varargin)
%
% INPUTS:
%    xk:     current point;
%    opt:    structure includes required parameters;
%
%              * .FR - concatenation of forward and reverse stoichiometric matrix
%              * .A - Reduced forward stoichiometric matrix
%              * .B - Reduced reverse stoichiometric matrix
%              * .k - initial kinetic
%
% OUTPUT:
%    fxk:    the vector `f(xk)`

if nargin ~= 2
    error('The number of input arguments is not valid');
end

if nargout >= 2
    error('The number of output arguments is not valid');
end

FR    = opt.FR;
FR_RF = opt.FR_RF;
k     = opt.k;
xk    = varargin{1};

% ======================== Function evaluation =========================

fxk = FR_RF*exp(k+FR'*xk);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% End of Rate_function.m %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
function fiber_torsion_vals=dtiFiberTorsion(fiber)

%Compute fiber torsion (at every segment long the fiber)
%For 3xn fiber as defined here
%http://en.wikipedia.org/wiki/Torsion_of_curves

%
%ER 12/2007
xprime=gradient(fiber(1, :));
yprime=gradient(fiber(2, :));
zprime=gradient(fiber(3, :));

xpp=gradient(xprime);
ypp=gradient(yprime);
zpp=gradient(zprime);

fiber_torsion_vals=(gradient(zpp).*(xprime.*ypp-yprime.*xpp)+zpp.*(gradient(xpp).*yprime-xprime.*gradient(ypp))+zprime.*(xpp.*gradient(ypp)-gradient(xpp).*ypp))./((xprime.^2+yprime.^2+zprime.^2).*(xpp.^2+ypp.^2+zpp.^2))

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function cmap = cmapExtendedHSV(numGrays,numColors,range)
%
% cmap = cmapExtendedHSV([numGrays=128],[numColors=96],[range=query user])
%
%The map created here are hsv maps where all of the colors
% can be placed in a subsection of the full color map.  In this way, the
% full range of colors spans less than 2pi, like in the double color map.
% Rather than compressing by a complete factor of 2, like the double color
% map, the compression factor can be a bit smaller.
%
%   There are numGrays gray scale entries.  They occupy the first part of the cmap, 1:numGrays
%   There are numColors hsv colors.  They fill the map entries following the gray,
%   from (numGrays+1):numGrays+numColors
%   When the range is 1, the hsv map is hsv(numColors) and we insert it
%   into the cmap.
%   When the range is 1.5, we compute tmp = hsv(numColors/1.5) and we
%   create [tmp,tmp(1:needed)] to fill up numColors entries.
%
% Examples:
%   cmap = cmapExtendedHSV(128,96,1.1);
%   cmap = cmapExtendedHSV(128,96);    -- 128 gray levels, 96 color levels,
%          query use for compression
%   cmap = cmapExtendedHSV;
%   cmap = cmapExtendedHSV(128,96,2);  -- Same as hsvDoubleCmap
%

function vw = cmapExtended(vw,range)
%
%   vw = cmapExtended(vw,range)
%
% Author: AAB, BW
% Purpose:
%   Compress the current color map to smaller range and add gray at the end.
%
%  vw = FLAT{1};

numGrays = viewGet(vw,'cmapcurnumgrays');
numColors = viewGet(vw,'cmapcurnumcolors');
mpColors = viewGet(vw,'cmapcurrent');

if ieNotDefined('vw'),  error('Must pass in the view.'); end
if ieNotDefined('range'), range = 1.2; end

if (range < 1) | (range > 2),  error('Range must be betweem 1 and 2.'); end

desiredColors = round(numColors/range);
extraColors = numColors - desiredColors;
newMap = [round(interp1(mpColors',(1:range:numColors)')')];
newMap = [newMap,ones(3,extraColors)*128];

vw = viewSet(vw,'cmapcurrent',newMap);

% cmap = zeros(numGrays+numColors,3);
%
% % If you want the map symmetric at the boundary, you should do this.
% % We could trap range == 2 and do it then ... which would be backwards
% % compatible?
% % cmap = [gray(numGrays); hsvMap; flipud(hsvMap(1:hsvColorsExtra,:))];
% cmap = [hsvMap; hsvMap(1:hsvColorsExtra,:)];
% shiftSize = round(hsvColorsExtra/2);
% hsvMap = circshift(cmap,shiftSize);
%
% cmap = [gray(numGrays); cmap];

return;

%----------------------------------------
function range = readRange

prompt={'Enter compression range for the hsv map (2 = double color map)'};
def={'1.2'};
dlgTitle='Color map compression factor';
lineNo=1;
range=inputdlg(prompt,dlgTitle,lineNo,def);
range = str2num(range{1});

return;


text:
function [ n_data, x, y, r, fxy ]  = bivariate_normal_cdf_values ( n_data )

%*****************************************************************************80
%
%% BIVARIATE_NORMAL_CDF_VALUES returns some values of the bivariate normal CDF.
%
%  Discussion:
%
%    FXY is the probability that two variables A and B, which are
%    related by a bivariate normal distribution with correlation R,
%    respectively satisfy A <= X and B <= Y.
%
%    Mathematica can evaluate the bivariate normal CDF via the commands:
%
%      <<MultivariateStatistics`
%      cdf = CDF[MultinormalDistribution[{0,0}{{1,r},{r,1}}],{x,y}]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 November 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    National Bureau of Standards,
%    Tables of the Bivariate Normal Distribution and Related Functions,
%    NBS, Applied Mathematics Series, Number 50, 1959.
%
%  Parameters:
%
%    Input, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.
%
%    Output, integer N_DATA, the routine increments the input value of N_DATA
%    by 1, and returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, Y, the parameters of the function.
%
%    Output, real R, the correlation value.
%
%    Output, real FXY, the value of the function.
%
  n_max = 41;

  fxy_vec(1:n_max) = [ ...
  0.02260327218569867, ...
  0.1548729518584100, ...
  0.4687428083352184, ...
  0.7452035868929476, ...
  0.8318608306874188, ...
  0.8410314261134202, ...
  0.1377019384919464, ...
  0.1621749501739030, ...
  0.1827411243233119, ...
  0.2010067421506235, ...
  0.2177751155265290, ...
  0.2335088436446962, ...
  0.2485057781834286, ...
  0.2629747825154868, ...
  0.2770729823404738, ...
  0.2909261168683812, ...
  0.3046406378726738, ...
  0.3183113449213638, ...
  0.3320262544108028, ...
  0.3458686754647614, ...
  0.3599150462310668, ...
  0.3742210899871168, ...
  0.3887706405282320, ...
  0.4032765198361344, ...
  0.4162100291953678, ...
  0.6508271498838664, ...
  0.8318608306874188, ...
  0.0000000000000000, ...
  0.1666666666539970, ...
  0.2500000000000000, ...
  0.3333333333328906, ...
  0.5000000000000000, ...
  0.7452035868929476, ...
  0.1548729518584100, ...
  0.1548729518584100, ...
  0.06251409470431653, ...
  0.7452035868929476, ...
  0.1548729518584100, ...
  0.1548729518584100, ...
  0.06251409470431653, ...
  0.6337020457912916 ];
  r_vec(1:n_max) = [ ...
     0.500,  0.500,  0.500,  0.500,  0.500, ...
     0.500, -0.900, -0.800, -0.700, -0.600, ...
    -0.500, -0.400, -0.300, -0.200, -0.100, ...
     0.000,  0.100,  0.200,  0.300,  0.400, ...
     0.500,  0.600,  0.700,  0.800,  0.900, ...
     0.673,  0.500, -1.000, -0.500,  0.000, ...
     0.500,  1.000,  0.500,  0.500,  0.500, ...
     0.500,  0.500,  0.500,  0.500,  0.500, ...
     0.500 ];
  x_vec(1:n_max) = [ ...
    -2.0, -1.0,  0.0,  1.0,  2.0, ...
     3.0, -0.2, -0.2, -0.2, -0.2, ...
    -0.2, -0.2, -0.2, -0.2, -0.2, ...
    -0.2, -0.2, -0.2, -0.2, -0.2, ...
    -0.2, -0.2, -0.2, -0.2, -0.2, ...
     1.0,  2.0,  0.0,  0.0,  0.0, ...
     0.0,  0.0,  1.0,  1.0, -1.0, ...
    -1.0,  1.0,  1.0, -1.0, -1.0, ...
     0.7071067811865475  ];
  y_vec(1:n_max) = [ ...
     1.0,  1.0,  1.0,  1.0,  1.0, ...
     1.0,  0.5,  0.5,  0.5,  0.5, ...
     0.5,  0.5,  0.5,  0.5,  0.5, ...
     0.5,  0.5,  0.5,  0.5,  0.5, ...
     0.5,  0.5,  0.5,  0.5,  0.5, ...
     0.5,  1.0,  0.0,  0.0,  0.0, ...
     0.0,  0.0,  1.0, -1.0,  1.0, ...
    -1.0,  1.0, -1.0,  1.0, -1.0, ...
     0.7071067811865475 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    r = 0.0;
    x = 0.0;
    y = 0.0;
    fxy = 0.0;
  else
    r = r_vec(n_data);
    x = x_vec(n_data);
    y = y_vec(n_data);
    fxy = fxy_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef StatusState
    %STATUSSTATE Represent status selections
    %   Enumerates a list of choices

    % Copyright 2020-2021 The MathWorks Inc.

    %% Enumerations
    enumeration
        none
        complete
        running
        warning
        error
        info
        unknown
    end %enumeration

end %classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function x = p02_x ( n )

%*****************************************************************************80
%
%% P02_X returns the least squares solution X for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Output, real X(N,1), the least squares solution.
%
  x = [ 5.7013; 121.1341; 152.4745 ];

  return
end