text:
% This file contains lists of ograns as they are used in the whole-body
% metabolic model.
% Please note that the desired sex has to be specified ('male' or 'female')
%
% Ines Thiele 2016 - 2017

OrgansListShort={ 'Heart'; 'Muscle'; 'Lung';'Skin'; 'Stomach'; 'sIEC'; 'Colon';...
    'Urinarybladder'; 'Retina'; 'Scord';'Brain';'Adipocytes'; ...
    'Liver';'Gall'; 'Kidney';'Pancreas'; 'Spleen';'Agland';'Thyroidgland';'Pthyroidgland';  ...
    'Ovary';'Uterus';...
    'Breast';'Cervix';...
    'Testis';'Prostate';...
    'Bcells';'CD4Tcells'; ...
    'Nkcells';'Monocyte';'Platelet'; 'RBC';%'Salvgland';
    %'Esophagus';'CD8Tcells';'Rectum';
    };
if strcmp(lower(sex),'female')
    OrgansListExt={'Heart'; 'Muscle'; 'Lung';'Skin'; 'Stomach'; 'sIEC'; 'Colon';...
        'Urinarybladder'; 'Retina'; 'Scord';'Brain';'Adipocytes'; ...
        'Liver';'Gall'; 'Kidney';'Pancreas'; 'Spleen';'Agland';'Thyroidgland';'Pthyroidgland';  ...
        'Ovary';'Uterus';...
        'Breast';'Cervix';...
        'Bcells';'CD4Tcells'; ...
        'Nkcells';'Monocyte';'Platelet'; 'RBC';...% 'Salvgland';
        'BBB';'Diet';'SI';'GI';'LI';'BileDuct';'Excretion'
        % 'Rectum';'CD8Tcells';'Esophagus';
        };
    OrgansList={'Heart'; 'Muscle'; 'Lung';'Skin'; 'Stomach'; 'sIEC'; 'Colon';...
        'Urinarybladder'; 'Retina'; 'Scord';'Brain';'Adipocytes'; ...
        'Liver';'Gall'; 'Kidney';'Pancreas'; 'Spleen';'Agland';'Thyroidgland';'Pthyroidgland';  ...
        'Ovary';'Uterus';...
        'Breast';'Cervix';...
        'Bcells';'CD4Tcells';...
        'Nkcells';'Monocyte';'Platelet'; 'RBC';...% 'Salvgland';
        %'Esophagus'; 'CD8Tcells';'Rectum';
        };% without 'BBB';'Diet';'SI';'GI';'LI';'BileDuct';'Excretion'
elseif strcmp(lower(sex),'male')
    OrgansListExt={'Heart'; 'Muscle'; 'Lung';'Skin'; 'Stomach'; 'sIEC'; 'Colon';...
        'Urinarybladder'; 'Retina'; 'Scord';'Brain';'Adipocytes'; ...
        'Liver';'Gall'; 'Kidney';'Pancreas'; 'Spleen';'Agland';'Thyroidgland';'Pthyroidgland';  ...
        'Testis';'Prostate';...
        'Bcells';'CD4Tcells'; ...
        'Nkcells';'Monocyte';'Platelet'; 'RBC';...% 'Salvgland';
        'BBB';'Diet';'SI';'GI';'LI';'BileDuct';'Excretion'
        %'Esophagus';'Rectum';'CD8Tcells';
        };
    OrgansList={'Heart'; 'Muscle'; 'Lung';'Skin'; 'Stomach'; 'sIEC'; 'Colon';...
        'Urinarybladder'; 'Retina'; 'Scord';'Brain';'Adipocytes'; ...
        'Liver';'Gall'; 'Kidney';'Pancreas'; 'Spleen';'Agland';'Thyroidgland';'Pthyroidgland';  ...
        'Testis';'Prostate';...
        'Bcells';'CD4Tcells'; ...
        'Nkcells';'Monocyte';'Platelet'; 'RBC';...% 'Salvgland';
        %'Esophagus';'Rectum'; 'CD8Tcells';
        };
end


text:
function [cMap]=iwarmcold(varargin)

% function [cMap]=iwarmcold(n)
% ------------------------------------------------------------------------
% Creates the colormap data for n levels for the  inverse warm-cold
% colormap.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2018/03/21
%------------------------------------------------------------------------

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

cf=[213 15  37; 238 178 17;]./255;
c=rot90(cf,2);
h=cf;

h=resampleColormap(h,3);
c=resampleColormap(c,3);

cMap=[c; [0 0 0]; h];

[cMap]=resampleColormap(cMap,n);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef OptimizerDualNestedInPrimal < Optimizer_PrimalDual

    properties (GetAccess = public, SetAccess = protected)
        type = 'DualNestedInPrimal'
    end

    properties (Access = private)
        desVarChangedValue
        costIncrease
        constraintProjector
    end

    methods (Access = public)

        function obj = OptimizerDualNestedInPrimal(cParams)
            obj.initOptimizer(cParams);
            obj.createLagrangian();
            obj.createOptimizerUnconstrained(cParams.uncOptimizerSettings)
            obj.createConstraintProjector();
        end

        function solveProblem(obj)
            obj.hasFinished = false;

            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();

            obj.updateOldValues();
            obj.unconstrainedOptimizer.startLineSearch();
            obj.unconstrainedOptimizer.updateConvergenceParams();
            obj.refreshMonitoring();
            obj.printOptimizerVariable();
       %     obj.printHistory();
            obj.nIter = obj.nIter+1;

            obj.designVariable.updateOld();
            obj.computeFeasibleDesignVariable();

            obj.cost.computeFunctionAndGradient();

            obj.updateOldValues();

            obj.unconstrainedOptimizer.updateConvergenceParams();
            obj.unconstrainedOptimizer.updateLineSearch();
            obj.refreshMonitoring();
            obj.printHistory();
            obj.saveDesignVariable();
            obj.printOptimizerVariable();
       %     obj.printHistory();
            obj.nIter = obj.nIter+1;

            %obj.hasFinished = false;
            obj.updateStatus();

            while ~obj.hasFinished

                obj.unconstrainedOptimizer.tryLineSearch();
                while ~obj.hasUnconstraintedOptimizerConverged()
                    obj.restartValues();
                    obj.computeValue();
                    if ~obj.hasUnconstraintedOptimizerConverged()
                        obj.unconstrainedOptimizer.updateLineSearch();
                    end
                end

                obj.cost.computeFunctionAndGradient();

                obj.unconstrainedOptimizer.updateConvergenceParams();

                obj.updateConvergenceStatus();
                obj.updateStatus();

                obj.updateOldValues();

                obj.refreshMonitoring();
                obj.printOptimizerVariable();
                obj.printHistory();
                obj.saveDesignVariable();
                obj.nIter = obj.nIter+1;
            end
            %obj.printOptimizerVariable();
            %obj.printHistory();
            obj.hasConverged = 0;
            obj.printHistoryFinalValues();
        end

        function saveDesignVariable(obj)
            x = obj.designVariable.value;
            mesh = obj.designVariable.mesh.innerMeshOLD;
            path = 'Output/CantileverTetraPerimeterTotal/DesignVariable';
            %save([path,num2str(obj.nIter)],'x','mesh');
        end

        function restartValues(obj)
            obj.designVariable.restart();
            obj.dualVariable.restart();
            obj.cost.restart();
            obj.constraint.restart();
            obj.updateLagrangian();
        end

        function computeFeasibleDesignVariable(obj)
            if obj.isNotFeasible()
                obj.constraint.updateOld();
                obj.updateLagrangian();
                obj.lagrangian.updateOld();

                obj.unconstrainedOptimizer.tryLineSearch();
                %kappa = obj.unconstrainedOptimizer.lineSearch.kappa;
                %obj.unconstrainedOptimizer.lineSearch.kappa = kappa;

                obj.constraintProjector.project();
            end
        end

        function itIsNot = isNotFeasible(obj)
            itIsNot = ~obj.isFeasible();
        end

        function update(obj)
            while ~obj.hasUnconstraintedOptimizerConverged()
                obj.computeValue();
                obj.unconstrainedOptimizer.updateLineSearch();
            end
        end

    end

    methods (Access = protected)

        function createLagrangianSettings(obj)
            cParams.type         = 'Lagrangian';
            cParams.cost         = obj.cost;
            cParams.constraint   = obj.constraint;
            cParams.dualVariable = obj.dualVariable;
            obj.lagrangianSettings = cParams;
        end

        function updateLagrangian(obj)
            obj.lagrangian.computeFunction();
            obj.lagrangian.computeGradient();
        end

        function updateOldValues(obj)
            obj.designVariable.updateOld();
            obj.dualVariable.updateOld();
            obj.cost.updateOld();
            obj.constraint.updateOld();
            obj.updateLagrangian();
            obj.lagrangian.updateOld();
        end

    end

    methods (Access = private)

        function createConstraintProjector(obj)
            cParams.cost           = obj.cost;
            cParams.constraint     = obj.constraint;
            cParams.designVariable = obj.designVariable;
            cParams.dualVariable   = obj.dualVariable;
            cParams.lagrangian     = obj.lagrangian;
            cParams.targetParameters = obj.targetParameters;
            cParams.unconstrainedOptimizer = obj.unconstrainedOptimizer;
            obj.constraintProjector = ConstraintProjector(cParams);
        end

        function computeValue(obj)
            obj.constraintProjector.project();
            obj.cost.computeFunction();
            obj.updateLagrangian();
        end

        function itHas = hasUnconstraintedOptimizerConverged(obj)
            itHas = obj.isStepAcceptable() || obj.isLineSeachTooSmall();
        end

        function itIs = isStepAcceptable(obj)
            incr = obj.lagrangian.computeIncrement();
            costHasDecreased = incr < 0;
            itIs = costHasDecreased;
        end

        function itIs = isLineSeachTooSmall(obj)
            itIs = obj.unconstrainedOptimizer.isLineSearchTooSmall();
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function group = GroupDV(Problem,DV,PV,nPerGroup)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Huangke Chen

   ub  = Problem.upper(DV);
   lb  = Problem.lower(DV);
   dim = length(DV);	% the dim of the distance variables

   fixPV  = (Problem.lower(PV)+Problem.upper(PV))/2;	% the value vector for position variables
   meanDV = (ub+lb)/2;

   f_archiveF = repmat(zeros(dim,dim),1,Problem.M);	% cell(dim, dim);

   fhatDVDecs     = repmat(lb,dim,1) + eye(dim).*repmat(meanDV,dim,1);
   fhatDecs       = [repmat(fixPV,dim,1) fhatDVDecs];
   fhatPopulation = Problem.Evaluation(fhatDecs);
   fhat_archiveF  = fhatPopulation.objs;

   lambdaF = repmat(zeros(dim,dim),1,Problem.M);	% cell(dim, dim);

   p1Dec   = [fixPV lb];
   tempFp1 = Problem.Evaluation(p1Dec);
   fp1F    = tempFp1.objs;

   for i = 1 : dim-1
       fp2 = fhat_archiveF(i,:);
       for j = i+1 : dim

           fp3 = fhat_archiveF(j,:);

           p4      = lb;
           p4(i)   = meanDV(i);	% temp;
           p4(j)   = meanDV(j);	% temp;
           p4Dec   = [fixPV p4];
           tempFp4 = Problem.Evaluation(p4Dec);
           fp4     = tempFp4.objs;

           f_archiveF(i,j:dim:size(f_archiveF,2)) = fp4;

           d1 = fp2 - fp1F;
           d2 = fp4 - fp3;

           lambdaF(i,j:dim:size(lambdaF,2)) = abs(d1-d2);
       end
   end

   % Check for each objective
   bigTheta = false(length(DV));

   for i = 1 : Problem.M

       lambda = lambdaF(:,(i-1)*dim+1:i*dim);
       fp1    = fp1F(i);

       fhat_archive = fhat_archiveF(:,i);

       tempF_archive = f_archiveF(:,(i-1)*dim+1:i*dim);
       f_archive     = tempF_archive + tempF_archive';

       F1 = ones(dim,dim)*fp1;
       F2 = repmat(fhat_archive',dim,1);
       F3 = repmat(fhat_archive,1,dim);
       F4 = f_archive;

       FS   = cat(3,F1,F2,F3,F4);
       Fmax = max(FS,[],3);

       FS       = cat(3,F1+F4,F2+F3);
       Fmax_inf = max(FS,[],3);

       theta = false(dim);

       muM       = eps/2;
       gamma     = @(n)((n.*muM)./(1-n.*muM));
       errlb     = gamma(2)*Fmax_inf;
       errub     = gamma(dim^0.5)*Fmax;
       I2        = lambda >= errub;
       theta(I2) = 1;

       % add, then find not less than 1
       bigTheta = bigTheta + theta;
   end

   bigTheta(bigTheta>0) = true;

   L = size(bigTheta,1);	% number of vertex

   % Breadth-first search
   labels = zeros(1,L);     % all vertex unexplored at the begining
   rts    = [];
   ccc    = 0;              % connected components counter
   while true
       ind = find(labels==0);
       if ~isempty(ind)
           fue  = ind(1);   % first unexplored vertex
           rts  = [rts fue];
           list = [fue];
           ccc  = ccc + 1;
           labels(fue) = ccc;
           while true
               list_new = [];
               for lc = 1 : length(list)
                   p   = list(lc);              % point
                   cp  = find(bigTheta(p,:));	% points connected to p
                   cp1 = cp(labels(cp)==0);     % get only unexplored vertecies
                   labels(cp1) = ccc;
                   list_new    = [list_new cp1];
               end
               list = list_new;
               if isempty(list)
                   break;
               end
           end
       else
           break;
       end
   end

   group_num = max(labels);
   allgroups = cell(1,group_num);
   for i = 1 : group_num
       allgroups{i} = find(labels==i);
   end

   h = @(x)(length(x)==1);
   sizeone = cellfun(h,allgroups);

   seps = allgroups(sizeone);
   seps = cell2mat(seps);

   allgroups(sizeone) = [];
   nonseps            = allgroups;

   group = {};
   for ns = 1 : length(nonseps)     % the non-seperate variables
       group = {group{1:end} DV(nonseps{ns})};
   end

   for g = 1 : nPerGroup : length(seps)
       index = seps(g:min(g+nPerGroup-1,length(seps)));
       group = {group{1: end} DV(index)};
   end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function n = type(d)

%tstoolbox/@description/type
%   return signal type (e.g. 'Correlation function', 'Spectrogram' etc.)
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

n = d.type;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ readerobj ] = open_generic_reader(filename, datasize,...
    datatype, complextype, data_offset, endian, symmetry, bands_ip,...
    blocksize)
%OPEN_GENERIC_READER Generic reader object for reading pixel interleaved
% data from a flat file
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if (nargin<7||isempty(symmetry)), symmetry=[0 0 0]; end; % Assume store azimuth-major order, viewed from above
if (nargin<8||isempty(bands_ip)), bands_ip=1; end; % Assume single band
if (nargin<9||isempty(blocksize)), blocksize=datasize; end; % Assume no blocking

[chipper_function, close_function]=generic_chipper( filename, datasize,...
    datatype, complextype, data_offset, endian, symmetry, bands_ip, blocksize);
if symmetry(3), datasize=datasize(end:-1:1); end;
readerobj = chipfun2readerobj(chipper_function, datasize);
readerobj.close=close_function;

end

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

author: ngageoint


text:
function res = odd(n);
%ODD          Boolean function: n odd?
%
%   res = 1    n odd
%   res = 0    n even
%
%   res = odd(n);
%

% written  10/29/97     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  res = mod(n,2);

repo: OpenMAS


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

%%
clear; close all; clc;

%% Syntax
% |dcmFolder2MATobject_MRE(PathName,MaxVarSize);|

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


text:
function [Q,Z,QAZ,QBZ]=HessenbergTriRed(A,B)
%%HESSENBERGTRIRED Perform a Hessenberg triangular reduction. Given A and
%             B, two nXn matrices, find orthogonal matrices Q and Z such
%             that Q'*A*Z is an upper Hessenberg matrix and Q'*B*Z is an
%             upper triangular matrix.  An upper Hessenberg matrix has zero
%             entries below the first subdisgonal and is thus almost upper
%             triangular. This implementation is only for real matrices.
%
%INPUTS: A, B Two nXn real matrices.
%
%OUTPUTS Q,Z Orthogonal nXn matrices such that Q'*A*Z is upper Hessenberg
%            and Q'*B*Z is lower triangular.
%        QAZ The matrix Q'*A*Z.
%        QBZThe matrix Q'*B*Z
%
%This implements algorithm 7.7.1 in Chapter 7.7.4 of [1], modified to
%obtain the Q and Z matrices explicitly. This is similar to the hess
%function in Matlab.
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(A,1);

%The initial upper triangularization. We are just using qr and not
%Algorithm 5.2.1.
%Make B upper triangular
[Q,B]=qr(B,0);
A=Q'*A;%Adjust A accordingly

Z=eye(n,n);
for j=1:(n-2)
    for i=n:-1:(j+2)
        %Zero A(i,j)
        [c,s]=GivensCS(A(i-1,j),A(i,j));
        A((i-1):i,j:n)=[c, s;-s,c]'*A(i-1:i,j:n);
        B((i-1):i,(i-1):n)=[c, s;-s,c]'*B(i-1:i,i-1:n);
        %Store the transformation
        Q(:,(i-1):i)=Q(:,(i-1):i)*[c,s;-s,c];
        %Zero B(i,j)
        [c,s]=GivensCS(-B(i,i),B(i,i-1));
        B(1:i,i-1:i) = B(1:i,i-1:i)*[c,s;-s,c];
        A(1:n,i-1:i) = A(1:n,i-1:i)*[c,s;-s,c];
        Z(:,(i-1):i) = Z(:,(i-1):i)*[c,s;-s,c];
    end
end

%Save the results.
QAZ=A;
QBZ=B;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function checkTrajectoryLogs(data, ID)
%_________________________________________________________________
%% Documentation
%
% Authors:      Alexander Wischnewski (alexander.wischnewski@tum.de)
% %
% Description:  used to recalculate the checks from the logs
%
% Inputs:
%   File                Log file loaded into data struct
%   ID                  Trajectory ID to be checked

% set parameters
drag_coefficient = 0.725;
roh_air = 1.22;
vehiclemass_kg  = 750;
EmergencyLine_b = 0;

TargetTrajectory.LapCnt = data.debug_slow.debug_slow_tartraj_LapCnt.Data(ID);
TargetTrajectory.TrajCnt = data.debug_slow.debug_slow_tartraj_TrajCnt.Data(ID);
TargetTrajectory.s_loc_m = data.debug_slow.debug_slow_tartraj_s_loc_m.Data(ID, :)';
TargetTrajectory.s_glob_m = data.debug_slow.debug_slow_tartraj_s_glob_m.Data(ID, :)';
TargetTrajectory.x_m = data.debug_slow.debug_slow_tartraj_x_m.Data(ID, :)';
TargetTrajectory.y_m = data.debug_slow.debug_slow_tartraj_y_m.Data(ID, :)';
TargetTrajectory.psi_rad = data.debug_slow.debug_slow_tartraj_psi_rad.Data(ID, :)';
TargetTrajectory.kappa_radpm = data.debug_slow.debug_slow_tartraj_kappa_radpm.Data(ID, :)';
TargetTrajectory.v_mps = data.debug_slow.debug_slow_tartraj_v_mps.Data(ID, :)';
TargetTrajectory.ax_mps2 = data.debug_slow.debug_slow_tartraj_ax_mps2.Data(ID, :)';
TargetTrajectory.banking_rad = data.debug_slow.debug_slow_tartraj_banking_rad.Data(ID, :)';
TargetTrajectory.ax_lim_mps2 = data.debug_slow.debug_slow_tartraj_ax_lim_mps2.Data(ID, :)';
TargetTrajectory.ay_lim_mps2 = data.debug_slow.debug_slow_tartraj_ay_lim_mps2.Data(ID, :)';

traj_ok = checkTrajectory(TargetTrajectory, drag_coefficient, roh_air, vehiclemass_kg, EmergencyLine_b);
if traj_ok
    disp('The trajectory was labelled ok');
else
    disp('The trajectory was labelled not ok');
end

author: TUMFTM


text:
function [mpe, ll] = calc_mpe_bucket(bnet, new_evidence, max_over)
%
% PURPOSE:
%       CALC_MPE Computes the most probable explanation to the network nodes
%       given the evidence.
%
%       [mpe, ll] = calc_mpe(engine, new_evidence, max_over)
%
% INPUT:
%       bnet  - the bayesian network
%       new_evidence - optional, if specified - evidence to be incorporated [cell(1,n)]
%       max_over - optional, if specified determines the variable elimination order [1:n]
%
% OUTPUT:
%       mpe - the MPE assignmet for the net variables (or [] if no satisfying assignment)
%       ll - log assignment probability.
%
% Notes:
% 1. Adapted from '@var_elim_inf_engine\marginal_nodes' for MPE by Ron Zohar, 8/7/01
% 2. Only discrete potentials are supported at this time.
% 3. Complexity: O(nw*) where n is the number of nodes and w* is the induced tree width.
% 4. Implementation based on:
%  - R. Dechter, "Bucket Elimination: A Unifying Framework for Probabilistic Inference",
%                 UA1 96, pp. 211-219.

ns = bnet.node_sizes;
n = length(bnet.dag);
evidence = cell(1,n);
if (nargin<2)
    new_evidence = evidence;
end

onodes = find(~isemptycell(new_evidence));  % observed nodes
hnodes = find(isemptycell(new_evidence));  % hidden nodes
pot_type = determine_pot_type(bnet, onodes);

if pot_type ~= 'd'
  error('only disrete potentials supported at this time')
end

for i=1:n
  fam = family(bnet.dag, i);
  CPT{i} = convert_to_pot(bnet.CPD{bnet.equiv_class(i)}, pot_type, fam(:), evidence);
end

% handle observed nodes: set impossible cases' probability to zero
% rather than prun matrix (this makes backtracking easier)

for ii=onodes
  lIdx = 1:ns(ii);
  lIdx = setdiff(lIdx, new_evidence{ii});

  sCPT=struct(CPT{ii});  % violate object privacy

  sargs = '';
  for jj=1:(length(sCPT.domain)-1)
    sargs = [sargs, ':,'];
  end
  for jj=lIdx
    eval(['sCPT.T(', sargs, num2str(jj), ')=0;']);
  end
  CPT{ii}=dpot(sCPT.domain, sCPT.sizes, sCPT.T);
end

B = cell(1,n);
for b=1:n
  B{b} = mk_initial_pot(pot_type, [], [], [], []);
end

if (nargin<3)
  max_over = (1:n);
end
order = max_over; % no attempt to optimize this

% Initialize the buckets with the CPDs assigned to them
for i=1:n
  b = bucket_num(domain_pot(CPT{i}), order);
  B{b} = multiply_pots(B{b}, CPT{i});
end

% Do backward phase
max_over = max_over(length(max_over):-1:1); % reverse
for i=max_over(1:end-1)
  % max-ing over variable i which occurs in bucket j
  j = bucket_num(i, order);
  rest = mysetdiff(domain_pot(B{j}), i);
  %temp = marginalize_pot_max(B{j}, rest);
  temp = marginalize_pot(B{j}, rest, 1);
  b = bucket_num(domain_pot(temp), order);
  %        fprintf('maxing over bucket %d (var %d), putting result into bucket %d\n', j, i, b);
  sB=struct(B{b});  % violate object privacy
  if ~isempty(sB.domain)
    B{b} = multiply_pots(B{b}, temp);
  else
    B{b} = temp;
  end
end
result = B{1};
marginal = pot_to_marginal(result);
[prob, mpe] = max(marginal.T);

% handle impossible cases
if ~(prob>0)
  mpe = [];
  ll = -inf;
  %warning('evidence has zero probability')
  return
end

ll = log(prob);

% Do forward phase
for ii=2:n
  marginal = pot_to_marginal(B{ii});
  mpeidx = [];
  for jj=order(1:length(mpe))
    assert(ismember(jj, marginal.domain)) %%% bug
    temp = find_equiv_posns(jj, marginal.domain);
    mpeidx = [mpeidx, temp] ;
    if isempty(temp)
      mpeidx = [mpeidx, Inf] ;
    end
  end
  [mpeidxsorted sortedtompe] = sort(mpeidx) ;

  % maximize the matrix obtained from assigning values from previous buckets.
  % this is done by building a string and using eval.

  kk=1;
  sargs = '(';
  for jj=1:length(marginal.domain)
    if (jj~=1)
      sargs = [sargs, ','];
    end
    if (mpeidxsorted(kk)==jj)
      sargs = [sargs, num2str(mpe(sortedtompe(kk)))];
      if (kk<length(mpe))
	kk = kk+1 ;
      end
    else
      sargs = [sargs, ':'];
    end
  end
  sargs = [sargs, ')'] ;
  eval(['[val, loc] = max(marginal.T', sargs, ');'])
  mpe = [mpe loc];
end
[I,J] = sort(order);
mpe = mpe(J);

%%%%%%%%%

function b = bucket_num(domain, order)

b = max(find_equiv_posns(domain, order));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function varargout = rdivide(varargin)
%./   Pointwise right divide of SPHEREFUN objects.
%   F./G if F is a SPHEREFUN and G is a double this returns (1/G)*F
%
%   F./G if F is a double and G is a v this returns F/G, but this doe
%   not work if G becomes numerically close to zero.
%
%   F./G we do not allow F and G to both be SPHEREFUN object.
%
%   F./G is the same as the command rdivide(F,G)
%
% See also SPHEREFUN/LDIVIDE.

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

[varargout{1:nargout}] = rdivide@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%Voice Based Biometric System
%By Ambavi K. Patel.

%finds 20 centroids from equally partitioned mfcc matrix
function cb2=centr(mfcc1)
%partitioning mfcc1 matrix into 20 region and taking centroid from each
[r,c]=size(mfcc1);
dc=floor(c/10);
dr=floor(r/2);
for i=1:10
    for j=1:2
        x=(dr*(j-1))+1:dr*j;
        y=(dc*(i-1))+1:dc*i;
        t=dc+dr;
        cb1(j,i)=sum(sum(mfcc1(x,y)))/t;
    end;
end;
k=1;
for i=1:10
    for j=1:2
        cb2(1,k)=cb1(j,i);
        k=k+1;
    end;
end;

author: Sable

repo: mcbench-benchmarks


text:
function c = cond1est (A)
%COND1EST 1-norm condition estimate.
% Example:
%   c = cond1est(A)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

[m n] = size (A) ;
if (m ~= n | ~isreal (A))                                                   %#ok
    error ('A must be square and real') ;
end
if isempty(A)
    c = 0 ;
    return ;
end
[L,U,P,Q] = lu (A) ;
if (~isempty (find (abs (diag (U)) == 0)))                                  %#ok
    c = Inf ;
else
    c = norm (A,1) * norm1est (L,U,P,Q) ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function line_nco_rule_test01 ( )

%*****************************************************************************80
%
%% LINE_NCO_RULE_TEST01 computes and prints NCO rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2014
%
%  Author:
%
%    John Burkardt
%
  a = -1.0;
  b = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_NCO_RULE_TEST01\n' );
  fprintf ( 1, '  LINE_NCO_RULE computes the Newton-Cotes Open rule\n' );
  fprintf ( 1, '  using N equally spaced points for an interval [A,B].\n' );

  for n = 1 : 12

    [ x, w ] = line_nco_rule ( n, a, b );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Newton-Cotes Open Rule #%d\n', n );
    fprintf ( 1, '   I       X(I)            W(I)\n' );
    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %14.6g  %14.6g\n', i, x(i), w(i) );
    end
    fprintf ( 1, '        Sum(|W)|) =  %14.6g\n', sum ( abs ( w(1:n) ) ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [baseScanNum,filtScanNum,planC] = createFilteredScanForDLS(identifierS,filtS,planC)
% createFilteredScanForDLS.m
%------------------------------------------------------------------------
% INPUTS
% identifierS : Dictionary of scan identifiers.
% filtS       : Dictionary of filter parameters acccepted by processImage.m
% planC
%------------------------------------------------------------------------
% AI 1/5/22

indexS = planC{end};

%% Get scan array from identifiers
%identifierS.filtered = 0;           %Identify base scan
%identifierS.warped = 0;             %Identify base scan
origFlag = 1;                       %Identify base scan
baseScanNum = getScanNumFromIdentifiers(identifierS,planC,origFlag);
scan3M = double(getScanArray(baseScanNum,planC));
scan3M = scan3M - double(...
    planC{indexS.scan}(baseScanNum).scanInfo(1).CTOffset);
mask3M = ones(size(scan3M));

%% Apply filter
imType = fieldnames(filtS.imageType);
imType = imType{1};
filterParS = filtS.imageType.(imType);
outS = processImage(imType, scan3M, mask3M, filterParS);
fieldName = fieldnames(outS);
filtScan3M = outS.(fieldName{1});

%% Add filtered scan to planC
filtScanNum = length(planC{indexS.scan}) + 1;
[xValsV, yValsV, zValsV] = getScanXYZVals(planC{indexS.scan}(baseScanNum));
if yValsV(1) > yValsV(2)
    yValsV = fliplr(yValsV);
end
dx = median(diff(xValsV));
dy = median(diff(yValsV));
sliceThicknessV = diff(zValsV);
scanInfoS.horizontalGridInterval = dy;
scanInfoS.verticalGridInterval = dx;
scanInfoS.coord1OFFirstPoint = xValsV(1);
scanInfoS.coord2OFFirstPoint = yValsV(1);
scanInfoS.zValues = zValsV;
scanInfoS.sliceThickness = [sliceThicknessV,sliceThicknessV(end)];
scanType = ['Filt_scan',num2str(baseScanNum)];
planC = scan2CERR(filtScan3M,scanType,'',scanInfoS,'',planC);
planC{indexS.scan}(filtScanNum).assocBaseScanUID = ...
    planC{indexS.scan}(baseScanNum).scanUID;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [s, cfg] = ft_statfun_actvsblT(cfg, dat, design)

% FT_STATFUN_ACTVSBLT calculates the activation-versus-baseline T-statistic on the
% biological data (the dependent variable), using the information on the independent
% variable (ivar) in design.
%
% Note that it does not make sense to use this test statistic when baseline
% correction was performed by subtracting the mean of the baseline period from the
% whole data (for ERP data) or by dividing by the mean (for TFR data). If baseline
% correction is desired, you should subtract the full baseline and activation period.
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option
%   cfg.statistic = 'ft_statfun_actvsblT'
%
% Configuration options
%   cfg.computestat    = 'yes' or 'no', calculate the statistic (default='yes')
%   cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no')
%   cfg.computeprob    = 'yes' or 'no', calculate the p-values (default='no')
%
% The following options are relevant if cfg.computecritval='yes' and/or
% cfg.computeprob='yes'.
%   cfg.alpha = critical alpha-level of the statistical test (default=0.05)
%   cfg.tail  = -1, 0, or 1, left, two-sided, or right (default=1)
%               cfg.tail in combination with cfg.computecritval='yes'
%               determines whether the critical value is computed at
%               quantile cfg.alpha (with cfg.tail=-1), at quantiles
%               cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at
%               quantile (1-cfg.alpha) (with cfg.tail=1).
%
% Design specification
%   cfg.ivar  = row number of the design that contains the labels of the conditions that must be
%               compared (default=1). The first condition, indicated by 1, corresponds to the
%               activation period and the second, indicated by 2, corresponds to the baseline period.
%   cfg.uvar  = row number of design that contains the labels of the units-of-observation (subjects or trials)
%               (default=2). The labels are assumed to be integers ranging from 1 to
%               the number of units-of-observation.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

% Copyright (C) 2006, Eric Maris
%
% 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$

% set defaults
if ~isfield(cfg, 'computestat'),       cfg.computestat='yes';     end
if ~isfield(cfg, 'computecritval'),    cfg.computecritval='no';   end
if ~isfield(cfg, 'computeprob'),       cfg.computeprob='no';      end
if ~isfield(cfg, 'alpha'),             cfg.alpha=0.05;            end
if ~isfield(cfg, 'tail'),              cfg.tail=1;                end

% perform some checks on the configuration
if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no')
  ft_error('P-values can only be calculated if the test statistics are calculated.');
end

% calculate the number of time samples
switch cfg.dimord
  case 'chan_freq_time'
    nchan = cfg.dim(1);
    nfreq = cfg.dim(2);
    ntime = cfg.dim(3);
    [nsmpls,nrepl] = size(dat);
  otherwise
    ft_error('Inappropriate dimord for the statistics function FT_STATFUN_ACTVSBLT.');
end

sel1 = find(design(cfg.ivar,:)==1);
sel2 = find(design(cfg.ivar,:)==2);
n1  = length(sel1);
n2  = length(sel2);
if (n1+n2)<size(design,2) || (n1~=n2)
  ft_error('Invalid specification of the design array.');
end
nunits = max(design(cfg.uvar,:));
df = nunits - 1;
if nunits<2
  ft_error('The data must contain at least two units (trials or subjects).')
end
if (nunits*2)~=(n1+n2)
  ft_error('Invalid specification of the design array.');
end

if strcmp(cfg.computestat,'yes')
  % compute the statistic
  % calculate the time averages of the activation and the baseline period
  % for all units-of-observation.
  meanreshapeddat=nanmean(reshape(dat,nchan,nfreq,ntime,nrepl),3);
  timeavgdat=repmat(reshape(meanreshapeddat,(nchan*nfreq),nrepl),ntime,1);

  % store the positions of the 1-labels and the 2-labels in a nunits-by-2 array
  poslabelsperunit=zeros(nunits,2);
  poslabel1=find(design(cfg.ivar,:)==1);
  poslabel2=find(design(cfg.ivar,:)==2);
  [dum,i]=sort(design(cfg.uvar,poslabel1),'ascend');
  poslabelsperunit(:,1)=poslabel1(i);
  [dum,i]=sort(design(cfg.uvar,poslabel2),'ascend');
  poslabelsperunit(:,2)=poslabel2(i);

  % calculate the differences between the conditions
  diffmat=dat(:,poslabelsperunit(:,1))-timeavgdat(:,poslabelsperunit(:,2));

  % calculate the dependent samples t-statistics
  avgdiff=nanmean(diffmat,2);
  vardiff=nanvar(diffmat,0,2);
  s.stat=sqrt(nunits)*avgdiff./sqrt(vardiff);
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  s.df      = df;
  if cfg.tail==-1
    s.critval = tinv(cfg.alpha,df);
  elseif  cfg.tail==0
    s.critval = [tinv(cfg.alpha/2,df),tinv(1-cfg.alpha/2,df)];
  elseif cfg.tail==1
    s.critval = tinv(1-cfg.alpha,df);
  end
end

if strcmp(cfg.computeprob,'yes')
  % also compute the p-values
  s.df      = df;
  if cfg.tail==-1
    s.prob = tcdf(s.stat,s.df);
  elseif  cfg.tail==0
    s.prob = 2*tcdf(-abs(s.stat),s.df);
  elseif cfg.tail==1
    s.prob = 1-tcdf(s.stat,s.df);
  end
end

author: spm


text:
%                                     [L.L, L.d, L.skip, L.add] = blkchol(L,X,pars,absd)
% BLKCHOL Fast block sparse Cholesky factorization.
%    The sparse Cholesky factor will be placed in the fields L.L, L.d;
%    the symbolic factorization fields remain unchanged.
%    On input, L should be the symbolic factorization structure of X,
%    as created by SYMBCHOL.
%    Performs LDL' factorization of X(L.perm,L.perm).
%
%   There are important differences with CHOL(X(L.perm,L.perm))':
%
%   -  BLKCHOL uses the supernodal partition L.XSUPER,
%    to use dense linear algebra on dense subblocks.
%    This explains the performance benefit of BLKCHOL over CHOL.
%
%   -  BLKCHOL never fails.
%
%   - To solve "X*y = b", use
%    >> [L.L,L.d,L.skip,L.add] = blkchol(symbchol(X),X);
%    >> L.d(find(L.skip)) = inf;
%    >> y = sparbwslv(L, sparfwslv(L,b) ./ L.d);
%
% For numerical reasons, "blkchol" may skip unstable pivots,
% or add on the diagonal. Such pivots are then listed in L.{skip,add}.
%
% See also symbchol, sparfwslv, sparbwslv, [symbfact, symmmd, chol].

function [LL, Ld, Lskip, Ladd] = blkchol(L,X,pars,absd) %#ok

 %
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
 %

error('Build the SeDuMi binaries by typing make at the command prompt.');

repo: REST


text:
function snr = get_snr(data)
% Data is a matrix whos columns index voxels, and rows index subjects (or trials, etc.)
%
% :Usage:
% ::
%
%     snr = get_snr(data)
%
% ..
%    Tor Wager
% ..

mystd = nanstd(data);
mystd(mystd == 0) = NaN;
snr = nanmean(data) ./ mystd;

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%*****************************************************************
%% sqlptermcode.m: explains the termination code in sqlp.m
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function termcode = sqlptermcode;

   fprintf('\n   3: norm(X) or norm(Z) diverging');
   fprintf('\n   2: dual   problem is suspected to be infeasible')
   fprintf('\n   1: primal problem is suspected to be infeasible')
   fprintf('\n   0: max(relative gap,infeasibility) < gaptol');
   fprintf('\n  -1: relative gap < infeasibility');
   fprintf('\n  -2: lack of progress in predictor or corrector');
   fprintf('\n  -3: X or Z not positive definite');
   fprintf('\n  -4: difficulty in computing predictor or corrector direction');
   fprintf('\n  -5: progress in relative gap or infeasibility is bad');
   fprintf('\n  -6: maximum number of iterations reached');
   fprintf('\n  -7: primal infeasibility has deteriorated too much');
   fprintf('\n  -8: progress in relative gap has deteriorated');
   fprintf('\n  -9: lack of progress in infeasibility');
   fprintf('\n')
%%*************************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function out = isfinite(f)
%ISFINITE   Test if a DELTAFUN is bounded.
%   ISFINITE(F) returns FALSE if F has any non trivial delta functions or
%   if the smooth part is infinite. It is TRUE otherwise.
%
% See also ISINF, ISNAN.

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

% Check the smooth part:
if ( ~isfinite(f.funPart) )
    out = 0;
    return
end

% Smooth part is finite, check the distributional part:
if ( isempty(f.deltaLoc) || isempty(f.deltaMag) )
    out = 1;
    return
end

% Get the tolerance:
pref = chebfunpref();
deltaTol = pref.deltaPrefs.deltaTol;

if ( max(abs(f.deltaMag)) < deltaTol )
    out = 1;
else
    out = 0;
end

end


text:
function f = f08_f0 ( n, x, y )

%*****************************************************************************80
%
%% F08_F0 returns the value of function 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real F(N,1), the function values.
%
  t1(1:n,1) = 5.0 - 10.0 * x(1:n,1);
  t2(1:n,1) = 5.0 - 10.0 * y(1:n,1);
  t3(1:n,1) = exp ( - 0.5 * t1(1:n,1) .* t1(1:n,1) );
  t4(1:n,1) = exp ( - 0.5 * t2(1:n,1) .* t2(1:n,1) );
  f(1:n,1) = t3(1:n,1) + 0.75 * t4(1:n,1) .* ( 1.0 + t3(1:n,1) );

  return
end


text:
function a = mid(a)
%MID          Hessian midpoint
%
%  r = mid(a)
%

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

  if isa(a.x,'intval')
    a.x = mid(a.x);
    a.dx = mid(a.dx);
    a.hx = mid(a.hx);
  end


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

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

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

predictions =  X*theta;
sqerrors = (predictions - y).^2;
J = 1/(2*m)* sum(sqerrors);

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

end

author: vugsus

repo: coursera-machine-learning


text:
function bf = const_bf(PATH)

% Path (x,y,z) by nodes
m = size(PATH,1);

bf = zeros(10*(m-1),size(PATH,2));
bf(1,:) = PATH(1,:);
bf(end-4,:) = PATH(end,:);
idx = 1;
for i=1:m-2
    bf(5*(idx)+1,:) = PATH(i+1,:);
    idx = idx+1;
    bf(5*(idx)+1,:) = PATH(i+1,:);
    idx = idx+1;
end
idx = [];
for i=1:m-2
    idx = [idx, 10*i+2:10*i+5];
end
bf(idx,:) = [];

author: yorgoon

repo: minimum-snap-geometric-control


text:
% 2012-fall by Greg Handy, based on cbct_back.m
% 2013-04-04 refined by Rebecca Malinas
% 2013-04-07 refined by Jeff Fessler

function img = cbct_back_mat_par(interpRect,rectGrid,proj, ns, nt, na, ...
	ds, dt, offset_s, offset_t, offset_source, ...
	dsd, dso, dfs, orbit, orbit_start, ...
	source_zs, ...
	mask, nz, dx, dy, dz, offset_xyz, ia_skip, scale_dang, extrapolate_t)

if any(source_zs ~= 0)
	warn('helix not yet tested')
end

[nx ny] = size(mask);
betas = deg2rad(orbit_start + orbit * [0:na-1] / na); % [na] source angles

% precompute as much as possible
wx = (nx-1)/2;
wy = (ny-1)/2;
wz = (nz-1)/2;
[xc yc] = ndgrid(([0:nx-1] - wx) * dx, ([0:ny-1] - wy) * dy);
zc = ([0:nz-1] - wz) * dz;

clear wx wy wz rr smax rmax

xc = xc(mask); % [np] pixels within mask
yc = yc(mask);

ws = (ns+1)/2 + offset_s; % trick: +1 because matlab starts from 1
wt = (nt+1)/2 + offset_t; % offset_t should be zero for cone-par and tent geometry

% loop over slices
img = zeros([size(mask) nz]);
sdim = [ns+1 nt]; % trick: extra zeros saves indexing in loop
proj1 = zeros(sdim);
ticker reset
for iz=1:nz

	ia_min = 1;
	ia_max = na;

	% loop over each projection angle
	img2 = 0;
	for ia=ia_min:ia_skip:ia_max
		ticker(mfilename, [iz ia], [nz na])
		beta = betas(ia);

		x_beta = +xc * cos(beta) + yc * sin(beta);
		y_betas = (-xc * sin(beta) + yc * cos(beta));

		if interpRect
			mag = sqrt(dso^2-(x_beta).^2) ./ (sqrt(dso^2-(x_beta).^2)+y_betas); %T-FDK?
		else
			mag = dsd ./ (sqrt(dso^2-(x_beta).^2)+y_betas); %P-FDK?
		end

		sprime = x_beta;

		tprime = mag .* (zc(iz)-source_zs(ia));

			bs = sprime / ds + ws;

		if interpRect
			newdt = rectGrid(2)-rectGrid(1);
			bt = tprime / newdt + wt;
		else
			bt = tprime / dt + wt;
		end

		bt = max(bt, 1);
		bt = min(bt, nt);

		% bi-linear interpolation:
		is = floor(bs); % left bin
		it = floor(bt);

		itbad = (it==nt);
		it(itbad)=nt-1;

		isbad = (bs < 1) | (bs > ns);
		is(isbad) = ns;

		wr = bs - is;	% left weight
		wr(isbad)=1;

		wl = 1 - wr;	% right weight
		wu = bt - it;	% upper weight
		wd = 1 - wu;	% lower weight

%		ibad = (is < 0) | (is > ns) | (it < 0) | (it > nt);
%		is(ibad) = ns+1; % trick! point at harmless zeros
%		it(ibad) = nt+1;

 %		proj1(1+[1:ns],1+[1:nt]) = proj(:,:,ia); % trick: left side

%		p1 =	wl .* proj1(sub2ind(sdim, is+1,it+1)) + ...
%			wr .* proj1(sub2ind(sdim, is+2,it+1));
%		p2 =	wl .* proj1(sub2ind(sdim, is+1,it+2)) + ...
%			wr .* proj1(sub2ind(sdim, is+2,it+2));

									proj1(1:ns,1:nt) = proj(:,:,ia);%now only have extra row of zeros at ns+1 % trick: left side

		p1 =	wl .* proj1(sub2ind(sdim, is, it)) + ...
			wr .* proj1(sub2ind(sdim, is+1,it));

		p2 =	wl .* proj1(sub2ind(sdim, is, it+1)) + ...
			wr .* proj1(sub2ind(sdim, is+1,it+1));
		p0 = wd .* p1 + wu .* p2; % vertical interpolation

		% no backprojection weighting needed
		img2 = img2 + p0;
	end % ia

	img(:,:,iz) = embed(img2, mask);
end % iz

if scale_dang % final "\der angle" scale:
	img = (0.5 * deg2rad(abs(orbit)) / (na/ia_skip)) * img;
end

end % cbct_back_mat_par()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%ISTEREO Stereo matching
%
% D = ISTEREO(LEFT, RIGHT, RANGE, H, OPTIONS) is a disparity image computed
% from the epipolar aligned stereo pair: the left image LEFT (HxW) and the
% right image RIGHT (HxW).  D (HxW) is the disparity and the value at each
% pixel is the horizontal shift of the corresponding pixel in IML as observed
% in IMR. That is, the disparity d=D(v,u) means that the pixel at RIGHT(v,u-d)
% is the same world point as the pixel at LEFT(v,u).
%
% RANGE is the disparity search range, which can be a scalar for disparities in
% the range 0 to RANGE, or a 2-vector [DMIN DMAX] for searches in the range
% DMIN to DMAX.
%
% H is the half size of the matching window, which can be a scalar for NxN or a
% 2-vector [N,M] for an NxM window.
%
% [D,SIM] = ISTEREO(LEFT, RIGHT, RANGE, H, OPTIONS) as above but returns SIM
% which is the same size as D and the elements are the peak matching score
% for the corresponding elements of D.  For the default matching metric ZNCC
% this varies between -1 (very bad) to +1 (perfect).
%
% [D,SIM,DSI] = ISTEREO(LEFT, RIGHT, RANGE, H, OPTIONS) as above but returns DSI
% which is the disparity space image (HxWxN) where N=DMAX-DMIN+1. The I'th
% plane is the similarity of IML to IMR shifted to the left by DMIN+I-1.
%
% [D,SIM,P] = ISTEREO(LEFT, RIGHT, RANGE, H, OPTIONS) if the 'interp' option is
% given then disparity is estimated to sub-pixel precision using quadratic
% interpolation.  In this case D is the interpolated disparity and P is
% a structure with elements A, B, dx.  The interpolation polynomial is
% s = Ad^2 + Bd + C where s is the similarity score and d is disparity relative
% to the integer disparity at which s is maximum.  P.A and P.B are matrices the
% same size as D whose elements are the per pixel values of the interpolation
% polynomial coefficients.  P.dx is the peak of the polynomial with respect
% to the integer disparity at which s is maximum (in the range -0.5 to +0.5).
%
% Options::
% 'metric',M   string that specifies the similarity metric to use which is
%              one of 'zncc' (default), 'ncc', 'ssd' or 'sad'.
% 'interp'     enable subpixel interpolation and D contains non-integer
%              values (default false)
% 'vshift',V   move the right image V pixels vertically with respect to left.
%
% Example::
%
% Load the left and right images
%         L = iread('rocks2-l.png', 'reduce', 2);
%         R = iread('rocks2-r.png', 'reduce', 2);
% then compute stereo disparity and display it
%         d = istereo(L, R, [40, 90], 3);
%         idisp(d);
%
% References::
%  - Robotics, Vision & Control, Section 14.3,
%    P. Corke, Springer 2011.
%
% Notes::
% - Images must be greyscale.
% - Disparity values pixels within a half-window dimension (H) of the edges
%   will not be valid and are set to NaN.
% - The C term of the interpolation polynomial is not computed or returned.
% - The A term is high where the disparity function has a sharp peak.
% - Disparity and similarity score can be obtained from the disparity space
%   image by [SIM,D] = max(DSI, [], 3)
%
%
% See also IRECTIFY, STDISP.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function [disp,sim, o3] = istereo(L, R, drange, h, varargin)

% TODO: score cube is float, can return it

    opt.metric = 'zncc';
    opt.interp = false;
    opt.vshift = 0;

    opt = tb_optparse(opt, varargin);

    % ensure images are greyscale
    L = imono(L);
    R = imono(R);

    opt.vshift = round(opt.vshift);
    if opt.vshift ~= 0
        if opt.vshift > 0
            L = L(1:end-opt.vshift,:);
            R = R(opt.vshift:end,:);
        else
            vshift = -vshift;
            L = L(opt.vshift:end,:);
            R = R(1:end-opt.vshift,:);
        end
    end

    % compute the score cube, 3rd dimension is disparity
    DSI = stereo_match(L, R, 2*h+1, drange(1:2), opt.metric);

    % best value along disparity dimension is the peak
    %   s best score
    %   d disparity at which it occurs
    %
    % both s and d are matrices same size as L and R.
    if strcmp(opt.metric, 'sad') | strcmp(opt.metric, 'ssd')
        [s,d] = min(DSI, [], 3);
    else
        [s,d] = max(DSI, [], 3);
    end

    d(isnan(s)) = NaN;

    if opt.interp
        % interpolated result required

        % get number of pixels and disparity range
        npix = prod(size(L));

        if length(drange) == 1,
            ndisp = drange + 1;
        else
            dmin = min(drange);
            dmax = max(drange);
            ndisp = dmax - dmin + 1;
        end

        % find all disparities that are not at either end of the range, we need
        % a point on either side to interpolate them
        valid = (d>1) & (d<ndisp);
        valid = valid(:);

        % make a vector of consecutive pixel indices (1 to width*height)
        ci = [1:npix]';
        % turn disparities into a column vector
        dcol = d(:);

        % remove all entries that are not valid
        ci(~valid) = [];
        dcol(~valid) = [];

        % both ci and dcol have the same number of entries

        % for every valid pixel and disparity, find the index into the 3D score
        % array.  We cheat and consider that array WxHxD as a 2D array (WxH)xD
        %
        % We compute the indices for the best score and one each side of it
        k_m = sub2ind([npix ndisp], ci, dcol-1);
        k_0 = sub2ind([npix ndisp], ci, dcol);
        k_p = sub2ind([npix ndisp], ci, dcol+1);

        % initialize matrices (size of L and R) to hold the the best score
        % and the one each side of it
        y_m = ones(size(L))*NaN;
        y_0 = ones(size(L))*NaN;
        y_p = ones(size(L))*NaN;

        % now copy over the valid scores into these arrays.  What doesnt
        % get copies is a NaN
        y_m(ci) = DSI(k_m);
        y_0(ci) = DSI(k_0);
        y_p(ci) = DSI(k_p);

        % figure the coefficients of the peak fitting parabola:
        %    y = Ax^2 + Bx + C
        % Each coefficient is a matrix same size as (L and R)
        % We don't need to compute C
        A = 0.5*y_m - y_0 + 0.5*y_p;
        B = -0.5*y_m + 0.5*y_p;

        % now the position of the peak is given by -B/2A
        dx = -B ./ (2*A);

        % and we add this fractional part to the integer value obtained
        % from the max/min function
        d = d + dx;

   end
   d = d + drange(1)-1;

    if nargout > 0,
        disp = d;
    end
    if nargout > 1,
        sim = s;
    end

    if nargout > 2
        if opt.interp
            o3.A = A;
            o3.B = B;
            o3.dx = dx;
        else
            o3 = DSI;
        end
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [ a, seed ] = markov_random ( n, key )

%*****************************************************************************80
%
%% MARKOV_RANDOM returns a random Markov matrix.
%
%  Discussion:
%
%    A Markov matrix, also called a "stochastic" matrix, is distinguished
%    by two properties:
%
%    * All matrix entries are nonnegative;
%    * The sum of the entries in each row is 1.
%
%    A "transition matrix" is the transpose of a Markov matrix, and
%    has column sums equal to 1.
%
%  Example:
%
%    N = 4
%
%    1/10  2/10  3/10  4/10
%    1     0     0     0
%    5/10  2/10  3/10  0
%    2/10  2/10  2/10  4/10
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    0 <= A(I,J) <= 1.0 for every I and J.
%
%    The sum of the entries in each row of A is 1.
%
%    Because it has a constant row sum of 1,
%    A has an eigenvalue of 1, and
%    a (right) eigenvector of ( 1, 1, 1, ..., 1 ).
%
%    All the eigenvalues of A have modulus no greater than 1.
%
%    The eigenvalue 1 lies on the boundary of all the Gerschgorin rowsum disks.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, integer KEY, a positive value that selects the data.
%
%    Output, real A(N,N), the matrix.
%
  seed = key;

  [ a, seed ] = r8mat_uniform_01 ( n, n, seed );

  for i = 1 : n

    row_sum = sum ( a(i,1:n) );

    a(i,1:n) = a(i,1:n) / row_sum;

  end

  return
end


text:
function value = i4vec_any_nonzero ( n, a )

%*****************************************************************************80
%
%% I4VEC_ANY_NONZERO: ( any A nonzero ) for I4VEC's.
%
%  Discussion:
%
%    An I4VEC is a vector of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer A(N), the vector.
%
%    Output, logical VALUE is TRUE if any entry is nonzero.
%
  value = any ( a(1:n) ~= 0 );

  return
end


text:
function [ r, center ] = tetrahedron_circumsphere_3d ( tetra )

%*****************************************************************************80
%
%% TETRAHEDRON_CIRCUMSPHERE_3D computes the circumsphere of a tetrahedron in 3D.
%
%  Discussion:
%
%    The circumsphere, or circumscribed sphere, of a tetrahedron is the sphere that
%    passes through the four vertices.  The circumsphere is not necessarily
%    the smallest sphere that contains the tetrahedron.
%
%    Surprisingly, the diameter of the sphere can be found by solving
%    a 3 by 3 linear system.  This is because the vectors P2 - P1,
%    P3 - P1 and P4 - P1 are secants of the sphere, and each forms a
%    right triangle with the diameter through P1.  Hence, the dot product of
%    P2 - P1 with that diameter is equal to the square of the length
%    of P2 - P1, and similarly for P3 - P1 and P4 - P1.  This determines
%    the diameter vector originating at P1, and hence the radius and
%    center.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real TETRA(3,4) the tetrahedron vertices.
%
%    Output, real R, CENTER(3), the center of the
%    circumscribed sphere, and its radius.  If the linear system is
%    singular, then R = -1, CENTER(1:3) = 0.
%
  dim_num = 3;
  nrhs = 1;
%
%  Set up the linear system.
%
  a(1:dim_num,1:3) = ( tetra(1:dim_num,2:4) )';

  for j = 1 : dim_num
    a(1:dim_num,j) = a(1:dim_num,j) - tetra(j,1);
  end

  for i = 1 : dim_num
    a(i,4) = sum ( a(i,1:3).^2 );
  end
%
%  Solve the linear system.
%
  [ a, info ] = r8mat_solve ( dim_num, nrhs, a );
%
%  If the system was singular, return a consolation prize.
%
  if ( info ~= 0 )
    r = -1.0;
    center(1:dim_num) = 0.0;
    return
  end
%
%  Compute R, X, Y, Z.
%
  r = 0.5 * sqrt ( sum ( a(1:dim_num,dim_num+1).^2 ) );

  center(1:dim_num) = tetra(1:dim_num,1) + 0.5 * a(1:dim_num,dim_num+1);

  return
end


text:
function varargout = process_evt_extended( varargin )
% PROCESS_EVT_EXTENDED: Convert simple events to extended events.
%
% USAGE:  OutputFiles = process_evt_extended('Run', sProcess, sInput)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Convert to extended event';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 62;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/EventMarkers#Other_menus';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Event name
    sProcess.options.eventname.Comment = 'Event names: ';
    sProcess.options.eventname.Type    = 'text';
    sProcess.options.eventname.Value   = '';
    % Time offset
    sProcess.options.timewindow.Comment = 'Time window around event:';
    sProcess.options.timewindow.Type    = 'range';
    sProcess.options.timewindow.Value   = {[-0.2, 0.2], 's', []};
end

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

%% ===== RUN =====
function OutputFile = Run(sProcess, sInput) %#ok<DEFNU>
    % Return all the input files
    OutputFile = {sInput.FileName};

    % ===== GET OPTIONS =====
    % Time offset
    TimeWindow = sProcess.options.timewindow.Value{1};
    % Event names
    EvtNames = strtrim(str_split(sProcess.options.eventname.Value, ',;'));
    if isempty(EvtNames)
        bst_report('Error', sProcess, [], 'No events selected.');
        return;
    end

    % ===== LOAD FILE =====
    % Get file descriptor
    isRaw = strcmpi(sInput.FileType, 'raw');
    % Load the raw file descriptor
    if isRaw
        DataMat = in_bst_data(sInput.FileName, 'F');
        sEvents = DataMat.F.events;
        sFreq = DataMat.F.prop.sfreq;
    else
        DataMat = in_bst_data(sInput.FileName, 'Events', 'Time');
        sEvents = DataMat.Events;
        sFreq = 1 ./ (DataMat.Time(2) - DataMat.Time(1));
    end
    % If no markers are present in this file
    if isempty(sEvents)
        bst_report('Error', sProcess, sInput, 'This file does not contain any event.');
        return;
    end
    % Find event names
    iEvtList = [];
    for i = 1:length(EvtNames)
        iEvt = find(strcmpi(EvtNames{i}, {sEvents.label}));
        if isempty(iEvt)
            bst_report('Warning', sProcess, sInput, 'This file does not contain any event.');
        elseif isempty(sEvents(iEvt).times)
            bst_report('Warning', sProcess, sInput, ['Event category "' sEvents(iEvt).label '" is empty.']);
        elseif (size(sEvents(iEvt).times,1) > 1)
            bst_report('Warning', sProcess, sInput, ['Event category "' sEvents(iEvt).label '" already contains extended events.']);
        else
            iEvtList(end+1) = iEvt;
        end
    end
    % No events to process
    if isempty(iEvtList)
        bst_report('Error', sProcess, sInput, 'No events to process.');
        return;
    end

    % ===== PROCESS EVENTS =====
    for i = 1:length(iEvtList)
        sEvents(iEvtList(i)).times = round([...
            sEvents(iEvtList(i)).times + TimeWindow(1);
            sEvents(iEvtList(i)).times + TimeWindow(2)] .* sFreq) ./ sFreq;
    end

    % ===== SAVE RESULT =====
    % Report results
    if isRaw
        DataMat.F.events = sEvents;
    else
        DataMat.Events = sEvents;
    end
    % Only save changes if something was change
    bst_save(file_fullpath(sInput.FileName), DataMat, 'v6', 1);
end

author: brainstorm-tools


text:
function doinvers_michael(mycat)
    % doinvers_michael calculates orientation of the stress tensor based on Gephard's algorithm.
    % stress tensor orientation. The actual calculation is done using a call to a fortran program.
    %
    % Stefan Wiemer 03/96

    % todo: HOW DOES THIS DIFFER FRON DOINVERSE?

    global mif1 mif2
    global tmpi cumu2
    % Save path to current folder
    sPath = pwd;
    ZG = ZmapGlobal.Data;
    fs = filesep;

    report_this_filefun();

    % Number of bootstraps
    rep = 2000;
    %if isunix ~= 1
    % errordlg('Misfit calculation only implemented for UNIX version! ');
    % return
    %end

    % making assumption that columns 10-12 were dip, dip idirection, etc.
    if isempty(mycat.Dip) || isempty(mycat.DipDirection) || isempty(mycat.Rake)
        errordlg('You do not have the required fields to calculate misfit');
        return
    end

    hodis = fullfile(hodi, 'external');

    tmpi = [ZG.mycat.Dip , ZG.mycat.DipDirection , ZG.mycat.Rake];

    cd(hodis);

    fid = fopen('data2','w');
    str = ['Inversion data'];
    str = str';

    fprintf(fid,'%s  \n',str');
    fprintf(fid,'%7.3f  %7.3f  %7.3f \n',tmpi');

    fclose(fid);

    delete data2.slboot Xtemp.slboot
    disp('Now doing inversion ... ')
    watchon
    helpdlg('The inversion is running right now ... it will take a few seconds ... please wait until results appear ');

    % slfast calculates the best solution for the stress tensor according to
    % Michael(1987): creates data2.slboot for bootslickw
    %[stat, res] = unix(['.' fs 'slfast data2 ']);

    slick_program = append_system_specific_postfix('slick');
    slfast_program = append_system_specific_postfix('slfast');
    bootslickw_program = append_system_specific_postfix('bootslickw');

    [stat, res] =  system([fullfile(ZG.hodi, slfast_program), ' data2'])

    if strcmp(res,'') == 0
        helpdlg('It seems that the inversion did not run, because the command slick_* data2 could not be executed. Is the executebale slick in the directory extrenal? ','error inverting');
        return
    end
    % slick calculates the best solution for the stress tensor according to
    % Michael(1987): creates data2.oput
    system([fullfile(ZG.hodi, slick_program), ' data2'])

    % Get data from data2.oput
    sFilename = 'data2.oput';
    [fBeta, fStdBeta, fTauFit, fAvgTau, fStdTau] = import_slickoput(sFilename);

    load data2.slboot
    % d0: Best solution for stress tensor
    d0 = data2;

    % bootslickw resamples the data and uses slfast to calculate the best fitting stress tensor
    system(sprintf('"%s" data2 %d 0.5', fullfile('.',bootslickw_program), rep));

    disp(' Done !  ')
    watchoff

    load Xtemp.slboot
    l = 2:2:2*rep;
    % d1: Bootstrap stress tensor solutions
    d1 = Xtemp(l,:);

    % Back to the original folder
    cd(sPath);

    % Now plot the results
    figure
    wulff
    set(gca,'NextPlot','add')

    %gridsize
    radius=0.05;

    % Plotting the bootstrap solutions of S1
    X = [d1(:,3) d1(:,2) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin

    colormap(hot);
    brighten(0.5);
    densiter=denserfocalv2(rho,theta,radius);
    [xvec, yvec] =meshgrid(-1:0.01:1,-1:0.01:1);
    zvec=griddata(densiter(:,1),densiter(:,2),densiter(:,3),xvec,yvec);
    p11=contourf(xvec,yvec,zvec);

    %calculate 95% confidence interval
    bpoints(:,1)= rho .* cos(theta); %xval bootstrap
    bpoints(:,2)= rho .* sin(theta); %yval bootstrap

    X = [d0(2,3) d0(2,2) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin
    solpoint(:,1)= rho .* cos(theta); %xval solution
    solpoint(:,2)= rho .* sin(theta); %yval solution

    cline = contour_dense(bpoints, solpoint, 0.95, 36);
    plot(cline(:,1),cline(:,2),'k', 'LineWidth',3,'LineStyle','--',...
        'DisplayName','95%');

    % Plotting the bootstrap solutions of S2
    X = [d1(:,5) d1(:,4) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin

    densiter=denserfocalv2(rho,theta,radius);
    [xvec, yvec] =meshgrid(-1:0.01:1,-1:0.01:1);
    zvec=griddata(densiter(:,1),densiter(:,2),densiter(:,3),xvec,yvec);
    p12=contourf(xvec,yvec,zvec);

    %calculate 95% confidence interval
    bpoints(:,1)= rho .* cos(theta); %xval bootstrap
    bpoints(:,2)= rho .* sin(theta); %yval bootstrap

    X = [d0(2,5) d0(2,4) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin
    solpoint(:,1)= rho .* cos(theta); %xval solution
    solpoint(:,2)= rho .* sin(theta); %yval solution

    cline = contour_dense(bpoints, solpoint, 0.95, 36);
    plot(cline(:,1),cline(:,2),'k', 'LineWidth',3,'LineStyle','--',...
        'DisplayName','95%');

    % Plotting the bootstrap solutions of S3
    X = [d1(:,7) d1(:,6) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin

    densiter=denserfocalv2(rho,theta,radius);
    [xvec, yvec] =meshgrid(-1:0.01:1,-1:0.01:1);
    zvec=griddata(densiter(:,1),densiter(:,2),densiter(:,3),xvec,yvec);
    p13=contourf(xvec,yvec,zvec);

    %calculate 95% confidence interval
    bpoints(:,1)= rho .* cos(theta); %xval bootstrap
    bpoints(:,2)= rho .* sin(theta); %yval bootstrap

    X = [d0(2,7) d0(2,6) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin
    solpoint(:,1)= rho .* cos(theta); %xval solution
    solpoint(:,2)= rho .* sin(theta); %yval solution

    cline = contour_dense(bpoints, solpoint, 0.95, 36);
    plot(cline(:,1),cline(:,2),'k', 'LineWidth',3,'LineStyle','--',...
        'DisplayName','95%');

    % Plot the best solution
    % Replaced d1 with d0, JW
    X = [d0(2,3) d0(2,2) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin
    xp = rho .* cos(theta);
    yp = rho .* sin(theta);
    pl1a = plot(xp,yp,'ks'); %s
    set(pl1a,'LineWidth',2,'MarkerSize',12,'Markerfacecolor','w')
    set(gca,'NextPlot','add')

    X = [d0(2,5) d0(2,4) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin
    xp = rho .* cos(theta);
    yp = rho .* sin(theta);
    pl2a = plot(xp,yp,'k^'); %^
    set(pl2a,'LineWidth',2,'MarkerSize',12,'Markerfacecolor','w')

    X = [d0(2,7) d0(2,6) ];
    theta = pi*(90-X(:,2))/180;      %az converted to MATLAB angle
    rho = tan(pi*(90-X(:,1))/360);   %projected distance from origin
    xp = rho .* cos(theta);
    yp = rho .* sin(theta);
    pl3a = plot(xp,yp,'ok'); %o
    set(pl3a,'LineWidth',2,'MarkerSize',12,'Markerfacecolor','w')
    set(gcf,'color','w');

    %doing the legend
    le = legend([pl1a pl2a pl3a],'S1','S2','S3');
    set(le,'pos',[0.1 0.8 0.15 0.1]);
    set(le,'Xcolor','w','ycolor','w','box','off');

    % Compute standard deviation of phi
    fStdPhi = std(d1(:,1),1,'omitnan');

    axes('pos',[0 0 1 1 ]);
    axis off
    text(0.01,0.22,['Variance: ' num2str(d0(1,1),2) ]);
    text(0.01,0.18,['Phi: ' num2str(d0(2,1),2) ' \pm ' num2str(fStdPhi)]);
    text(0.01,0.14,['S1: trend: ' num2str(d0(2,2),4) '; plunge: '  num2str(d0(2,3),4) ]);
    text(0.01,0.1,['S2: trend: ' num2str(d0(2,4),4) '; plunge: '  num2str(d0(2,5),4) ]);
    text(0.01,0.06,['S3: trend: ' num2str(d0(2,6),4) '; plunge: '  num2str(d0(2,7),4) ]);
    text(0.78,0.22,['Tau: ' num2str(fAvgTau,4) ' \pm ' num2str(fStdTau,4)]);
    text(0.78,0.18,['Ratio: ' num2str(fTauFit,4)]);
    text(0.78,0.14,['Beta: ' num2str(fBeta,4) ' \pm ' num2str(fStdBeta,4)]);

    % Determine the faulting style based on Zoback, 1992
    ste = [d0(2,3) d0(2,2)+180 d0(2,5) d0(2,4)+180 d0(2,7) d0(2,6)+180];

    faulting_type = faultType(ste(1),ste(3),ste(5));
    if strcmpi(faulting_type,'Thrust to Strike Slip')
        l4 = pl; %?
    end

    text(0.01,0.02,['Faulting style: ' faulting_type]);
    % Link to World stress map
    uicontrol('Units','normal',...
        'Position',[.4 .0 .1 .04],'String','Info ',...
        'callback',@callbackfun_001);

    % Plot CDF and histogram of phi-values
    figure_w_normalized_uicontrolunits('Name','Phi - ratio of relative stress magnitude','Position',[700   700   472   232])
    % CDF
    subplot(1,2,1)
    vDistribution = d1(:,1);
    % Remove NaN-values
    vSel = ~isnan(vDistribution);
    vPlotDist = vDistribution(vSel,:);
    % Plot the cumulative density function
    nLen = length(vPlotDist);
    vIndices = [1:nLen]/nLen;
    vDist = sort(vPlotDist);
    plot(vDist, vIndices,'Linewidth',1.5,'Color',[0 0 0]);
    set(gca,'NextPlot','add')
    xlabel('\phi','FontSize',12,'Fontweight','bold')
    ylabel('Fraction of cases','FontSize',12,'Fontweight','bold')
    set(gca,'Box','on','Linewidth',1.5,'Xlim',[0 1])
    plot([d0(2,1) d0(2,1)],[0 1],'Color',[0.8 0 0],'Linewidth',1.5)
    legend(' \phi CDF bootstrap','\phi org. data')
    % Histogram
    subplot(1,2,2)
    histogram(vDistribution,0:0.05:1)
    set(gca,'Xlim',[0 1]);
    xlabel('\phi','FontSize',12,'Fontweight','bold')
    ylabel('Frequency','FontSize',12,'Fontweight','bold')

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        web http://www.world-stress-map.org/;
    end

    function ftype = faultType(s1_east, s2_west, s3_north)
        % faultTypes determine faulting ftype from the s1_east,s2_west,s3_north plunge values
        % returns a string description
        ftype = 'Unknow';

        if s1_east>52                         && s3_north<35
            ftype = 'Normal';
        end

        if 40<s1_east  && s1_east<52               && s3_north<20
            ftype = 'Normal to Strike Slip';
        end

        if s1_east<40            && s2_west>45     && s3_north<20
            ftype = 'Strike Slip';
        end

        if s1_east<20            && s2_west>45     && s3_north<40
            ftype = 'Strike Slip';
        end

        if s1_east<20                         && 40<s3_north && s3_north<52
            ftype = 'Thrust to Strike Slip';
        end

        if s1_east<35                         && s3_north>52
            ftype = 'Thrust';
        end
    end

end

author: CelsoReyes


text:
function M = ft_omri_volume_to_mosaic(V)

% function M = ft_omri_volume_to_mosaic(V)
%
% Reshuffle [MxNxS] volume to [XxY] mosaic

% Copyright (C) 2010, Stefan Klanke
%
% 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$

[h,w,ns] = size(V);
nn = ceil(sqrt(ns));
mm = ceil(ns/nn);

M = zeros(w*mm, h*nn, class(V));
row = 1;
col = 1;
for n = 1:ns
	is = 1 + (row-1)*h;
	ie = row*h;
	js = 1 + (col-1)*w;
	je = col*w;

	M(is:ie, js:je) = rot90(V(:,:,n));
	col = col + 1;
	if col>nn
	   row = row + 1;
	   col = 1;
	end
end


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

%%Transformation of the end effector frame following the RPY angles in
%%closed form
%%input: rotation angles along y,z,x axis expressed in degrees

function [RPYmat]=RPY(fia,fio,fin)

RPYmat=[cosd(fia)*cosd(fio),cosd(fia)*sind(fin)*sind(fio)-sind(fia)*cosd(fin),cosd(fia)*sind(fio)*cosd(fin)+sind(fia)*sind(fin),0;
        sind(fia)*cosd(fio),sind(fia)*sind(fio)*sind(fin)+cosd(fia)*cosd(fin),sind(fia)*sind(fio)*cosd(fin)-cosd(fia)*sind(fin),0;
        -sind(fio),cosd(fio)*sind(fin),cosd(fio)*cosd(fin),0;
        0,0,0,1];

end


text:
classdef TextDetectorCNN < handle
    %TEXTDETECTORCNN  Class providing functionality of text detection
    %
    % This class provides the functionality of text bounding box detection
    % (finds bounding boxes of text words given an input image).
    %
    % This class uses OpenCV dnn module to load pre-trained model described in
    % [LiaoSBWL17].
    % The original repository with the modified SSD Caffe version:
    % [GitHub](https://github.com/MhLiao/TextBoxes).
    %
    % Model (90.68 MB) can be downloaded from
    % [DropBox](https://www.dropbox.com/s/g8pjzv2de9gty8g/TextBoxes_icdar13.caffemodel?dl=0).
    %
    % Modified `.prototxt` file with the model description can be obtained from
    % [GitHub](https://github.com/opencv/opencv_contrib/raw/3.4.1/modules/text/samples/textbox.prototxt).
    %
    % ## References
    % [LiaoSBWL17]:
    % > Minghui Liao, Baoguang Shi, Xiang Bai, Xinggang Wang, and Wenyu Liu.
    % > "Textboxes: A fast text detector with a single deep neural network".
    % > In AAAI, 2017.
    %
    % See also: cv.TextDetectorCNN.TextDetectorCNN, cv.Net
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    %% TextDetectorCNN
    methods
        function this = TextDetectorCNN(varargin)
            %TEXTDETECTORCNN  Creates an instance using the provided parameters
            %
            %     obj = cv.TextDetectorCNN(modelArchFilename, modelWeightsFilename)
            %     obj = cv.TextDetectorCNN(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __modelArchFilename__ path to the prototxt file describing the
            %   classifiers architecture (the `.prototxt` file).
            % * __modelWeightsFilename__ path to the file containing the
            %   pretrained weights of the model in caffe-binary form
            %   (the `.caffemodel` file).
            %
            % ## Options
            % * __DetectionSizes__ a list of sizes for multiscale detection.
            %   The values
            %   `{[300,300], [700,500], [700,300], [700,700], [1600,1600]}`
            %   are recommended in [LiaoSBWL17] to achieve the best quality.
            %   default `{[300,300]}`.
            %
            % See also: cv.TextDetectorCNN
            %
            this.id = TextDetectorCNN_(0, 'new', varargin{:});
        end

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

    %% TextDetector
    methods
        function [bbox, conf] = detect(this, img)
            %DETECT  Detect text inside an image
            %
            %     [bbox, conf] = obj.detect(img)
            %
            % ## Input
            % * __img__ an image to process, expected to be a 8-bit color
            %   image of any size.
            %
            % ## Output
            % * __bbox__ a vector of rectangles of the detected word bounding
            %   boxes.
            % * __conf__ a vector of float of the confidences the classifier
            %   has for the corresponding bounding boxes.
            %
            % See also: cv.TextDetectorCNN.TextDetectorCNN
            %
            [bbox, conf] = TextDetectorCNN_(this.id, 'detect', img);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [c,ceq] = mycon(x,a2)
c = a2/x(1) - x(2);
ceq = x(1)-x(2);


text:
clc
clear all
names = {'DenseFuse', 'RFN-Nest', 'FusionGAN', 'SeAFusion', 'PIAFusion', 'IFCNN', 'PMGI', 'SDNet', 'U2Fusion'};
rows = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I'] ;
easy = 1; %% easy=1 用于测试：EN, SF,SD,PSNR,MSE, MI, VIF, AG, CC, SCD, Qabf等指标； easy=0 用于测试：Nabf, SSIM, MS_SSIM, FMI_pixel, FMI_dct, FMI_w等指标
dataset = 'TNO'
row_name1 = 'row1';
row_data1 = 'row2';
for i = 1 : length(names)
    method_name = cellstr(names(i));
    row = rows(i);
    row_name = strrep(row_name1, 'row', row);
    row_data = strrep(row_data1, 'row', row);
    fileFolder=fullfile('../Image/Source-Image', dataset, 'ir'); % 源图像A所在文件夹 此处是'Evaluation\Image\Source-Image\TNO\ir'
    dirOutput=dir(fullfile(fileFolder,'*.*'));
    fileNames = {dirOutput.name};
    [m, num] = size(fileNames);
    ir_dir = fullfile('../Image/Source-Image', dataset, 'ir'); % 源图像A所在文件夹 此处是'Evaluation\Image\Source-Image\TNO\ir'
    vi_dir = fullfile('../Image/Source-Image', dataset, 'vi'); % 源图像B所在文件夹 此处是'Evaluation\Image\Source-Image\TNO\vi'
    Fused_dir = '../';
    Fused_dir = fullfile(Fused_dir, 'Image', 'Algorithm', strcat(cell2mat(names(i)), '_', dataset)); % 融合结果所在文件夹 此处是 'Evaluation\Image\Algorithm\DenseFuse_TNO'
    EN_set = [];    SF_set = [];SD_set = [];PSNR_set = [];
    MSE_set = [];MI_set = [];VIF_set = []; AG_set = [];
    CC_set = [];SCD_set = []; Qabf_set = [];
    SSIM_set = []; MS_SSIM_set = [];
    Nabf_set = [];FMI_pixel_set = [];
    FMI_dct_set = []; FMI_w_set = [];
     for j = 1:num
        if (isequal(fileNames{j}, '.') || isequal(fileNames{j}, '..'))
            continue;
        else
            fileName_source_ir = fullfile(ir_dir, fileNames{j});
            fileName_source_vi = fullfile(vi_dir, fileNames{j});
            fileName_Fusion = fullfile(Fused_dir, fileNames{j});
            ir_image = imread(fileName_source_ir);
            vi_image = imread(fileName_source_vi);
            fused_image   = imread(fileName_Fusion);
            if size(ir_image, 3)>2
                ir_image = rgb2gray(ir_image);
            end

            if size(vi_image, 3)>2
                vi_image = rgb2gray(vi_image);
            end

            if size(fused_image, 3)>2
                fused_image = rgb2gray(fused_image);
            end

            [m, n] = size(fused_image);
        %     fused_image = fused_image(7:m-6, 7:n-6);
            ir_size = size(ir_image);
            vi_size = size(vi_image);
            fusion_size = size(fused_image);
            if length(ir_size) < 3 && length(vi_size) < 3
                [EN, SF,SD,PSNR,MSE, MI, VIF, AG, CC, SCD, Qabf, Nabf, SSIM, MS_SSIM, FMI_pixel, FMI_dct, FMI_w] = analysis_Reference(fused_image,ir_image,vi_image, easy);
                EN_set = [EN_set, EN];SF_set = [SF_set,SF];SD_set = [SD_set, SD];PSNR_set = [PSNR_set, PSNR];
                MSE_set = [MSE_set, MSE];MI_set = [MI_set, MI]; VIF_set = [VIF_set, VIF];
                AG_set = [AG_set, AG]; CC_set = [CC_set, CC];SCD_set = [SCD_set, SCD];
                Qabf_set = [Qabf_set, Qabf]; Nabf_set = [Nabf_set, Nabf];
                SSIM_set = [SSIM_set, SSIM]; MS_SSIM_set = [MS_SSIM_set, MS_SSIM];
                FMI_pixel_set = [FMI_pixel_set, FMI_pixel]; FMI_dct_set = [FMI_dct_set,FMI_dct];
                FMI_w_set = [FMI_w_set, FMI_w];
            else
                disp('unsucessful!')
                disp( fileName_Fusion)
            end

            fprintf('Fusion Method:%s, Image Name: %s\n', cell2mat(names(i)), fileNames{j})
        end
    end
    save_dir = '../Metric'; %存放Excel结果的文件夹
    if exist(save_dir,'dir')==0
        mkdir(save_dir);
    end
    %% 将测试结果写入 Excel， 此处采用writetable， 第一行可能会有问题，算法名在第二行，评估结果从第三行开始
    file_name = fullfile(save_dir, strcat('Metric_', dataset, '.xlsx')); %存放Excel文件的文件名
    if easy ==1
        SD_table = table(SD_set');
        PSNR_table = table(PSNR_set');
        MSE_table = table(MSE_set');
        MI_table = table(MI_set');
        VIF_table = table(VIF_set');
        AG_table = table(AG_set');
        CC_table = table(CC_set');
        SCD_table = table(SCD_set');
        EN_table = table(EN_set');
        Qabf_table = table(Qabf_set');
        SF_table = table(SF_set');
        method_table = table(method_name);

        writetable(SD_table,file_name,'Sheet','SD','Range',row_data);
        writetable(PSNR_table,file_name,'Sheet','PSNR','Range',row_data);
        writetable(MSE_table,file_name,'Sheet','MSE','Range',row_data);
        writetable(MI_table,file_name,'Sheet','MI','Range',row_data);
        writetable(VIF_table,file_name,'Sheet','VIF','Range',row_data);
        writetable(AG_table,file_name,'Sheet','AG','Range',row_data);
        writetable(CC_table,file_name,'Sheet','CC','Range',row_data);
        writetable(SCD_table,file_name,'Sheet','SCD','Range',row_data);
        writetable(EN_table,file_name,'Sheet','EN','Range',row_data);
        writetable(Qabf_table,file_name,'Sheet','Qabf','Range',row_data);
        writetable(SF_table,file_name,'Sheet','SF','Range',row_data);

        writetable(method_table,file_name,'Sheet','SD','Range',row_name);
        writetable(method_table,file_name,'Sheet','PSNR','Range',row_name);
        writetable(method_table,file_name,'Sheet','MSE','Range',row_name);
        writetable(method_table,file_name,'Sheet','MI','Range',row_name);
        writetable(method_table,file_name,'Sheet','VIF','Range',row_name);
        writetable(method_table,file_name,'Sheet','AG','Range',row_name);
        writetable(method_table,file_name,'Sheet','CC','Range',row_name);
        writetable(method_table,file_name,'Sheet','SCD','Range',row_name);
        writetable(method_table,file_name,'Sheet','EN','Range',row_name);
        writetable(method_table,file_name,'Sheet','Qabf','Range',row_name);
        writetable(method_table,file_name,'Sheet','SF','Range',row_name);
    else
        Nabf_table = table(Nabf_set');
        SSIM_table = table(SSIM_set');
        MS_SSIM_table = table(MS_SSIM_set');
        FMI_pixel_table = table(FMI_pixel_set');
        FMI_dct_table = table(FMI_dct_set');
        FMI_w_table = table(FMI_w_set');
        method_table = table(method_name);

        writetable(Nabf_table,file_name,'Sheet','Nabf','Range',row_data);
        writetable(SSIM_table,file_name,'Sheet','SSIM','Range',row_data);
        writetable(MS_SSIM_table,file_name,'Sheet','MS_SSIM','Range',row_data);
        writetable(FMI_pixel_table,file_name,'Sheet','FMI_pixel','Range',row_data);
        writetable(FMI_dct_table,file_name,'Sheet','FMI_dct','Range',row_data);
        writetable(FMI_w_table,file_name,'Sheet','FMI_w','Range',row_data);

        writetable(method_table,file_name,'Sheet','Nabf','Range',row_name);
        writetable(method_table,file_name,'Sheet','SSIM','Range',row_name);
        writetable(method_table,file_name,'Sheet','MS_SSIM','Range',row_name);
        writetable(method_table,file_name,'Sheet','FMI_pixel','Range',row_name);
        writetable(method_table,file_name,'Sheet','FMI_dct','Range',row_name);
        writetable(method_table,file_name,'Sheet','FMI_w','Range',row_name);

    end
end

author: Linfeng-Tang

repo: Image-Fusion


text:
function rect = genbinarymap(imglabel, salmap)
imx = processing(salmap);
imx = imresize(imx,[size(imglabel,1) size(imglabel,2)]);
stats = regionprops(imx, 'BoundingBox');
rect = stats.BoundingBox;
figure;
imshow(imglabel);
hold on;
rectangle('Position',rect,'EdgeColor','r', 'LineWidth',4);
F = getframe;
imwrite(F.cdata, 'Object.jpg', 'jpg');
hold off;

% binarization
function X = processing(X, minarea)
level = graythresh(X);
X1 = im2bw(X,level);
X1 = imfill(X1,'holes');
if nargin<2,
    minarea = length(find(X1))*0.01;   %
end
BW1 = bwlabel(X1,4);
nregions = max(BW1(:));
for i = 1:nregions,
    idxi = logical(BW1==i);
    idxi_n = find(BW1==i);
    %X_gain=sum(X(idxi))/length(idxi_n);
    X_gain=max(X(idxi));
    if i ==1
        X_gain_old=X_gain;
        idxi_n_old=idxi_n;
    else
        if X_gain > X_gain_old
            X1(idxi_n_old) = 0;
            X_gain_old=X_gain;
            idxi_n_old=idxi_n;
        else
            X1(idxi_n) = 0;
        end
    end
end
X=X1;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ r, seed ] = r8vec_uniform_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real R(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  r = zeros ( n, 1 );

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM_01 - Fatal error!' );
  end

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + 2147483647;
    end

    r(i) = seed * 4.656612875E-10;

  end

  return
end


text:
function tno_order_test ( )

%*****************************************************************************80
%
%% TNO_ORDER_TEST tests TNO_ORDER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TNO_ORDER_TEST:\n' );
  fprintf ( 1, '  TNO_ORDER is given a level L, and returns the order N of\n' );
  fprintf ( 1, '  a Truncated Normal Odd (TNO) quadrature rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   Level   Order\n' );
  fprintf ( 1, '\n' );

  for l = 1 : 20

    n = tno_order ( l );

    fprintf ( 1, '  %2d     %6d\n', l, n )

  end

  return
end


text:
%% Recover function p from P0 to P1
% [rp,nodenew,elemnew] = recoverP02P1(node,elem,p,recoverMethod)
% The recover methods
% 1) 'LS': least square
% 2) 'LA': solve the Laplacian problem, only this method
%           the new node and new elem will be constructed.
%
% Created by Lin Zhong April, 2013.

%% Node and Elem
clear all; close all;
[node,elem] = squaremesh([0,1,0,1],0.5);
bdFlag = setboundary(node,elem,'Dirichlet');
figure(1)
showmesh(node,elem);
findnode(node);
findelem(node,elem);
display(elem);

NT = size(elem,1);
p = ones(NT,1);
[rp,nodenew,elemnew] = recoverP02P1(node,elem,p,'LA');
figure(2)
showmesh(nodenew,elemnew);
findnode(nodenew);
findelem(nodenew,elemnew);
display(elemnew);

%% LS and LA
% For LS, the value at the nodes are constructed by least square linear
% fitting of the barycenters of the adjacent elemts.
%
% For LA, we construct new nodes and new elems. The new nodes are the
% vertices of the triangles and the barycenters of the elemenst. The new
% elems are obtained by connecting two barycenters of the interior edges,
% and each end point of that edge. Please check the graph.
% We get the values of the interior vertices by solving the Laplacian
% problem. To improve the efficiency, we only construct part of the
% stiffness matrix.
% $$ A_{N,N+NT}(rp;p)^t = 0$$
% we solve $$rp = -A_{N,N}\A(N,N+NT)*p$$
% We notice that $A_{N,N}$ is diagonal, so this process should be fast.
% For the boundary nodes, we use least square.
%
% Both methods achieve the same accuracy.
% Because of the large size for loop, the LS method is slower than LA.
[node,elem] = squaremesh([0,1,0,1],0.05);
NT = size(elem,1);
p = ones(NT,1);
profile on
recoverP02P1(node,elem,p,'LA');
profile viewer

repo: Aero-matlab


text:
%  Utility function that transforms a string whose evaluation would give a
%  matrix, into a cell array.
%
%  designed for use with Ben Petschel stuff on Grobner basis
%

author: jmaih


text:
function determ = kms_determinant ( alpha, n )

%*****************************************************************************80
%
%% KMS_DETERMINANT returns the determaint of the KMS matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar that defines A.
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  if ( n == 1 )
    determ = 1.0
  else
    determ = ( 1.0 - alpha * alpha )^( n - 1 );
  end

  return
end


text:
function days = month_length_eg_civil ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_EG_CIVIL returns the number of days in an Egyptian Civil month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer DAYS, the number of days
%    in the month.
%
  if ( m < 1 )
    days = 0;
  elseif ( m <= 12 )
    days = 30;
  elseif ( m == 13 )
    days = 5;
  else
    days = 0;
  end

  return
end


text:
%%% ERROR_ANALYSIS
%%% This simulation helps coputing the acturacies of calibration
%%% Run it after the main calibration

N_runs = 200;

%N_ima_active = 4;

saving = 1;

if 1, %~exist('fc_list'), % initialization

   % Initialization:

   load Calib_Results;
   check_active_images;

	fc_list = [];
	cc_list = [];
	kc_list = [];
   active_images_list = [];

	for kk=1:n_ima,

   	eval(['omc_list_' num2str(kk) ' = [];']);
   	eval(['Tc_list_' num2str(kk) ' = [];']);

   end;

	%sx = median(abs(ex(1,:)))*1.4836;
	%sy = median(abs(ex(2,:)))*1.4836;

	sx = std(ex(1,:));
	sy = std(ex(2,:));

	% Saving the feature locations:

	for kk = 1:n_ima,

   	eval(['x_save_' num2str(kk) ' = x_' num2str(kk) ';']);
   	eval(['y_save_' num2str(kk) ' = y_' num2str(kk) ';']);

	end;

   active_images_save = active_images;
   ind_active_save = ind_active;

   fc_save = fc;
   cc_save = cc;
   kc_save = kc;
   KK_save = KK;

end;

%%% The main loop:

for ntrial = 1:N_runs,

   fprintf(1,'\nRun number: %d\n',ntrial);
   fprintf(1,  '----------\n');

   for kk = 1:n_ima,

      eval(['y_kk = y_save_' num2str(kk) ';'])

      if active_images(kk) & ~isnan(y_kk(1,1)),

         Nkk = size(y_kk,2);

         x_kk_new = y_kk + [sx * randn(1,Nkk);sy*randn(1,Nkk)];

         eval(['x_' num2str(kk) ' = x_kk_new;']);

      end;

   end;

   N_active = length(ind_active_save);
   junk = randn(1,N_active);
   [junk,junk2] = sort(junk);

   active_images = zeros(1,n_ima);
   active_images(ind_active_save(junk2(1:N_ima_active))) = ones(1,N_ima_active);

   fc = fc_save;
   cc = cc_save;
   kc = kc_save;
   KK = KK_save;

   go_calib_optim;

   fc_list = [fc_list fc];
   cc_list = [cc_list cc];
   kc_list = [kc_list kc];
   active_images_list = [active_images_list active_images'];

   for kk=1:n_ima,

   	eval(['omc_list_' num2str(kk) ' = [ omc_list_' num2str(kk) ' omc_' num2str(kk) ' ];']);
   	eval(['Tc_list_' num2str(kk) ' = [ Tc_list_' num2str(kk) ' Tc_' num2str(kk) ' ];']);

	end;

end;

if 0,

% Restoring the feature locations:

for kk = 1:n_ima,

   eval(['x_' num2str(kk) ' = x_save_' num2str(kk) ';']);

end;

fprintf(1,'\nFinal run (with the real data)\n');
fprintf(1,  '------------------------------\n');

active_images = active_images_save;
ind_active = ind_active_save;

go_calib_optim;

fc_list = [fc_list fc];
cc_list = [cc_list cc];
kc_list = [kc_list kc];
active_images_list = [active_images_list active_images'];

for kk=1:n_ima,

   eval(['omc_list_' num2str(kk) ' = [ omc_list_' num2str(kk) ' omc_' num2str(kk) ' ];']);
   eval(['Tc_list_' num2str(kk) ' = [ Tc_list_' num2str(kk) ' Tc_' num2str(kk) ' ];']);

end;

end;

if saving,

disp(['Save Calibration accuracy results under Calib_Accuracies_' num2str(N_ima_active) '.mat']);

string_save = ['save Calib_Accuracies_' num2str(N_ima_active) ' active_images n_ima N_ima_active N_runs active_images_list fc cc kc fc_list cc_list kc_list'];

for kk = 1:n_ima,
   string_save = [string_save ' Tc_list_' num2str(kk) ' omc_list_' num2str(kk)  ' Tc_' num2str(kk) ' omc_' num2str(kk) ];
end;

eval(string_save);

end;

return;

std(fc_list')

std(cc_list')

std(kc_list')

for kk = 1:n_ima,

   eval(['std(Tc_list_'  num2str(kk) ''')'])

end;

repo: ekfmonoslam


text:
function [cmap] = makeColorGradient(col1, col2, ncol)
% Generates a color gradient in hex format, based on color 1 and 2
%
% USAGE:
%
%    [cmap] = makeColorGradient(col1, col2, ncol)
%
% INPUTS:
%    col1:                      First color
%    col2:                      Last color
%    ncol:                      Number of colors in between
%
% OUTPUT:
%    cmap:                      Color gradient map with the corresponding
%                               colors in hex format
%
% .. Author: - Nicolas Mendoza-Mejia May/2021

hexformat = "#%02x%02x%02x";

rgb1 = sscanf(col1, hexformat);
rgb2 = sscanf(col2, hexformat);

T = (rgb2 -rgb1)/(ncol - 1);
rgb = zeros(3, ncol);

for i=1:3
   if T(i) ~= 0
    rgb(i, :) =  rgb1(i):T(i):rgb2(i);
   end
end

rgb = round(rgb);
cmap = cell(1, ncol);
for i=1:ncol
   cmap{i} = sprintf(hexformat, rgb(1, i), rgb(2, i), rgb(3, i));
end

end

author: opencobra

repo: cobratoolbox