text:
function [warped]= sn2individual(P, input)

% SN2INDIVIDUAL warps the input coordinates (defined as Nx3 matrix) from
% normalised MNI coordinates to individual headspace coordinates, using the
% warp parameters defined in the structure spmparams.
%
% this is modified from code from nutmeg: nut_mni2mri, which was itself
% modified from code originally written by John Ashburner:
% http://www.sph.umich.edu/~nichols/JG2/get_orig_coord2.m

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

if isfield(P, 'Tr')
  % this is an old-style representation of the parameters, so it uses the
  % code adjusted from nut_mni2mri

  if numel(P.Tr)==0
    % only an affine transformation has been done
    T      = P.VF.mat*P.Affine/(P.VG.mat);
    warped = ft_warp_apply(T, input);

  else
    % we need the spm_dctmtx function for the nonlinear case
    if ~ft_hastoolbox('spm')
      % add SPM8 or later to the path
      ft_hastoolbox('spm8up', 1);
    end

    dim  = P.VG.dim(1:3);
    xyz  = ft_warp_apply(inv(P.VG.mat), input); % goes into voxel coordinates

    basX = spm_dctmtx(dim(1), size(P.Tr,1), xyz(:,1)-1);
    basY = spm_dctmtx(dim(2), size(P.Tr,2), xyz(:,2)-1);
    basZ = spm_dctmtx(dim(3), size(P.Tr,3), xyz(:,3)-1);

    siz = size(P.Tr);
    Tr1 = reshape(P.Tr(:,:,:,1),siz(1)*siz(2),siz(3));
    Tr2 = reshape(P.Tr(:,:,:,2),siz(1)*siz(2),siz(3));
    Tr3 = reshape(P.Tr(:,:,:,3),siz(1)*siz(2),siz(3));

    xyztmp = zeros(size(xyz));
    for i=1:size(xyz,1)
      bx = basX(i,:);
      by = basY(i,:);
      bz = basZ(i,:);
      tx = reshape(Tr1*bz', siz(1), siz(2) );
      ty = reshape(Tr2*bz', siz(1), siz(2) );
      tz = reshape(Tr3*bz', siz(1), siz(2) );
      xyztmp(i,:) = [bx*tx*by' bx*ty*by' bx*tz*by'];
    end

    T      = P.VF.mat*P.Affine;
    warped = ft_warp_apply(T, xyz+xyztmp);
  end

else
  % the only way I can come up with to do this, is to write a deformation
  % to disk, and to sample this one. This requires spm12 on the path
  ft_hastoolbox('spm12', 1);

  fprintf('creating the deformation field and writing it to a temporary file\n');

  fname  = [tempname,'.nii'];
  V      = nifti;
  V.dat  = file_array(fname, P.image.dim(1:3), [spm_type('float32') spm_platform('bigend')], 0, 1, 0);
  V.mat  = P.image.mat;
  if isfield(P.image, 'mat0')
    V.mat0 = P.image.mat0;
  end
  V.descrip = 'deformation field';
  create(V);
  V.dat(:) = 0; % this is necessary, otherwise SPM fails: image too small

  P.image  = spm_vol(fname);
  spm_preproc_write8(P, zeros(6,4), [0 0], [0 1], 1, 1, nan(2,3), nan);

  [pth,nam,ext] = fileparts(fname);
  V  = nifti(fullfile(pth,['y_',nam,ext]));
  y  = squeeze(V.dat(:,:,:,:,:));

  siz = size(y);
  VT.dim = siz(1:3);
  VT.mat = P.tpm(1).mat;

  % 2b: write the deformation fields in x/y/z direction to temporary files
  V1.fname     = [tempname '.img'];
  V1.dim(1:3)  = VT.dim(1:3);
  V1.pinfo     = [1 0 0]';
  V1.mat       = VT.mat;
  V1.dt        = [64 0];
  V1.descrip   = 'Deformation field';
  spm_write_vol(V1,y(:,:,:,1));

  V2.fname     = [tempname '.img'];
  V2.dim(1:3)  = VT.dim(1:3);
  V2.pinfo     = [1 0 0]';
  V2.mat       = VT.mat;
  V2.dt        = [64 0];
  V2.descrip   = 'Deformation field';
  spm_write_vol(V2,y(:,:,:,2));

  V3.fname     = [tempname '.img'];
  V3.dim(1:3)  = VT.dim(1:3);
  V3.pinfo     = [1 0 0]';
  V3.mat       = VT.mat;
  V3.dt        = [64 0];
  V3.descrip   = 'Deformation field';
  spm_write_vol(V3,y(:,:,:,3));

  % first warp to voxel coordinates
  input_vox = ft_warp_apply(inv(VT.mat), input);  % Express as voxel indices

  % apply the non-linear warp
  warped = cat(2, spm_sample_vol(V1,input_vox(:,1),input_vox(:,2),input_vox(:,3),1), ...
    spm_sample_vol(V2,input_vox(:,1),input_vox(:,2),input_vox(:,3),1), ...
    spm_sample_vol(V3,input_vox(:,1),input_vox(:,2),input_vox(:,3),1));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function newName = rename(oldName , format, varargin)
% rename(fileTest, '%s/<name>--DRIM.png', dirDrimResult)

[path, name, ext] = fileparts(oldName);
if numel(varargin) ~= 0
    format = sprintf(format, varargin{:});
end
newName = strrep(format, '<path>', path);
newName = strrep(newName, '<name>', name);
newName = strrep(newName, '<ext>', ext(2:end)); % no . mark
end

author: baidut

repo: BIMEF

save_path: github-repos/MATLAB/baidut-BIMEF


text:
%FEATSELM Feature selection map
%
% [W,R] = FEATSELM(A,CRIT,METHOD,K,T,PAR1,...)
%
% INPUT
%   A      	Training dataset
%   CRIT   	Name of criterion: 'in-in', 'maha-s', 'NN' or others
%           (see FEATEVAL) or an untrained classifier V (default: 'NN')
%   METHOD  - 'forward' : selection by featself (default)
% 	        - 'float'   : selection by featselp
% 	        - 'backward': selection by featselb
% 	        - 'b&b'     : branch and bound selection by featselo
% 	        - 'ind'     : individual
% 	        - 'lr'      : plus-l-takeaway-r selection by featsellr
%           - 'sparse'  : use sparse untrained classifier CRIT
%   K      	Desired number of features (default: K = 0, return optimal set)
%   T      	Tuning set to be used in FEATEVAL (optional)
%   PAR1,.. Optional parameters:
% 	        - L,R       : for 'lr' (default: L = 1, R = 0)
%
% OUTPUT
%   W       Feature selection mapping
%   R       Matrix with step by step results
%
% DESCRIPTION
% Computation of a mapping W selecting K features. This routines offers a
% central interface to all other feature selection methods. W can be used
% for selecting features in a dataset B using B*W.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FEATEVAL, FEATSELO, FEATSELB, FEATSELI,
% FEATSELP, FEATSELF, FEATSELLR

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

function [w,res] = featselm(a,crit,arg3,ksel,t,par1,par2)

	if (nargin < 2 | isempty(crit))
		prwarning(2,'criterion not specified, assuming NN');
		crit = 'NN';
	end
	if (nargin < 3 | isempty(arg3))
		prwarning(2,'method not specified, assuming forward');
		arg3 = 'forward';
	end
	if (nargin < 4)
		ksel = [];
	end
	if (nargin < 5)
		prwarning(3,'no tuning set supplied (risk of overfit)');
		t = [];
	end
	if (nargin < 6), par1 = []; end;
	if (nargin < 7), par2 = []; end;

	% If no arguments are supplied, return an untrained mapping.

	if (nargin == 0) | (isempty(a))
		w = prmapping('featselm',{crit,arg3,ksel,t,par1,par2});
		w = setname(w,'Feature Selection');
		return
	end

	a = testdatasize(a);
	[m,k] = size(a);

	if (isstr(arg3))
		method = arg3;												% If the third argument is a string,
		switch (method)												%   it specifies the method to use.
		 case {'forward','featself'}
		  [w,res] = featself(a,crit,ksel,t);
		 case {'float','featselp'}
		  [w,res] = featselp(a,crit,ksel,t);
		 case {'backward','featselb'}
		  [w,res] = featselb(a,crit,ksel,t);
		 case {'b&b','featselo'}
		  [w,res] = featselo(a,crit,ksel,t);
		 case {'ind','featseli'}
		  [w,res] = featseli(a,crit,ksel,t);
		 case {'lr','featsellr'}
		  [w,res] = featsellr(a,crit,ksel,par1,par2,t);
		 case {'sparse'}
		  v = a*crit;
		  if isaffine(v)
		  	v = getdata(v,'rot');
		  	w = featsel(size(a,2),find(v(:,1) == 0));
			else
				v = getdata(v,'beta');
		  end
		  w = featsel(size(a,2),find(v(:,1) ~= 0));
		 otherwise
		  error('Unknown method specified.')
		end
	elseif (ismapping(arg3))
		w = arg3;												% If the third argument is a mapping,
		isuntrained(w);									%  assert it is untrained and train
		[w,res] = feval(mfilename,a,crit,w.mapping_file,ksel,t,par1,par2);
	else
		error('Illegal method specified.')
	end

return

author: marianux

repo: ecg-kit

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


text:
% SHELLGUI
%
% Files
%   shellgui - GUI interface for seashell function
%   seashell - draws a pretty Florida seashell, using a 3D parametric surface.
%

% Example
%   shellgui


text:
function [axis, axisLim] = extractAxisDataMultipleYAxes(obj,parentAxisData,yaxIndex)

    childAxisData = get(parentAxisData.YAxis(yaxIndex));

    %-------------------------------------------------------------------------%

    %-axis-side-%
    if yaxIndex == 1
        axis.side = 'left';
    elseif yaxIndex == 2
        axis.side = 'right';
    end

    %-------------------------------------------------------------------------%

    %-y-axis initializations-%
    axis.zeroline = false;
    axis.autorange = false;
    axis.exponentformat = obj.PlotlyDefaults.ExponentFormat;

    %-------------------------------------------------------------------------%

    %-y-axis tick font size-%
    axis.tickfont.size = childAxisData.FontSize;

    %-------------------------------------------------------------------------%

    %-y-axis tick font family-%
    axis.tickfont.family = matlab2plotlyfont(childAxisData.FontName);

    %-------------------------------------------------------------------------%

    ticklength = min(obj.PlotlyDefaults.MaxTickLength,...
        max(childAxisData.TickLength(1)*parentAxisData.Position(3)*obj.layout.width,...
        childAxisData.TickLength(1)*parentAxisData.Position(4)*obj.layout.height));
    %-y-axis ticklen-%
    axis.ticklen = ticklength;

    %-------------------------------------------------------------------------%

    %-y-axis coloring-%
    axiscol = sprintf('rgb(%f,%f,%f)', 255*childAxisData.Color);

    axis.linecolor = axiscol;
    axis.tickcolor = axiscol;
    axis.tickfont.color = axiscol;

    try
        axis.gridcolor = sprintf('rgba(%f,,%f,%f,%f)', 255*parentAxisData.GridColor, parentAxisData.GridAlpha);
    catch
        axis.gridcolor = axiscol;
    end

    %-------------------------------------------------------------------------%

    %-axis show grid-%
    if strcmp(parentAxisData.YGrid, 'on')
        axis.showgrid = true;
    else
        axis.showgrid = false;
    end

    %-------------------------------------------------------------------------%

    % grid = eval(['parentAxisData.' axisName 'Grid;']);
    % minorGrid = eval(['parentAxisData.' axisName 'MinorGrid;']);

    % if strcmp(grid, 'on') || strcmp(minorGrid, 'on')
    %     %-axis show grid-%
    %     axis.showgrid = true;
    % else
    %     axis.showgrid = false;
    % end

    %-------------------------------------------------------------------------%

    %-line widths-%
    linewidth = max(1,childAxisData.LineWidth*obj.PlotlyDefaults.AxisLineIncreaseFactor);

    axis.linewidth = linewidth;
    axis.tickwidth = linewidth;
    axis.gridwidth = linewidth;

    %-------------------------------------------------------------------------%

    %-axis type-%
    axis.type = childAxisData.Scale;

    %-------------------------------------------------------------------------%

    %-axis showtick labels / ticks-%
    tickValues = childAxisData.TickValues;

    if isempty(tickValues)

        %-axis ticks-%
        axis.ticks = '';
        axis.showticklabels = false;

        %-axis autorange-%
        axis.autorange = true;

        %---------------------------------------------------------------------%

        % switch parentAxisData.Box
        %     case 'on'
        %         %-axis mirror-%
        %         axis.mirror = true;
        %     case 'off'
        %         axis.mirror = false;
        % end

        %---------------------------------------------------------------------%

    else

        %-get axis limits-%
        axisLim = childAxisData.Limits;

        %-axis tick direction-%
        switch childAxisData.TickDirection
            case 'in'
                axis.ticks = 'inside';
            case 'out'
                axis.ticks = 'outside';
        end

        %---------------------------------------------------------------------%
        % switch parentAxisData.Box
        %     case 'on'
        %         %-axis mirror-%
        %         axis.mirror = 'ticks';
        %     case 'off'
        %         axis.mirror = false;
        % end

        %---------------------------------------------------------------------%

        %-LOG TYPE-%
        if strcmp(axis.type, 'log')

            %-axis range-%
            axis.range = log10(axisLim);
            %-axis autotick-%
            axis.autotick = true;
            %-axis nticks-%
            axis.nticks = length(tickValues) + 1;

        %---------------------------------------------------------------------%

        %-LINEAR TYPE-%
        elseif strcmp(axis.type, 'linear')

            %-----------------------------------------------------------------%

            %-get tick label mode-%
            tickLabelMode = childAxisData.TickLabelsMode;

            %-----------------------------------------------------------------%

            %-AUTO MODE-%
            if strcmp(tickLabelMode, 'auto')

                %-------------------------------------------------------------%

                if isnumeric(axisLim)
                    axis.range = axisLim;

                %-------------------------------------------------------------%

                elseif isduration(axisLim)
                   [temp,type] = convertDuration(axisLim);

                   if (~isduration(temp))
                       axis.range = temp;
                       axis.type = 'duration';
                       axis.title = type;
                   else
                       nticks = length(tickValues) + 1;
                       delta = 0.1;
                       axis.range = [-delta nticks+delta];
                       axis.type = 'duration - specified format';
                   end

                %-------------------------------------------------------------%

                elseif isdatetime(axisLim)
                    axis.range = convertDate(axisLim);
                    axis.type = 'date';
                else
                    % data is a category type other then duration and datetime
                end

                %-axis autotick-%
                axis.autotick = true;
                %-axis numticks-%
                axis.nticks = length(tickValues) + 1;
                axis.showticklabels = true;

            %-----------------------------------------------------------------%

            %-CUSTOM MODE-%
            else

                %-------------------------------------------------------------%

                %-get tick labels-%
                tickLabels = childAxisData.TickLabels;

                %-------------------------------------------------------------%

                %-hide tick labels as lichkLabels field is empty-%
                if isempty(tickLabels)

                    %-------------------------------------------------------------%

                    %-hide tick labels-%
                    axis.showticklabels = false;

                    %-------------------------------------------------------------%

                    %-axis autorange-%
                    axis.autorange = true;

                    %-------------------------------------------------------------%

                %-axis show tick labels as tickLabels matlab field-%
                else

                    %-------------------------------------------------------------%

                    axis.showticklabels = true;
                    % axis.type = 'linear';

                    %-------------------------------------------------------------%

                    if isnumeric(axisLim)
                        axis.range = axisLim;
                    else
                        axis.autorange = true;
                    end

                    %-------------------------------------------------------------%

                    axis.tickvals = tickValues;
                    axis.ticktext = tickLabels;

                    %-------------------------------------------------------------%

                    %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%
                    % NOTE:
                    %    The next piece of code was replaced by the previous one.
                    %    I think that the new piece of code is better, optimal and
                    %    extends to all cases. However, I will leave this piece of
                    %    code commented in case there is a problem in the future.
                    %
                    %    If there is a problem with the new piece of code, please
                    %    comment and uncomment the next piece of code.
                    %
                    %    If everything goes well with the new gripping piece, at
                    %    the end of the development we will be able to remove the
                    %    commented lines
                    %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%

                    %-axis labels
                    % labels = str2double(tickLabels);
                    % try
                    %     %find numbers in labels
                    %     labelnums = find(~isnan(labels));
                    %     %-axis type linear-%
                    %     axis.type = 'linear';
                    %     %-range (overwrite)-%
                    %     delta = (labels(labelnums(2)) - labels(labelnums(1)))/(labelnums(2)-labelnums(1));
                    %     axis.range = [labels(labelnums(1))-delta*(labelnums(1)-1) labels(labelnums(1)) + (length(labels)-labelnums(1))*delta];
                    %     %-axis autotick-%
                    %     axis.autotick = true;
                    %     %-axis numticks-%
                    %     axis.nticks = eval(['length(parentAxisData.' axisName 'Tick) + 1;']);
                    % catch
                    %     %-axis type category-%
                    %     axis.type = 'category';
                    %     %-range (overwrite)-%
                    %     axis.autorange = true;
                    %     %-axis autotick-%
                    %     % axis.autotick = true;
                    % end
                end
            end
        end
    end

    %-------------------------------------------------------------------------%

    %-scale direction-%
    if strcmp(childAxisData.Direction, 'reverse')
        axis.range = [axis.range(2) axis.range(1)];
    end

    %-------------------------------------------------------------------------%

    %-y-axis label-%
    label = childAxisData.Label;
    labelData = get(label);

    %STANDARDIZE UNITS
    fontunits = get(label,'FontUnits');
    set(label,'FontUnits','points');

    %-------------------------------------------------------------------------%

    %-title settings-%
    if ~isempty(labelData.String)
        axis.title = parseString(labelData.String,labelData.Interpreter);
    end

    axis.titlefont.color = sprintf('rgb(%f,%f,%f)', 255*labelData.Color);
    axis.titlefont.size = labelData.FontSize;
    axis.titlefont.family = matlab2plotlyfont(labelData.FontName);

    %-------------------------------------------------------------------------%

    %REVERT UNITS
    set(label,'FontUnits',fontunits);

    %-------------------------------------------------------------------------%

    if strcmp(childAxisData.Visible,'on')
        %-axis showline-%
        axis.showline = true;
    else
        %-axis showline-%
        axis.showline = false;
        %-axis showticklabels-%
        axis.showticklabels = false;
        %-axis ticks-%
        axis.ticks = '';
        %-axis showline-%
        axis.showline = false;
        %-axis showticklabels-%
        axis.showticklabels = false;
        %-axis ticks-%
        axis.ticks = '';
    end

    %-------------------------------------------------------------------------%
end

author: plotly


text:
function [T, Y]=FOChuaNR(parameters, orders, TSim, Y0)
%
% Numerical Solution of the Fractional-Order Chua's System
% with piecewise linear nonlinearity defined in function f_x()
%
%   D^q1 x(t) = alpha(y(t) - x(t) - f(x))
%   D^q2 y(t) = x(t) - y(t) + z(t)
%   D^q3 z(t) = -beta y(t) - gamma z(t)
%   where  f(x) is defined as f_x(m0, m1, x);
%
% function [T, Y] = FOChuaNR(parameters, orders, TSim, Y0)
%
% Input:    parameters - model parameters [alpha, beta, gamma, m0, m1]
%           orders - derivatives orders [q1, q2, q3]
%           TSim - simulation time (0 - TSim) in sec
%           Y0 - initial conditions [Y0(1), Y0(2), Y0(3)]
%
% Output:   T - simulation time (0 : Tstep : TSim)
%           Y - solution of the system (x=Y(1), y=Y(2), z=Y(3))
%
% Author:  (c) Ivo Petras (ivo.petras@tuke.sk), 2010.
%

% time step:
h=0.005;
% number of calculated mesh points:
n=round(TSim/h);
%orders of derivatives, respectively:
q1=orders(1); q2=orders(2); q3=orders(3);
% constants of Volta's system:
alpha=parameters(1); beta=parameters(2); gamma=parameters(3);
m0=parameters(4); m1=parameters(5);
% binomial coefficients calculation:
cp1=1; cp2=1; cp3=1;
for j=1:n
    c1(j)=(1-(1+q1)/j)*cp1;
    c2(j)=(1-(1+q2)/j)*cp2;
    c3(j)=(1-(1+q3)/j)*cp3;
    cp1=c1(j); cp2=c2(j); cp3=c3(j);
end
% initial conditions setting:
x(1)=Y0(1); y(1)=Y0(2); z(1)=Y0(3);
% calculation of phase portraits /numerical solution/:
for i=2:n
    x(i)=(alpha*(y(i-1)-x(i-1)-f_x(m1,m0,x(i-1))))*h^q1 - memo(x, c1, i);
    y(i)=(x(i)-y(i-1)+z(i-1))*h^q2 - memo(y, c2, i);
    z(i)=(-beta*y(i)-gamma*z(i-1))*h^q3 - memo(z, c3, i);
end
for j=1:n
    Y(j,1)=x(j);
    Y(j,2)=y(j);
    Y(j,3)=z(j);
end
T=h:h:TSim;
%

author: Sable

repo: mcbench-benchmarks


text:
%% Import Script for Tensor Data
%
% This script was automatically created by the import wizard. You should
% run the whoole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.

%% Specify Crystal and Specimen Symmetries

% crystal symmetry
CS = {crystal symmetry};

% plotting convention
setMTEXpref('xAxisDirection',{xAxisDirection});
setMTEXpref('zAxisDirection',{zAxisDirection});

%% Specify File Names

% path to files
pname = {path to files};

% which files to be imported
fname = {file names};

%% Import the Data

% create a Tensor variable containing the data
tensor = tensor.load(fname,CS,'interface',{interface},{options});

author: mtex-toolbox

repo: mtex

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


text:
classdef Max < dagnn.ElementWise
  %MAX DagNN max layer
  %   The Max layer takes the max of all its inputs and store the result
  %   as its only output.

  properties (Transient)
    numInputs
%     inputSize
  end

  methods
    function outputs = forward(obj, inputs, params)
      obj.numInputs = numel(inputs) ;
%       [h w d]=size(inputs{1});
%       obj.inputSize = ones(h,w,d);
      outputs{1} = inputs{1} ;
      for k = 2:obj.numInputs
          outputs{1} = max(outputs{1},inputs{k});

%           for hight = 1:h
%               for width = 1:w
%                   for depth = 1:d
%         if(outputs{1}(hight,width,depth) < inputs{k}(hight,width,depth))
%             outputs{1}(hight,width,depth) = inputs{k}(hight,width,depth);
%             obj.inputSize(hight,width,depth) = k;
%         end
%                   end
%               end
%           end
      end
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
        sz = size(inputs{1});
        if numel(sz) < 4, sz(4) = 1;
        sz(5) = obj.numInputs;
        tempinputs = single(zeros(sz));
        tempinputs = gpuArray(tempinputs);
        end
      for k = 1:obj.numInputs

          tempinputs(:,:,:,:,k) = inputs{k};
      end
         testinputs{1} = permute(tempinputs, [1 2 5 3 4]);
         sz_in = size(testinputs{1});
         [~, idx]  = mex_maxpool3d(testinputs{1},...
        'pool',[1,1,obj.numInputs], 'stride', [1,1,obj.numInputs], 'pad', [0,0,0,0,0,0]);

%          sz = size(derOutputs{1});
%          if numel(sz) < 4, sz(4) = 1;
%          sz(5) = obj.numInputs;
%          tempoutputs = single(zeros(sz));
% %          tempoutputs = gpuArray(tempoutputs);
%          end
%          for k = 1:obj.numInputs
%              tempoutputs(:,:,:,:,k) = derOutputs{1};
%          end

         derOutputs{1} = permute(derOutputs{1}, [1 2 5 3 4]);
         derInputs{1} = mex_maxpool3d(derOutputs{1}, idx, sz_in ,...
        'pool',[1,1,obj.numInputs], 'stride', [1,1,obj.numInputs], 'pad', [0,0,0,0,0,0]);
         tempoutputs = permute(derInputs{1}, [1 2 4 5 3]);
         for k = 1:obj.numInputs

          derInputs{k} = tempoutputs(:,:,:,:,k);
         end
%         for hight = 1:h
%               for width = 1:w
%                   for depth = 1:d
%                derInputs{obj.inputSize(hight,width,depth)}(hight,width,depth) = derOutputs{1}(hight,width,depth);
% %              if (obj.inputSize(hight,width,depth) ==k)
% %           derInputs{k}(hight,width,depth) = derOutputs{1}(hight,width,depth) ; %not sure whether times obj.numInputs
% %              end
%                   end
%               end
%         end

      derParams = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes{1} = inputSizes{1} ;
      for k = 2:numel(inputSizes)
        if all(~isnan(inputSizes{k})) && all(~isnan(outputSizes{1}))
          if ~isequal(inputSizes{k}, outputSizes{1})
            warning('Max layer: the dimensions of the input variables is not the same.') ;
          end
        end
      end
    end

    function rfs = getReceptiveFields(obj)
      numInputs = numel(obj.net.layers(obj.layerIndex).inputs) ;
      rfs.size = [1 1] ;
      rfs.stride = [1 1] ;
      rfs.offset = [1 1] ;
      rfs = repmat(rfs, numInputs, 1) ;
    end

    function obj = Max(varargin)
      obj.load(varargin) ;
    end
  end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function Tmodel = cleanTmodel(Tmodel, varargin)
% Reorganizes and checks `Tmodel` for completeness. To be used
% immediately after `addToTmodel`. Called by `addToTmodel`, `mergeModelsBorg`, `mergeTrxns`, calls `TmodelFields`, `buildRxnEquations`, `removeDuplicateNames`, `makeNamesUnique`, `charpos`, `removeProblematicCharacters`.
%
% USAGE:
%
%    Tmodel = cleanTmodel(Tmodel)
%
% INPUTS:
%    Tmodel:       Template model.
%
% OPTIONAL INPUTS:
%    'Verbose':    Ask when performing certain operations.
%
% OUTPUTS:
%    Tmodel:       Cleaned template model
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

verbose = false ; % Declare variables.
if ~isempty(varargin)
    if sum(strcmp('Verbose', varargin))
        verbose = true ;
    end
end

% Grab fields names.
fields = TmodelFields ;

% Reaction related cell field names.
rxnFields = fields{1} ;
% Reaction related double Fields (which are kept model specific).
rNumFields = fields{2} ;
% Metabolite related cell field names.
metFields = fields{3} ;
% Metabolite related double field names.
mNumFields = fields{4} ;

% Model names.
moNames = fieldnames(Tmodel.Models) ;

% Number of reactions and number of metabolites.
nRxns = length(Tmodel.rxns) ;
nMets = length(Tmodel.mets) ;

% Set flags for if Verbose mode is on
if isempty(varargin)
    mergeRxns = 'y' ;
    askeverytime = 0 ;
    beautifynames = 'y' ;
    askeverytimeb = 0 ;
else
    if verbose
        mergeRxns = input(['Should I try merging identical reactions? ' char(10) ...
                               'y...yes in all cases' char(10) ...
                               'a...ask me every time' char(10) ...
                               'n...no thanks (default)' char(10)], 's') ;
        askeverytime = strcmpi(mergeRxns, 'a') ;

        beautifynames = input(['Should I try to replace cryptic names by nicer ones? ' char(10) ...
                           'y...yes in all cases' char(10) ...
                           'a...ask me every time' char(10) ...
                           'n...no thanks (default)' char(10)], 's') ;
        askeverytimeb = strcmpi(beautifynames,'a') ;
    else
        mergeRxns = 'y' ;
        askeverytime = 0 ;
        beautifynames = 'y' ;
        askeverytimeb = 0 ;
    end
end

%% Lengthen rev, lb, ub, and c arrays for old models.
for iField = 1:length(rNumFields)
    for iMo = 1:length(moNames)
        nowArray = Tmodel.(rNumFields{iField}).(moNames{iMo}) ;
        if length(nowArray) ~= nRxns
            nowArray(end:nRxns) = 0 ;
            Tmodel.(rNumFields{iField}).(moNames{iMo}) = nowArray;
        end
    end
end

%% Lengthen reaction and metabolite identity arrays, convert to logicals.
for iMo = 1:length(moNames)
    % Reactions.
    rxnsArray = Tmodel.Models.(moNames{iMo}).rxns ;
    if length(rxnsArray) ~= nRxns
        rxnsArray(end+1:nRxns) = 0 ;
    end
    Tmodel.Models.(moNames{iMo}).rxns = logical(rxnsArray) ;

    % Metabolites.
    metsArray = Tmodel.Models.(moNames{iMo}).mets ;
    if length(metsArray) ~= nMets
        metsArray(end+1:nMets) = 0 ;
    end
    Tmodel.Models.(moNames{iMo}).mets = logical(metsArray) ;
end

%% Check for accidently duplicated reactions and metabolites
Sfind = zeros(size(Tmodel.S)) ;
Sfind(Tmodel.S ~= 0) = 1 ;
metpos = cell(1, nRxns) ;
metposstr = cell(1, nRxns) ;
for i = 1:nRxns
    metpos{i} = find(Sfind(:, i)) ;
    metposstr{i} = num2str(metpos{i}') ;
end
[~, a, b] = unique(metposstr) ;
mergerxns = {} ; mergerxnratio = [] ;
modelnames = fieldnames(Tmodel.Models)' ;
for im = 1:length(modelnames)
    rxnNowPresent(:, im) = Tmodel.Models.(modelnames{im}).rxns ;
end

if ismember(mergeRxns, {'y', 'a', 'Y', 'A'})
    for i = 1:nRxns
        if ~ismember(i, a) % rxns that use the same metabolites
            eqset = find(b == b(i)) ;
            stoicheck = 1 ;
            for ie = 2:length(eqset)
                rxnratio = unique(Tmodel.S(metpos{i}, eqset(1)) ./ Tmodel.S(metpos{i}, eqset(ie))) ;
                if length(rxnratio) == 1 && ~isnan(rxnratio)  && ...
                        sum(rxnNowPresent(eqset(1), :) .* rxnNowPresent(eqset(ie), :)) == 0
                    if askeverytime
                    % if reactions have the same stoichiometry, maybe with a factor
                    disp([Tmodel.rxns{eqset(1)} '(' strjoin(modelnames(rxnNowPresent(eqset(1), :)), '|') ')' ' and ' ...
                          Tmodel.rxns{eqset(ie)}  '(' strjoin(modelnames(rxnNowPresent(eqset(ie), :)), '|') ')'  ...
                          ' seem to have the same stoichiometry' char(10) ...
                         Tmodel.rxnEquations{eqset(1)} '  VS  ' Tmodel.rxnEquations{eqset(ie)} char(10) ])
                    shouldmerge = input('Do you want to merge them? [Y/n] ' ,'s') ;
                    elseif strcmpi(mergeRxns, 'y')
                        shouldmerge = 'y' ;
                    end
                    if ~strcmp(shouldmerge, 'n')
                        mergerxns{end + 1} = eqset ;
                        mergerxnratio(end + 1) = rxnratio ;
                        disp(['merging ' Tmodel.rxns{eqset(1)} '(' strjoin( modelnames(rxnNowPresent(eqset(1), :)), '|') ')' ' and ' ...
                             Tmodel.rxns{eqset(ie)}  '(' strjoin( modelnames(rxnNowPresent(eqset(ie), :)), '|') ')' ])
                    end
                end
            end
        end
    end
    if length(mergerxns) ~= 0
        Tmodel = mergeTrxns(Tmodel, mergerxns, mergerxnratio) ;
        disp(char(10))
    end
end

% Re-count number of reactions and number of metabolites.
nRxns = length(Tmodel.rxns) ;
nMets = length(Tmodel.mets) ;

%% Prefer biologist-friendly names for reactions and metabolites
if ismember(beautifynames,{'y', 'a', 'Y', 'A'})
    for im = 1:nMets
        nowIDsep = {Tmodel.mets{im}(1:end-3)} ;
        nowcomp = Tmodel.mets{im}(end-2:end) ;
        % parse alternative metIDs
        nowIDs = Tmodel.metID{im};
        nowdelim1 = strfind(nowIDs,':') ;
        nowdelim2 = [strfind(nowIDs,'|') strfind(nowIDs,':') length(nowIDs) + 1] ;
        for ii = 1:length(nowdelim1)
            nowIDsep{ii + 1} = nowIDs(nowdelim1(ii) + 1:min(nowdelim2(nowdelim2 > nowdelim1(ii))) - 1) ;
            nowIDsep{ii + 1} = lower(nowIDsep{ii + 1}(1:min([strfind(nowIDsep{ii + 1},'[') - 1, length(nowIDsep{ii + 1})]))) ;
            nowIDsep{ii + 1} = removeProblematicCharacters(nowIDsep{ii + 1}) ;
        end
        [nowIDalt, ~, freq] = unique(nowIDsep) ;
        nowIDalt = nowIDalt(~cellfun(@isempty, nowIDalt)) ;
        tailnum = [] ; idfreq = [] ;
        if length(nowIDalt) == 1
            % no alternative name found, keep the current one
        else
            for ii = 1:length(nowIDalt)
                tailnum(ii) = length(nowIDalt{ii}(max(find(charpos(nowIDalt{ii}))):end)) - 1 ;
                idfreq(ii) = sum(freq == ii) ;
            end
            [~, goodids] = min((tailnum + 2)./(2*idfreq)) ;
            % if leading character is a number in original name but "a" in new
            % name, then do nothing
            nowcharpos = charpos(Tmodel.mets{im}) ;
            if ~strcmp(nowIDalt{goodids}, Tmodel.mets{im}(1:end - 3)) && ...
                    ~(~nowcharpos(1) && strcmp(nowIDalt{goodids}(1), 'a')) && ...
                   length(Tmodel.mets{im}(max(find(charpos(Tmodel.mets{im}(1:end - 3)) ) ):end - 3)) > 4
                if askeverytimeb
                    replacethistime = input(['Should I rename ' Tmodel.mets{im} ' to '...
                                             nowIDalt{goodids} nowcomp ' (y/n)?'], 's') ;
                end
                if ~askeverytimeb || strcmpi(replacethistime, 'y')
                    disp(['renaming ' Tmodel.mets{im} ' to ' nowIDalt{goodids} nowcomp])
                    Tmodel.mets{im} = [nowIDalt{goodids} nowcomp] ;
                end
            end
        end
    end
    for ir = 1:nRxns
        % transport to extracellular reaction
        if numel(find(Tmodel.S(:, ir))) == 1
            if strcmp(Tmodel.mets{(Tmodel.S(:, ir))~=0}(end - 2:end), '[e]')
                newrxnname = ['ex_' lower(Tmodel.mets{(Tmodel.S(:, ir))~=0}(1:end - 3)) '_e'] ;
                if ~strcmp(newrxnname, Tmodel.rxns{ir})
                    if askeverytimeb
                        replacethistime = input(['Should I rename '  Tmodel.rxns{ir}  ' to '...
                                             newrxnname ' (y/n)?'], 's') ;
                    end
                    if ~askeverytimeb || strcmpi(replacethistime, 'y')
                        disp(['correcting ' Tmodel.rxns{ir} ' to ' newrxnname])
                        Tmodel.rxns{ir} = newrxnname ;
                    end
                end
            end
        else % normal reactions
            nowIDsep = Tmodel.rxns(ir) ;
            if strcmp(nowIDsep{1}(end - 1),'_')
                % this reaction has probably been already renamed to ensure
                % unique reaction names. Therefore do nothing.
                continue
            end
            % parse alternative rxnIDs
            nowIDs = Tmodel.rxnID{ir} ;
            nowdelim1 = strfind(nowIDs, ':') ;
            nowdelim2 = [strfind(nowIDs, '|') strfind(nowIDs, ':') length(nowIDs) + 1] ;
            for ii = 1:length(nowdelim1)
                nowIDsep{ii+1} = nowIDs(nowdelim1(ii) + 1:min(nowdelim2(nowdelim2 > nowdelim1(ii))) - 1) ;
                nowIDsep{ii+1} = lower(nowIDsep{ii + 1}(1:min([strfind(nowIDsep{ii + 1},'[') - 1, length(nowIDsep{ii + 1})]))) ;
                nowIDsep{ii+1} = removeProblematicCharacters(nowIDsep{ii + 1}) ;
            end
            [nowIDalt, ~, freq] = unique(nowIDsep) ;
            tailnum = [] ; idfreq = [] ;
            if length(nowIDalt) == 1
                % no alternative name found, keep the current one
            else
                for ii = 1:length(nowIDalt)
                    tailnum(ii) = length(nowIDalt{ii}(max(find(charpos(nowIDalt{ii}))):end)) - 1 ;
                    idfreq(ii) = sum(freq == ii) ;
                end
                [~, goodids] = min((tailnum + 2) ./ (2 * idfreq)) ;
                % check that the name is not already taken
                if sum(strcmp(Tmodel.rxns, nowIDalt{goodids})) == 0 && ...
                   length(Tmodel.rxns{ir}(max(find(charpos(Tmodel.rxns{ir}))):end)) > 4
                    if askeverytimeb
                        replacethistime = input(['Should I rename '  Tmodel.rxns{ir}  ' to '...
                                             nowIDalt{goodids} ' (y/n)?'], 's') ;
                    end
                    if ~askeverytimeb || strcmpi(replacethistime, 'y')
                        disp(['replacing ' Tmodel.rxns{ir} ' by ' nowIDalt{goodids}])
                        Tmodel.rxns{ir} = nowIDalt{goodids} ;
                    end
                end
            end
        end
    end
end

%% Make sure all names are unique before reorganizing.
% Rebuild reaction equations so there is a hint for renaming rxns.
Tmodel = buildRxnEquations(Tmodel) ;
if verbose ; fprintf('Checking if reaction IDs (.rxns) are unique.\n') ; end
Tmodel.rxns = makeNamesUnique(Tmodel.rxns, Tmodel.rxnEquations) ;
if verbose ; fprintf('Checking if metabolite IDs (.mets) are unique.\n') ; end
Tmodel.mets = makeNamesUnique(Tmodel.mets, Tmodel.metNames) ;

% Remove duplicate met IDs.
Tmodel.metID = removeDuplicateNames(Tmodel.metID) ;

% Rebuild reaction equations.
Tmodel = buildRxnEquations(Tmodel) ;

%% Reorder reactions and metabolites alphabetically.
% Reorder reaction related lists.
[~, rxnI] = sort(Tmodel.rxns) ;
for iField = 1:length(rxnFields)
    Tmodel.(rxnFields{iField}) = Tmodel.(rxnFields{iField})(rxnI) ;
end
for iField = 1:length(rNumFields)
    for iMo = 1:length(moNames)
        Tmodel.(rNumFields{iField}).(moNames{iMo}) = ...
            Tmodel.(rNumFields{iField}).(moNames{iMo})(rxnI) ;
    end
end

% Reorder metabolite related lists.
[~, metI] = sort(Tmodel.mets) ;
for iField = 1:length(metFields)
    Tmodel.(metFields{iField}) = Tmodel.(metFields{iField})(metI) ;
end
for iField = 1:length(mNumFields)
    Tmodel.(mNumFields{iField}) = ...
        Tmodel.(mNumFields{iField})(metI) ;
end

% Reorder identity arrays.
for iMo = 1:length(moNames)
    Tmodel.Models.(moNames{iMo}).rxns = ...
        Tmodel.Models.(moNames{iMo}).rxns(rxnI) ;
    Tmodel.Models.(moNames{iMo}).mets = ...
        Tmodel.Models.(moNames{iMo}).mets(metI) ;
end

% Reorder S matrix.
Tmodel.S = Tmodel.S(:, rxnI) ;
Tmodel.S = Tmodel.S(metI, :) ;

%% Share information between the same metabolite in different compartments.
% Metabolite fields to share information between
share = {'metNames' 'metFormulas' 'metKEGGID' 'metSEEDID' 'metChEBIID' ...
         'metPubChemID' 'metInChIString'} ;

metsNoComp = cell(nMets,1) ;
for iMet = 1:nMets
    metsNoComp{iMet} = Tmodel.mets{iMet}(1:end - 3) ;
end
uniqMets = unique(metsNoComp,'first') ;

% Create combined formation fields for each metabolite and redistribute.
for iMet = 1:length(uniqMets)
  % Find sister metabolites
  sMets = find(strcmp(uniqMets{iMet}, metsNoComp)) ;
  if length(sMets) > 1
    % Create combined information for each field.
    for iF = 1:length(share)
      info = {''} ;
      for iSis = 1:length(sMets)
        % If the field is not empty, break information into parts
        % (as seperated by pipes), and add that information in if
        % it is not already present.
        if ~isempty(Tmodel.(share{iF}){sMets(iSis)})
          pipePos = [0 ...
                     strfind(Tmodel.(share{iF}){sMets(iSis)}, '|') ...
                     length(Tmodel.(share{iF}){sMets(iSis)}) + 1] ;
          for j = 1:length(pipePos) - 1
            nowInfo = Tmodel.(share{iF}){sMets(iSis)} ...
                      (pipePos(j) + 1:pipePos(j + 1) - 1) ;
            % If the info does not already exist.
            if isempty(strfind(info{1}, nowInfo))
              % If it is the first piece of information or additional.
              if isempty(info{1})
                info{1} = nowInfo ;
              else
                info{1} = strcat(info{1}, '|', nowInfo) ;
              end
            end
          end
        end
      end
      for iSis = 1:length(sMets)
          Tmodel.(share{iF}){sMets(iSis)} = info{1} ;
      end
      clear info
    end
  end
end

% Now remove duplicate information which may have snuck in.
for iField = 1:length(share)
    Tmodel.(share{iField}) = removeDuplicateNames(Tmodel.(share{iField})) ;
end

%% Remove duplicate reaction longnames and IDs
Tmodel.rxnNames = removeDuplicateNames(Tmodel.rxnNames) ;
Tmodel.rxnID = removeDuplicateNames(Tmodel.rxnID) ;

%% Clear leftover fields from the matching process
deleteFields = {'rxnNamesFix' 'rxnComp' 'metNums' 'rxnMetNames'} ;
existFields = isfield(Tmodel, deleteFields) ;
for iF = 1:length(existFields)
    if existFields(iF)
        Tmodel = rmfield(Tmodel,deleteFields{iF}) ;
    end
end

author: opencobra

repo: cobratoolbox


text:
function [tt1]=conj(tt)
%Compute a complex conjugate of TT-tensor
%   [TT1]=CONJ(TT)
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

tt1=tt;
tt1.core=conj(tt.core);
return
end

author: oseledets

repo: TT-Toolbox


text:
function done = mission2(occupiers, player)

captured = zeros(1,6);
for i = 1:6
    captured(i) = all(occupiers{i} == player);
end

done = captured(3) && captured(5) && sum(captured) > 2;


text:
function hmObj = buildHeadModel(eeg,surfFile,templateFile,plotModel,scalingFactor)
% surfFile is the path to the output surface file that will be generated
% templateFile is the path to an existing template file

% electrode locations [X Y Z]
elocs   = [[eeg.chanlocs.X]' [eeg.chanlocs.Y]' [eeg.chanlocs.Z]']*scalingFactor;
labels  = {eeg.chanlocs.labels};

% remove any electrodes which lack locations
unknownLocs = find(arrayfun(@(x) isempty(x.X), eeg.chanlocs));
if ~isempty(unknownLocs)
    lbl = cellfun(@(x) [x ', '],labels(unknownLocs),'UniformOutput',false);
    lbl{end} = lbl{end}(1:end-2);
    fprintf('Electrodes [%s] do not have locations. Removing...\n',char(lbl{:})');
    labels(unknownLocs) = [];
end

%
% % build the head model object
% hmObj = headModel('surfaces',surfFile,'atlas',template.atlas,'fiducials',template.fiducials,'channelSpace',elocs,'label',labels);
% if plotModel
%     plotHeadModel(hmObj);
% end

hmObj = headModel('channelSpace',elocs,'label',labels);
hmObj.warpTemplate2channelSpace(templateFile,surfFile);
% hmObj.warpChannelSpace2Template(templateFile,surfFile);

if plotModel
    plotHeadModel(hmObj);
end

%% solving the forward problem with OpenMEEG
conductivity = [0.33 0.022 0.33]; % brain and scalp = 0.33 S/m, skull = 0.022 S/m; these conductivies were taken from
                                  % Valdes-Hernandez et al., 2009, Oostendrop TF, 2000; Wendel and Malmivuo, 2006
normal2surface = true;
hmObj.computeLeadFieldBEM(conductivity,normal2surface);

repo: REST


text:
function out = convpower( in, uin, uout )
%  CONVPOWER Convert from power units to desired power units.
%   OUT = CONVPOWER( IN, UI, UO ) converts the input power IN (floating
%   point array) from unit specified in UI (string) to unit specified in UO
%   (string). If UO is not specified, IN is converted to the SI unit, which
%   is Watts.
%
%   Allowable input/output strings (not case sensitive):
%      'hp'          :horsepower
%      'W'           :watt
%      'kW'          :kilowatt
%      'lbf-ft/s'    :pound-feet/second (also 'lb-ft/s' or 'ft-lb/s')
%
%   Example:
%
%   Convert a matrix of power values from lbf-ft/s to hp:
%       hp = convpower([2 3; 4 5], 'lbf-ft/s', 'hp')
%
%   See also in aerospace toolbox CONVACC, CONVANG, CONVANGACC, CONVANGVEL,
%   CONVDENSITY, CONVFORCE, CONVLENGTH, CONVMASS, CONVPRES, CONVTEMP,
%   CONVVEL.

%   Author: Sky Sartorius
%   http://www.mathworks.com/matlabcentral/fileexchange/authors/101715

if ~isfloat( in )
    error('Input is not floating point');
end
if nargin < 3
    uout = 'w';
end
uin = lower(uin);
uout = lower(uout);
%conversion to watts
if strcmp('hp',uin)
    slope = 745.6998716; %W/hp
elseif strcmp('kw',uin)
    slope = 1000; %W/kW
elseif strcmp('lbf-ft/s',uin) || strcmp('lb-ft/s',uin) || strcmp('ft-lb/s',uout)
    slope = 1.355817948363637; %lb-ft/s/W
elseif strcmp('w',uin)
    slope = 1;
else
    error('invalid input unit string')
end

%conversion from watts
if strcmp('hp',uout)
    slope = slope/745.6998716; %W/hp
elseif strcmp('kw',uout)
    slope = slope/1000; %W/kW
elseif strcmp('lbf-ft/s',uout) || strcmp('lb-ft/s',uout) || strcmp('ft-lb/s',uout)
    slope = slope/1.355817948363637; %lb-ft/s/W
elseif strcmp('w',uout)
%     slope = slope/1;
else
    error('invalid output unit string')
end

out = in.*slope;


text:
function P=caesar_D(C,k)
% P=plain text
% k=key (can be 0 to 26); k=0 mean No encryption
% Example: P='Hello World' k=3
% C=caesar(P,k), C=Khoor Zruog
%If you have any problem or feedback please contact me @
%%===============================================
% NIKESH BAJAJ
% Asst. Prof., Lovely Professional University, India
% Almameter: Aligarh Muslim University, India
% +919915522564, bajaj.nikkey@gmail.com
%%===============================================

P=double(C)-k;
l=find(P<65);
P(l)=P(l)+26;

l=find(P<97);
l=find(P(l)>90)
P(l)=P(l)+26;

l=find(C==32);
P(l)=32;

P=char(P);
disp(' ')
disp('Cipher Text C =')
disp(C)
disp(' ')
disp('Plain Text')


text:
function [ value, ifault ] = alogam ( x )

%*****************************************************************************80
%
%% ALOGAM computes the logarithm of the Gamma function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by Malcolm Pike, David Hill.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Malcolm Pike, David Hill,
%    Algorithm 291:
%    Logarithm of Gamma Function,
%    Communications of the ACM,
%    Volume 9, Number 9, September 1966, page 684.
%
%  Parameters:
%
%    Input, real X, the argument of the Gamma function.
%    X should be greater than 0.
%
%    Output, real ALOGAM, the logarithm of the Gamma
%    function of X.
%
%    Output, integer IFAULT, error flag.
%    0, no error.
%    1, X <= 0.
%
  if ( x <= 0.0 )
    ifault = 1;
    value = 0.0;
    return
  end

  ifault = 0;
  y = x;

  if ( x < 7.0 )

    f = 1.0;
    z = y;

    while ( z < 7.0 )
      f = f * z;
      z = z + 1.0;
    end

    y = z;
    f = - log ( f );

  else

    f = 0.0;

  end

  z = 1.0 / y / y;

  value = f + ( y - 0.5 ) * log ( y ) - y ...
    + 0.918938533204673 + ...
    ((( ...
    - 0.000595238095238   * z ...
    + 0.000793650793651 ) * z ...
    - 0.002777777777778 ) * z ...
    + 0.083333333333333 ) / y;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [MHz] = radps2MHz(radps)
% Convert frequency from radians per second to megahertz.
% Chad A. Greene 2012
MHz = radps/(2*pi)/1000000;


text:
function x = pad2minsize(x, minsize, v)
% PAD2MINSIZE pads INPUT with given VALUE so that both spatial sizes are at least MINSIZE.
%
%   X = pad2minsize(X, MINSIZE, V)
%     Default minimum image size MINSIZE is 70.
%     Default padding value V is 0.

	if ~exist('v'), v = 0; end
	if ~exist('minsize'), minsize = 70; end
	x = padarray(x, [max(minsize-size(x,1), 0), max(minsize-size(x,2), 0)], v, 'post');

author: filipradenovic

repo: cnnimageretrieval


text:
function centroids = computeCentroids(X, idx, K)
%COMPUTECENTROIDS returs the new centroids by computing the means of the
%data points assigned to each centroid.
%   centroids = COMPUTECENTROIDS(X, idx, K) returns the new centroids by
%   computing the means of the data points assigned to each centroid. It is
%   given a dataset X where each row is a single data point, a vector
%   idx of centroid assignments (i.e. each entry in range [1..K]) for each
%   example, and K, the number of centroids. You should return a matrix
%   centroids, where each row of centroids is the mean of the data points
%   assigned to it.
%

% Useful variables
[m n] = size(X);

% You need to return the following variables correctly.
centroids = zeros(K, n);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every centroid and compute mean of all points that
%               belong to it. Concretely, the row vector centroids(i, :)
%               should contain the mean of the data points assigned to
%               centroid i.
%
% Note: You can use a for-loop over the centroids to compute this.
%

for i = 1:K
    idx_i = idx == i;
    centroids(i, :) = (idx_i' * X) ./ sum(idx_i);
end

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

end

author: zlotus


text:
function g = onl_MetaCtrlPnlMiscOpts(varargin)

g = arg_define(0,varargin,...
        arg({'dispBenchmark','DisplayBenchmarking'},false,[],'Display Benchmarking','type','logical'), ...
        arg({'winLenSec','DataWindowLength'},2,[0 Inf],'EEG chunk size (sec)'), ...
        arg_subtoggle({'dispCSD','DisplayCSD'},[],@vis_csd,'Display Current Density. This only takes effect if the SourceLocalization filter is enabled'), ...
        arg_subtoggle({'doSIFT','RunSIFT'},[], ...
        {...
            arg({'channelSubset','ChannelSubset'},[], [], 'Cell array of channel names to retain.','type','cellstr','shape','row'), ...
            arg({'dispSpectrum','DisplaySpectrum'},false,[],'Display Spectral Density'), ...
            arg_subtoggle({'dispBrainMovie','DisplayBrainMovie'},[],...
            { ...
                arg({'adaptBrainMovieLimits','AdaptBrainMovieLimits'},true,[],'Adapt Color and Edge Scaling of Brainmovie 3D'), ...
                arg_subtoggle({'colorCortexByRoi','ColorCortexByROI'},{},@hlp_getROIVertexColorTable,'Color cortical surface by ROI'), ...
            }, ...
            'Display BrainMovie3D') ...
        },'Enable Source Information Flow Toolbox'), ...
        arg({'streamLSL','StreamLSL'},true,[],'Stream results via LSL'), ...
        arg({'suppress_console_output','SuppressConsoleOutput'},true,[],'Suppress most output to console') ...
        );


text:
function planC = dose2CERR(doseNew,doseError,fractionGroupID,doseEdition,description,register,regParamsS,overWrite,assocScanUID,planC)
%function dose2CERR(doseNew,doseError,fractionGroupID,doseEdition,description,register,regParamsS,overWrite,planC)
%Use this function to put a new dose distribution into CERR.
%See test_dose2CERR in the putNewDose subdirectory for an example of its use.
%doseNew - 3-D array of dose values.
%doseError - estimate of standard deviation for Monte Carlo plans (can be left empty if not applicable, i.e. [])
%fractionGroupID - unique ID as a function of the number of fractions.
%doseEdition - another ID, see the RTOG specification.
%description - Short description string if desired.
%register - 'CT' or blank (defaults to CT), or non-CT.  If CT, registration geomtrical
%information is taken from the CT scan.  Otherwise, registration data is taken from regParamsS:
%regParamsS should contain geometric registration data including the following fields:
%regParamsS.horizontalGridInterval = 0.2 (say)  (x voxel width)
%regParamsS.verticalGridInterval   = 0.2 (say)  (y voxel width)
%regParamsS.coord1OFFirstPoint     = 0.5 (say)  (x value of center of upper left voxel on all slices)
%regParamsS.coord2OFFirstPoint     = 25  (say)  (y value of center of upper left voxel on all slices
%regParamsS.zValues                = [0.5 1.0 1.5 2.0 ...] (say) (z values of all slices)
%
%If overWrite = 'yes', overwrites the last written dose distribution.
%(x,y,z) are in the AAPM/RTOG coordinate system.
%El Naqa & Deasy, Feb 03.
%Latest Modifications: 7 March 03, JOD; input can be non-registered to CT scan.  Now using regParamsS.
%                      31 July 03, JOD, corrected use of fraction group ID.
%                       5  Nov 03, JOD, added overwrite flag
%                      12  Oct 06, APA, corrected horizontal/vertical grid interval assignment
%
% 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/>.

if ~exist('planC','var')
    global planC
end
indexS=planC{end};
W = size(doseNew);
if ismatrix(doseNew) %Single slice
    temp = nan(W(1),W(2),3);
    temp(:,:,2) = doseNew;
    doseNew = temp;
    W(3) = 3;
end

%How many old dose distributions?
prevSets = length(planC{indexS.dose});

setIndex = prevSets + 1;

if nargin > 7 & strcmpi(overWrite,'yes')
  setIndex = prevSets;
end

%Get latest dose structure:
% [headerInitS, commentInitS, scanInitS, structureInitS, beamGeometryInitS, beamsInitS, digitalFilmInitS, doseInitS, ...
%            DVHInitS, seedGeometryInitS, RTTreatmentInitS, importLogInitS, planInitC, garbage] = initializeCERR;
doseInitS = initializeCERR('dose');
minDose = min(doseNew(:));
if minDose >= 0
    doseInitS(1).doseArray = doseNew;
else
    doseInitS(1).doseArray = doseNew - minDose;
    doseInitS(1).doseOffset = -minDose;
end
doseInitS(1).imageNumber=  '';
doseInitS(1).imageType='DOSE';
doseInitS(1).caseNumber=1;
doseInitS(1).patientName='';
doseInitS(1).doseNumber= setIndex;
doseInitS(1).doseType='PHYSICAL';
doseInitS(1).doseUnits='GRAYS';
doseInitS(1).doseScale=1;
doseInitS(1).fractionGroupID= fractionGroupID;
doseInitS(1).numberOfTx= '';
doseInitS(1).orientationOfDose='TRANSVERSE';
doseInitS(1).numberRepresentation=  '';
doseInitS(1).numberOfDimensions=length(W);
doseInitS(1).sizeOfDimension1=W(2);
doseInitS(1).sizeOfDimension2=W(1);
doseInitS(1).sizeOfDimension3=W(3);
if ~isempty(description)
  doseInitS(1).doseDescription=description;
end
if ~isempty(doseEdition)
  doseInitS(1).doseEdition= doseEdition;
end
doseInitS(1).unitNumber='';
doseInitS(1).writer='';
doseInitS(1).dateWritten='';
doseInitS(1).planNumberOfOrigin='';
doseInitS(1).planEditionOfOrigin='';
doseInitS(1).studyNumberOfOrigin='';
doseInitS(1).versionNumberOfProgram='';
doseInitS(1).xcoordOfNormaliznPoint='';
doseInitS(1).ycoordOfNormaliznPoint='';
doseInitS(1).zcoordOfNormaliznPoint='';
doseInitS(1).doseAtNormaliznPoint='';
doseInitS(1).coord3OfFirstPoint= '';   %Not used in CERR
doseInitS(1).depthGridInterval = '';
doseInitS(1).assocScanUID = assocScanUID;
doseInitS(1).doseUID = createUID('DOSE');

%Monte Carlo specific:
if ~isempty(doseError)
  doseInitS(1).doseError=doseError;
end

numDoses = length(planC{indexS.dose});
for i = 1 : numDoses
  origS = planC{indexS.dose}(i);
  tmpS(i) = backFillFieldNames(doseInitS,origS);
end

if numDoses > 0
  planC{indexS.dose} = tmpS;
end

if nargin > 5 & strcmpi(register,'CT')

    scanNum = getAssociatedScan(assocScanUID,planC);

    grid2Units = planC{indexS.scan}(scanNum).scanInfo(1).grid2Units;
    grid1Units = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
    doseInitS(1).horizontalGridInterval = grid2Units;
    doseInitS(1).verticalGridInterval= - abs(grid1Units);

    abGrid1Units = abs(grid1Units);
    abGrid2Units = abs(grid2Units);

    xOffset = planC{indexS.scan}(scanNum).scanInfo(1).xOffset;
    yOffset = planC{indexS.scan}(scanNum).scanInfo(1).yOffset;

    CTWidth = planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension2;
    CTHeight = planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension1;

    doseInitS(1).coord1OFFirstPoint=  xOffset - (CTWidth/2) * abGrid2Units + abGrid2Units/2;
    doseInitS(1).coord2OFFirstPoint=  yOffset + (CTHeight/2) * abGrid1Units - abGrid1Units/2;

    % get from CT info:

    zValues = [planC{indexS.scan}(scanNum).scanInfo(:).zValue];
    doseInitS(1).zValues = zValues;
    doseInitS(1).delivered='';

elseif nargin > 6 & strcmpi(register,'UniformCT')

    scanNum = getAssociatedScan(assocScanUID,planC);

    uniformInfoS = planC{indexS.scan}(scanNum).uniformScanInfo;
    %[CTUniform3D, uniformInfoS] = getUniformizedCTScan;

    grid2Units = uniformInfoS.grid2Units;
    grid1Units = uniformInfoS.grid1Units;
    doseInitS(1).horizontalGridInterval = grid2Units;
    doseInitS(1).verticalGridInterval= - abs(grid1Units);

    abGrid1Units = abs(grid1Units);
    abGrid2Units = abs(grid2Units);

    xOffset = uniformInfoS.xOffset;
    yOffset = uniformInfoS.yOffset;

    CTWidth = uniformInfoS.sizeOfDimension2;
    CTHeight = uniformInfoS.sizeOfDimension1;

    doseInitS(1).coord1OFFirstPoint=  xOffset - (CTWidth/2) * abGrid2Units + abGrid2Units/2;
    doseInitS(1).coord2OFFirstPoint=  yOffset + (CTHeight/2) * abGrid1Units - abGrid1Units/2;

    % get from CT info:

    %sizeArray = getUniformizedSize(planC)
    sizeArray = getUniformScanSize(planC{indexS.scan}(scanNum));
    numSlices = sizeArray(3);
%    numSlices = size(CTUniform3D,3);

    zFirst = uniformInfoS.firstZValue;

    sliceThickness = uniformInfoS.sliceThickness;

    zLast = (numSlices - 1) * sliceThickness + zFirst;

    zValues = zFirst : sliceThickness: zLast;

    doseInitS(1).zValues = zValues;
    doseInitS(1).delivered='';

elseif nargin > 6 & ~strcmpi(register,'UniformCT') & ~strcmpi(register,'CT')

  doseInitS(1).horizontalGridInterval = regParamsS.horizontalGridInterval;
  doseInitS(1).verticalGridInterval   = regParamsS.verticalGridInterval;
  doseInitS(1).coord1OFFirstPoint     = regParamsS.coord1OFFirstPoint;
  doseInitS(1).coord2OFFirstPoint     = regParamsS.coord2OFFirstPoint;
  doseInitS(1).zValues                = regParamsS.zValues;

else

  error('Inputs to dose2CERR are incorrect')

end

% %Update the CERR dose menu
% try
%   putDoseMenu(stateS.handle.CERRSliceViewer, planC, indexS);
% end

planC{indexS.dose} = dissimilarInsert(planC{indexS.dose},doseInitS, setIndex,[]);  %This way future field additions are automatically included.

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function show_mov(in, in2)
    % ZMAP script show_map.m. Creates Dialog boxes for Z-map calculation
    % does the calculation and makes displays the map
    % stefan wiemer 11/94
    %
    % make dialog interface and call maxzlta
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    % Input Rubberband
    %
    report_this_filefun();

    if in2 ~= 'calma'

        %initial values
        nustep = 10;
        ZG.compare_window_dur_v3 = years(1.5);
        it = t0b +1;
        figure(mess);
        clf
        set(gca,'visible','off')
        set(gcf,'Units','pixel','NumberTitle','off','Name','Input Parameters');

        set(gcf,'pos',[ ZG.welcome_pos, ZG.welcome_len +[200, -50]]);

        % creates a dialog box to input some parameters
        %

        inp2_field=uicontrol('Style','edit',...
            'Position',[.80 .80 .18 .15],...
            'Units','normalized','String',num2str(nustep),...
            'callback',@callbackfun_001);

        txt2 = text(...
            'Position',[0. 0.9 0 ],...
            'FontWeight','bold',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'String','Please input Number of Frames:');

        if in == 'rub' | in == 'lta'

            txt3 = text(...
                'Position',[0. 0.65 0 ],...
                'FontWeight','bold',...
                'FontSize',ZmapGlobal.Data.fontsz.m ,...
                'String','Please input window length in years (e.g. 1.5):');
            inp3_field=uicontrol('Style','edit',...
                'Position',[.80 .575 .18 .15],...
                'Units','normalized','String',num2str(years(ZG.compare_window_dur)),...
                'callback',@callbackfun_002);

        end   % if in = rub

        close_button=uicontrol('Style','Pushbutton',...
            'Position', [.60 .05 .15 .15 ],...
            'Units','normalized','Callback',@(~,~)close(),'String','Cancel');

        go_button=uicontrol('Style','Pushbutton',...
            'Position',[.25 .05 .15 .15 ],...
            'Units','normalized',...
            'callback',@callbackfun_003,...
            'String','Go');

        set(gcf,'visible','on');watchoff

        % do the calculations:
        %

    else     % if in2 ~=calma

        % check if time are with limits
        %

        % initial parameter
        winlen_days = ZG.compare_window_dur/ZG.bin_dur;
        ti = (it -t0b)/days(ZG.bin_dur);
        var1 = zeros(1,ncu);
        var2 = zeros(1,ncu);
        mean1 = zeros(1,ncu);
        mean2 = zeros(1,ncu);
        as = zeros(1,ncu);
        [len, ncu] = size(cumuall); len = len-2;
        len = len -2;
        step = len/nustep;

        % loop over all frames

        j = 0;
        figure
        rect = [0.10 0.30 0.55 0.50 ];
        rect1 = rect;
        axes('position',rect1)
        axis('off')
        m = moviein(length(1:step:len-winlen_days));
        for ti = winlen_days:step:len-winlen_days
            j = j+1;
            var1 = zeros(1,ncu);
            var2 = zeros(1,ncu);
            mean1 = zeros(1,ncu);
            mean2 = zeros(1,ncu);
            as = zeros(1,ncu);

            % loop over all grid points for percent
            %
            %
            if in =='per'

                for i = 1:ncu
                    mean1(i) = mean(cumuall(1:ti,i));
                    mean2(i) = mean(cumuall(ti:len,i));
                end    %for i
                as = -((mean1-mean2)./mean1)*100;

                strib = 'Change in Percent';
                stri2 = ['ti=' num2str(ti*days(ZG.bin_dur) + t0b)  ];

            end  % if in = = per

            % loop over all point for rubber band
            %
            if in =='rub'

                for i = 1:ncu
                    mean1(i) = mean(cumuall(1:ti,i));
                    mean2(i) = mean(cumuall(ti+1:ti+winlen_days,i));
                    var1(i) = cov(cumuall(1:ti,i));
                    var2(i) = cov(cumuall(ti+1:ti+winlen_days,i));
                end %  for i ;
                as = (mean1 - mean2)./(sqrt(var1/ti+var2/winlen_days));

            end % if in = rub

            % make the AST function map
            if in =='ast'
                for i = 1:ncu
                    mean1(i) = mean(cumuall(1:ti,i));
                    var1(i) = cov(cumuall(1:ti,i));
                    mean2(i) = mean(cumuall(ti+1:len,i));
                    var2(i) = cov(cumuall(ti+1:len,i));
                end    %for i
                as = (mean1 - mean2)./(sqrt(var1/ti+var2/(len-ti)));
            end % if in = ast

            if in =='lta'
                disp('Calculate LTA')
                %cu = [cumuall(1:ti-1,:) ; cumuall(ti+winlen_days+1:len,:)];
                mean1 = mean([cumuall(1:ti-1,:) ; cumuall(ti+winlen_days+1:len,:)]);
                mean2 = mean(cumuall(ti:ti+winlen_days,:));
                for i = 1:ncu
                    var1(i) = cov([cumuall(1:ti-1,i) ; cumuall(ti+winlen_days+1:len,i)]);
                    var2(i) = cov(cumuall(ti:ti+winlen_days,i));
                end     % for i
                as = (mean1 - mean2)./(sqrt(var1/(len-winlen_days)+var2/winlen_days));
            end % if in = lta

            normlap1=nan(length(tmpgri(:,1)),1)
            normlap2=nan(length(tmpgri(:,1)),1)
            normlap2(ll)= as(:);
            %construct a matrix for the color plot
            valueMap=reshape(normlap2,length(yvect),length(xvect));

            %plot imge
            % set values gretaer ZG.tresh_km = nan
            %
            re4 = valueMap;
            [len, ncu] = size(cumuall);
            [n1, n2] = size(cumuall);
            s = cumuall(n1,:);
            normlap2(ll)= s(:);
            r=reshape(normlap2,length(yvect),length(xvect));
            re4(r > ZG.tresh_km) = nan;

            orient landscape
            set(gcf,'PaperPosition',[ 0.1 0.1 8 6])
            axes('position',rect1)
            set(gca,'NextPlot','add')
            pco1 = pcolor(gx,gy,re4);
            caxis([ZG.minc ZG.maxc]);
            axis([ s2_west s1_east s4_south s3_north])
            set(gca,'NextPlot','add')
            %overlay
            if in == 'ast'
                tx2 = text(0.07,0.85 ,['AS; t=' num2str(ti*days(ZG.bin_dur)+t0b)  ] ,...
                    'Units','Norm','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k','FontWeight','bold');
            end

            if in == 'lta'
                tx2 = text(0.07,0.85 ,['LTA; t=' num2str(ti*days(ZG.bin_dur)+t0b)  ] ,...
                    'Units','Norm','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k','FontWeight','bold');
            end

            if in == 'rub'
                tx2 = text(0.07,0.85 ,['RUB; t=' num2str(ti*days(ZG.bin_dur)+t0b)  ] ,...
                    'Units','Norm','FontSize',ZmapGlobal.Data.fontsz.m,'Color','k','FontWeight','bold');
            end

            set(gca,'FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
                'FontWeight','bold','LineWidth',1.5,...
                'Box','on','SortMethod','childorder')

            shading interp
            has = gca;
            disp('now getting frame...')
            m(:,j) = getframe(has);
            delete(gca);
            delete(gca);
            delete(gca)
            fs_m = get(gcf,'pos');

        end  % loop over frames

        close(gcf)

        showmovi
    end   % if calma ~| in2

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        nustep=str2double(inp2_field.String);
        inp2_field.String=num2str(nustep);
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ZG.compare_window_dur=years(str2double(mysrc.String));
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        nustep=str2num(inp2_field.String);
        ZG.compare_window_dur=years(str2num(inp3_field.String));
        watchon;
        drawnow;
        in2 = 'calma';
        fixaxmo;
    end

end

author: CelsoReyes


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb = compactbit(b)
%
% b = bits array
% cb = compacted string of bits (using words of 'word' bits)

[nSamples nbits] = size(b);
nwords = ceil(nbits/8);
cb = zeros([nSamples nwords], 'uint8');

for j = 1:nbits
    w = ceil(j/8);
    cb(:,w) = bitset(cb(:,w), mod(j-1,8)+1, b(:,j));
end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function U = spm_mvb_U(Y,priors,X0,xyz,vox,nu)
% Constructs patterns U for Multivariate Bayesian inversion of a linear model
% FORMAT U = spm_mvb_U(Y,priors,X0,xyz,vox,nu)
% Y      - data-feature matrix
% priors - 'null'      % no patterns
%        - 'compact'   % reduced (ns/3); using SVD on local compact support
%        - 'sparse'    % a pattern is a voxel
%        - 'smooth'    % patterns are local Gaussian kernels
%        - 'singular'  % patterns are global singular vectors
%        - 'support'   % the patterns are the images
%
% X0     - confounds
% xyz    - location in mm
% vox    - voxel size in mm
% nu     - number of patterns (for 'compact')
%
% U      - pattern or mode weights
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mvb_U.m 7654 2019-08-25 20:09:35Z karl $

% defaults
%--------------------------------------------------------------------------
try, X0;  catch, X0  = [];   end
try, xyz; catch, xyz = [];   end

% get orders
%--------------------------------------------------------------------------
ns     = size(Y,1);                      % number of samples
nv     = size(Y,2);                      % number of voxels
if nargin < 6
    nu = min(ceil(ns/3),nv);             % number of patterns
end

% confounds
%--------------------------------------------------------------------------
if isempty(X0); X0 = zeros(ns,1); end

% get U: X = Y*P + X0*Q + R
%        P = U*E;
%--------------------------------------------------------------------------
% assemble empirical priors
%==========================================================================
switch priors

    case 'null'
        %------------------------------------------------------------------
        U     = sparse(nv,0);

    case 'sparse'
        %------------------------------------------------------------------
        U     = speye(nv,nv);

    case 'smooth'
        %------------------------------------------------------------------
        sm    = 4;                            % smoothness fixed at 4mm std
        dlim  = (4*sm)^2;                     % spatial limit (mm)^2
        s     = sm^2;                         % Smoothness variance
        xyz   = xyz';
        nlim  = 256;                          % voxel limit
        Vvx   = prod(vox);                    % volume of a voxel
        Vlr   = 4/3*pi*dlim^(3/2);            % VOI around a voxel
        Nlr   = round(Vlr/Vvx*.85);           % # of voxel in VOI; keep 85%
        U     = spalloc(nv,nv,nv*Nlr);        % pre-allocate memory
        fprintf('Creating smooth patterns - please wait\n')
        kk    = floor(nv/4);
        fprintf('\n0%%')
        for i = 1:nv
            if ~rem(i,kk), fprintf('.....%2i%%',25*i/kk); end
            u      = (xyz(:,1) - xyz(i,1)).^2;
            j      = find(u < dlim);
            u(j)   = u(j) + (xyz(j,2) - xyz(i,2)).^2;
            j      = j((u(j) < dlim));
            u(j)   = u(j) + (xyz(j,3) - xyz(i,3)).^2;
            j      = j((u(j) < dlim));
            if length(j)>nlim
                [q,k]  = sort(u(j));
                k      = k(1:nlim);
                j      = j(k);
            end
            U(j,i) = exp(-u(j)/(2*s));
        end
        fprintf('\nThank you\n')

    case 'singular'

        % get kernel (singular vectors)
        %------------------------------------------------------------------
        Y       = Y - X0*(pinv(X0)*Y);         % remove confounds
        [u,s,v] = spm_svd(Y,1/4);              % c.f., Kaiser criterion
        U       = v/s;

    case 'support'

        % get kernel (image vectors)
        %------------------------------------------------------------------
        if nv > ns
            R  = speye(size(X0,1)) - X0*pinv(X0);
            U  = (R*Y)';
        else
            U  = speye(nv,nv);
        end

    case 'compact'

        % get kernel (compact vectors)
        %------------------------------------------------------------------
        nc    = max(fix(nv/nu),1);           % voxels in compact support
        X0    = spm_svd(X0);                 % confounds
        Y     = Y - X0*(X0'*Y);              % remove confounds
        C     = sum(Y.^2);                   % variance of Y
        U     = spalloc(nv,nu,nc*nu);
        J     = 1:nv;
        for i = 1:nu

            % find maximum variance voxel
            %--------------------------------------------------------------
            [v,j] = max(C);
            d     = 0;
            for k = 1:size(xyz,1)
               d  = d + (xyz(k,:) - xyz(k,j)).^2;
            end
            [d,j] = sort(d);
            try
                j = j(1:nc);
            end

            % save principal eigenvector
            %--------------------------------------------------------------
            k        = J(j);
            u        = spm_svd(Y(:,k)');
            U(k,i)   = u(:,1);

            % remove compact support voxels and start again
            %--------------------------------------------------------------
            J(j)     = [];
            C(j)     = [];
            xyz(:,j) = [];

        end

    otherwise
        disp('unknown prior')
        return
end

author: spm


text:
function res = get_release_mode()

global POSE_RELEASE_MODE;

if isempty(POSE_RELEASE_MODE)
    POSE_RELEASE_MODE = false;
end

res = POSE_RELEASE_MODE;

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


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

clear a1 a2 a3 a4 a5
a1=zeros(size(newt2));
a2=zeros(size(newt2));
a3=zeros(size(newt2));
a4=zeros(size(newt2));
a5=zeros(size(newt2));
a6=zeros(size(newt2));
bins=160
c1=0.9;
c2=0.9;
c3=0.9;
a1=newt2;
nmb=sum(a1(:,6)==0);
n=size(a1,1);
str1=['Anteil mit Mb=0: ' num2str(nmb/n*100) '%'];
disp(str1)

disp(' ')
exfig=figure_exists('Mb=f(Ms)',0);

if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Mb=f(Ms)',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end

hold on
[p,r]=corrmbms2(a1(:,10),a1(:,6),1,'Ms','Mb');
hold off

a2=a1;

for i=1:size(a2,1)
    if and(a2(i,6)==0,not(a2(i,10)==0));
        a2(i,6)=round(10*(a2(i,10)*p(1)+p(2)))/10;
        a4(i,1)=a2(i,6);
    end
end

nmb=sum(a2(:,6)==0);
n=size(a2,1);
str1=['Anteil mit Mb=0: ' num2str(nmb/n*100) '%'];
disp(str1)
disp(' ')

for i=1:size(a2,1)
    if a2(i,11)>10;
        a2(i,11)=0;
    end
end

exfig=figure_exists('Mb=f(???)',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits('Name','Mb=f(???)',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end

hold on
[p,r]=corrmbms2(a2(:,11),a2(:,6),1,'???','mb');
hold off

a3=a2;

for i=1:size(a3,1)
    if and(a3(i,6)==0,not(a3(i,11)==0));
        a3(i,6)=round(10*(a3(i,11)*p(1)+p(2)))/10;
        a5(i,1)=a3(i,6);
    end
end

nmb=sum(a3(:,6)==0);
n=size(a3,1);
str1=['Anteil mit Mb=0: ' num2str(nmb/n*100) '%'];
disp(str1)
disp(' ')
disp('Unbenutzbare Eintr�ge')
str1=num2str(sum((a.Magnitude==0&a(:,10)==0)&a(:,11)==0)+sum((a.Magnitude==0&a(:,10)==0)&a(:,11)>10));
disp(str1)

bins=0.05:.1:9;

exfig=figure_exists('Mb Histogram',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Mb Histogram',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end
clf
%subplot(2,3,1)
lo=not(a1(:,6)==0);
ah=a1(lo,:);
bar(bins, histc(ah(:,6),bins))
Xlim([1 9.5])
Xlabel('Mb')

exfig=figure_exists('Ms Histogram',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Ms Histogram',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end
clf
%subplot(2,3,2)
lo=a1(:,6)==0&not(a1(:,10)==0);
ah=a1(lo,:);
bar(bins, histc(ah(:,10),bins))
Xlim([1 9.5])
Xlabel('Ms')

exfig=figure_exists('Mu Histogram',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Mu Histogram',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end
clf
%subplot(2,3,3)
lo=a2(:,6)==0&a2(:,10)==0&not(a2(:,11)==0);
ah=a1(lo,:);
bar(bins, histc(ah(:,11),bins))
Xlim([1 9.5])
Xlabel('Mu')

exfig=figure_exists('Mb + Histogram',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Mb + Histogram',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end
clf
%subplot(2,3,4)
bar( bins,histc(a3(:,6),bins))
Xlim([1 9.5])
Xlabel('Mb +')

exfig=figure_exists('Mb aus Ms Histogram',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Mb aus Ms Histogram',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end
clf
%subplot(2,3,5)
bar( bins,histc(a4(:,1),bins))
Xlabel('Mb aus Ms')
Xlim([1 9.5])

exfig=figure_exists('Mb aus Mu Histogram',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Mb aus Mu Histogram',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end
clf
%subplot(2,3,6)
bar( bins,histc(a5(:,1),bins))
Xlabel('Mb aus Mu')
Xlim([1 9.5])

%a1=a3;
%clear('a2')
%clear('a3')


text:
classdef SparsePyrLKOpticalFlow < handle
    %SPARSEPYRLKOPTICALFLOW  Class used for calculating a sparse optical flow
    %
    % The class can calculate an optical flow for a sparse feature set using
    % the iterative Lucas-Kanade method with pyramids.
    %
    % See also: cv.SparsePyrLKOpticalFlow.calc, cv.calcOpticalFlowPyrLK
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % default [21,21]
        WinSize
        % default 3
        MaxLevel
        % default `struct('type','Count+EPS', 'maxCount',30, 'epsilon',0.01)`
        TermCriteria
        % default 0
        Flags
        % default 1e-4
        MinEigThreshold
    end

    %% Constructor/destructor
    methods
        function this = SparsePyrLKOpticalFlow()
            %SPARSEPYRLKOPTICALFLOW  Creates instance of SparsePyrLKOpticalFlow
            %
            %     obj = cv.SparsePyrLKOpticalFlow()
            %
            % See also: cv.SparsePyrLKOpticalFlow.calc
            %
            this.id = SparsePyrLKOpticalFlow_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.SparsePyrLKOpticalFlow
            %
            if isempty(this.id), return; end
            SparsePyrLKOpticalFlow_(this.id, 'delete');
        end
    end

    %% SparseOpticalFlow
    methods
        function varargout = calc(this, prevImg, nextImg, prevPts, varargin)
            %CALC  Calculates a sparse optical flow
            %
            %     nextPts = obj.calc(prevImg, nextImg, prevPts)
            %     [nextPts, status, err] = obj.calc(...)
            %     [...] = obj.calc(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __prevImg__ First input image.
            % * __nextImg__ Second input image of the same size and the same
            %   type as `prevImg`.
            % * __prevPts__ Vector of 2D points for which the flow needs to be
            %   found.
            %
            % ## Output
            % * __nextPts__ Output vector of 2D points containing the
            %   calculated new positions of input features in the second image.
            % * __status__ Output status vector. Each element of the vector is
            %   set to 1 if the flow for the corresponding features has been
            %   found. Otherwise, it is set to 0.
            % * __err__ Optional output vector that contains error response
            %   for each point (inverse confidence).
            %
            % ## Options
            % * __InitialFlow__ Vector of 2D points to be used for the initial
            %   estimate of `nextPts`. Not set by default.
            %
            % See also: cv.SparsePyrLKOpticalFlow, cv.calcOpticalFlowPyrLK
            %
            [varargout{1:nargout}] = SparsePyrLKOpticalFlow_(this.id, 'calc', prevImg, nextImg, prevPts, varargin{:});
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.SparsePyrLKOpticalFlow.empty
            %
            SparsePyrLKOpticalFlow_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm is empty (e.g. in the very
            %   beginning or after unsuccessful read).
            %
            % See also: cv.SparsePyrLKOpticalFlow.clear
            %
            b = SparsePyrLKOpticalFlow_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.SparsePyrLKOpticalFlow.save, cv.SparsePyrLKOpticalFlow.load
            %
            name = SparsePyrLKOpticalFlow_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in a file storage.
            %
            % See also: cv.SparsePyrLKOpticalFlow.load
            %
            SparsePyrLKOpticalFlow_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from a file storage.
            % The previous model state is discarded.
            %
            % See also: cv.SparsePyrLKOpticalFlow.save
            %
            SparsePyrLKOpticalFlow_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.WinSize(this)
            value = SparsePyrLKOpticalFlow_(this.id, 'get', 'WinSize');
        end
        function set.WinSize(this)
            SparsePyrLKOpticalFlow_(this.id, 'set', 'WinSize', value);
        end

        function value = get.MaxLevel(this)
            value = SparsePyrLKOpticalFlow_(this.id, 'get', 'MaxLevel');
        end
        function set.MaxLevel(this)
            SparsePyrLKOpticalFlow_(this.id, 'set', 'MaxLevel', value);
        end

        function value = get.TermCriteria(this)
            value = SparsePyrLKOpticalFlow_(this.id, 'get', 'TermCriteria');
        end
        function set.TermCriteria(this)
            SparsePyrLKOpticalFlow_(this.id, 'set', 'TermCriteria', value);
        end

        function value = get.Flags(this)
            value = SparsePyrLKOpticalFlow_(this.id, 'get', 'Flags');
        end
        function set.Flags(this)
            SparsePyrLKOpticalFlow_(this.id, 'set', 'Flags', value);
        end

        function value = get.MinEigThreshold(this)
            value = SparsePyrLKOpticalFlow_(this.id, 'get', 'MinEigThreshold');
        end
        function set.MinEigThreshold(this)
            SparsePyrLKOpticalFlow_(this.id, 'set', 'MinEigThreshold', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%GETUSER Get the user field of a dataset
%
%    USERFIELD = GETUSER(A,FIELD)
%
% INPUT
%   A          Dataset
%   FIELD      Character string: name of structure field of USER field of A.
%
% OUTPUT
%   USERFIELD  Requested field
%
% If FIELD is is the empty string ('') the entire user field is returned.
%
% If the requested field does not exist USERFIELD = [];
%
% Note the the USER field of datasets was originally intended for a user
% defined description of datasets. Later its usage was extended to a field
% for storing general information on datasets. For that reason 'old'
% datasets without a structure in the user field are transformed such that
% this information is stored in a subfield USER in the user field. It can
% be retrieved by GETUSER(A).


text:
function [F] = CutOffFilter2D(Nc,frac)

% function [F] = CutOffFilter2D(Nc,frac)
% Purpose : Initialize 2D cut off filter matrix of order Norderin

Globals2D;

filterdiag = ones(Np,1);

% build exponential filter
sk = 1;
for i=0:N
  for j=0:N-i
    if (i+j>=Nc)
      filterdiag(sk) = frac;
    end
    sk = sk+1;
  end
end

F = V*diag(filterdiag)*invV;
return;

repo: Aero-matlab


text:
classdef CoreErrorUtilities < handle
    % CoreErrorUtilities. Error-related utility functions.
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    methods (Static)

        % Obtain the name of the function which signalled the error (useful in
        % error reporting)
        function [calling_function, stack] = GetCallingFunction(levels_to_ignore)
            caller_found = false;
            max_levels = 10;
            calling_function = [];

            while ~caller_found
                % Get the call stack, excluding the call to this function and the
                % caller of this function
                [stack, ~] = dbstack(levels_to_ignore, '-completenames');
                if isempty(stack)
                    caller_found = true;
                else
                    stack_top = stack(1);
                    calling_function = stack_top.name;
                    if length(calling_function) < 11 || ~strcmp(calling_function(1:11), 'CoreReporting')
                        caller_found = true;
                    end
                    levels_to_ignore = levels_to_ignore + 1;
                end
            end
        end

        function ThrowException(id, message)
            [stack, ~] = dbstack(2, '-completenames');

            msgStruct = [];
            msgStruct.message = message;
            msgStruct.identifier = id;
            msgStruct.stack = stack;
            error(msgStruct);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function model2text(data,model)
% This script takes in a reconstruction model and prints it out into a text
% file. It can be used standalone from the `rBioNet`.
%
% USAGE:
%
%    model2text(data, model)
%
% INPUTS:
%    data:    model or .mat datafile table as in ReconstructionCreator (RC)
%    model:   true if data is a model, false if RC table.
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is
if model
    data = model2data(data,1);
end

[FileName, PathName] = uiputfile('.txt');
if FileName == 0
    return;
end

fid = fopen(fullfile(PathName,FileName),'w');
rxns = data{1};
mets = data{5};
desc = data{2};

%	#Model - the model fields are in the following order. Number of rows are fixed.
%		Name
%		Organism
%		Author
%		Notes
%		Geneindex
%		Geneindex notes

%desc_fields = {'% Name','%Organism','Author','Notes','Geneindex','Geneindex notes'};
fprintf(fid,...
    ['%% This is a COBRA reconstruction model exported to a text file \n'...
    '%% with the rBioNet. These structures are intended to move models\n'...
    '%% between rBioNet in Matlab and the web version. Currently under\n'...
    '%% construction. June 28 2011, Stefan Gretar Thorleifsson.\n'...
    '%%\n'...
    '%% \t # - denotes a specific field is starting description/reaction/metabolites\n'...
    '%% \t %% - denotes comments\n']);

fprintf(fid,'%%\n%%\n%%\n');

fprintf(fid,['#Description\n'...
    'Name\t%s\n'...
    'Organism\t%s\n'...
    'Author\t%s\n'...
    'Notes\t%s\n'...
    'Geneindex_File\t%s\n'...
    'Geneindex_Date\t%s\n'...
    'Geneindex_Source\t%s\n'],...
    desc{1}, desc{2}, desc{3}, desc{4}, desc{5}, desc{6}, desc{7});

% for i = 1:length(desc)
%     fprintf(fid,'%s\n',desc{i});
% end

rxn_fields = {'% Abbreviation',' Description ', 'Formula', 'Reversible', 'GPR' ,...
    'LB' , 'UB' , 'Mechanism Confidence Score', 'Subsystem', 'References',...
    'Notes', 'EC number', 'KeggID'};

fprintf(fid,'%%\n%%\n%%\n');

for i = 1:length(rxn_fields)
    fprintf(fid, '%s\t',rxn_fields{i});
end
fprintf(fid,'\n#Reactions\n');

for i = 1:size(rxns,1)
    %Cannot print logical
    if rxns{i,5}
        rxns{i,5} = '1';
    else
        rxns{i,5} = '0';
    end
    rxns{i,7} = num2str(rxns{i,7});
    rxns{i,8} = num2str(rxns{i,8});
    rxns{i,9} = num2str(rxns{i,9});
    if ~isempty(rxns{i,3})
        rxns{i,3} = regexprep(rxns{i,3},'\t', '');
    end
    if ~isempty(rxns{i,12})
        rxns{i,12} = regexprep(rxns{i,12},'\t', '');
    end

    if ~isempty(rxns{i,13})
        rxns{i,13} = regexprep(rxns{i,13},'\t', '');
    end
     if ~isempty(rxns{i,14})
        rxns{i,14} = regexprep(rxns{i,14},'\t', '');
     end

    fprintf(fid, '%s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\n',...
        rxns{i,2}, rxns{i,3}, rxns{i,4}, rxns{i,5}, rxns{i,6},...
        rxns{i,7}, rxns{i,8}, rxns{i,9}, rxns{i,10}, rxns{i,11}, rxns{i,12},...
        rxns{i,13}, rxns{i,14});
end

%	Abbreviation | Description | Neutral formula | Charged formula | Charge | KeggID | PubChemID |
%	CheBIID | Inchi String | Smile | HMDB

met_fields = {'%Abbreviation', 'Description', 'Neutral formula', 'Charged formula',...
    'Charge', 'KeggID','PubChemID', 'CheBlID', 'Inchi String', 'Smile', 'HMDB',...
    'metHepatoNetID','metEHMNID'};

fprintf(fid,'%%\n%%\n%%\n');

for i = 1:length(met_fields)
    fprintf(fid, '%s\t',met_fields{i});
end
fprintf(fid,'\n#Metabolites\n');

for i = 1:size(mets,1)

    if isa(mets{i,4},'cell')
        a = mets{i,4};
        mets{i,4} = a{1};
    end
    mets{i,4} = num2str(mets{i,4});

    fprintf(fid, '%s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \n',...
        mets{i,1}, mets{i,2}, '',mets{i,3}, mets{i,4}, mets{i,5}, mets{i,6},...
        mets{i,7}, mets{i,8}, mets{i,9},mets{i,10}, mets{i,11}, mets{i,12});
end

fclose(fid);


text:
function [distri]=Cell_follicledistribution(segout, I1,coeff, n)
if(~exist('n','var'))
    n=13;
end
[boundary(:,2),boundary(:,1)]=find(segout==255);

[fo(:,2), fo(:,1)]=find(I1==1);
theta1=orderpoints(boundary, mean(boundary), coeff);
[trash, index1]=sort(theta1,'ascend');
theta1=theta1(index1);
boundary=boundary(index1,:);

theta2=orderpoints(fo, mean(boundary), coeff);
[trash, index2]=sort(theta2,'ascend');
theta2=theta2(index2);
fo=fo(index2,:);

center=mean(boundary);
KK=500;
interval=linspace(0,2*pi,n);
inter=1:floor(n/4):n;
figure;
imshow(I1)

for i=1:n-1
    hold on
    index_boundary=find(theta1>interval(i)&theta1<=interval(i+1));
    index_focell=find(theta2>interval(i)&theta2<=interval(i+1));
    length(i)=level_center_length(boundary(index_boundary,:)',size(boundary,2));
    number_focell(i)=numel(index_focell);
    temp=coeff*[cos(interval(i));sin(interval(i))];
    if ~ismember(i, inter)
        hold on
        plot([center(1)-KK*temp(1), center(1)+KK*temp(1)],...
            [center(2)-KK*temp(2), center(2)+KK*temp(2)],'r-','linewidth',2)
    end
end

hold on
scatter(center(1), center(2),100, 'ro','filled')
hold on
quiver(center(1)-KK*coeff(1,1),center(2)-KK*coeff(2,1),...
    2*KK*coeff(1,1), 2*KK*coeff(2,1),'g-','linewidth',2)
hold on
quiver(center(1)-KK*coeff(1,2),center(2)-KK*coeff(2,2),...
    2*KK*coeff(1,2), 2*KK*coeff(2,2),'g-','linewidth',2)
distri=number_focell./length;

p=[];
for i=1:n-1
    temp=linspace(interval(i),interval(i+1),100);
    ptemp=[distri(i)*sin(temp);distri(i)*cos(temp)];
    p=[p;ptemp'];
end
p=p*coeff';

KK1=25;
figure;
scatter(p(:,2),p(:,1),'filled')
for i=1:n-1
     temp=coeff*[cos(interval(i));sin(interval(i))];
    if ~ismember(i, inter)
        hold on
        plot([-KK1*temp(1), +KK1*temp(1)],...
            [+KK1*temp(2), -KK1*temp(2)],'r-','linewidth',2)
    end
end
hold on
scatter(0, 0,100, 'ro','filled')
hold on
quiver(-KK1*coeff(1,1),KK1*coeff(2,1),...
    2*KK1*coeff(1,1), -2*KK1*coeff(2,1),'g-','linewidth',2)
hold on
quiver(-KK1*coeff(1,2),KK1*coeff(2,2),...
    2*KK1*coeff(1,2), -2*KK1*coeff(2,2),'g-','linewidth',2)
axis equal


text:
function filtScan3M = filtImgGabor(scan3M,r, sig, lam, theta, omega)
% AI 03/22/18
% INPUTS
%  r       - final mask will be 2r+1 x 2r+1
%  sig     - standard deviation of Gaussian mask
%  lam     - elongation of Gaussian mask
%  theta   - orientation (in degrees)
%  omega   - [1] wavlength of underlying sine (sould be >=1)

hFilt = filterGabor2d( r, sig, lam, theta, omega, 0 );
filtScan3M = convn(scan3M,hFilt,'same');

end