text:
function hg=gauss_kernel(k,nd,f,m)

% function hg=gauss_kernel(k,nd,f,m)
% ------------------------------------------------------------------------
%
%
% ------------------------------------------------------------------------

%%
switch m
    case {'sigma',1}
        %NOTE: This method is equivalent to using hg = fspecial('gaussian',[k k], S)
        S=f3;
        switch nd
            case 1
                x = linspace(-((k-1)/2),((k-1)/2),k);
                hg=exp(-(x.^2)./(2*S^2));
            case 2
                [x,y] = meshgrid(linspace(-((k-1)/2),((k-1)/2),k));
                hg=exp(-(x.^2 + y.^2)./(2*S^2));
            case 3
                [x,y,z] = meshgrid(linspace(-((k-1)/2),((k-1)/2),k));
                hg=exp(-(x.^2 + y.^2 + z.^2)./(2*S^2));
        end
    case {'width',2}
        % Here "f" defines "where the bell curve is" at the edges of the
        % kernel e.g. f=2 results in twice the standard deviation.
        switch nd
            case 1
                x = linspace(-f,f,k);
                hg=exp(-(x.^2)./2);
            case 2
                [x,y] = meshgrid(linspace(-f,f,k));
                hg=exp(-(x.^2 + y.^2)./2);
            case 3
                [x,y,z] = meshgrid(linspace(-f,f,k));
                hg=exp(-(x.^2 + y.^2 + z.^2)./2);
        end
    otherwise
        warning('False input for method! Valid options are "sigma" and "width"');
end

hg(hg<eps*max(hg(:)))=0;

%Making sure mask sums to 1
hg=hg./sum(hg(:));

end

%%
% _*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 GraphCost
    properties
        nCost = 1;

        totalCost = [];     % total cost of all cost nodes and also regularization cost
        taskCost = [];      % total cost of all cost nodes without regularization cost
        reguCost = [];      % total regularization cost
        subCost = [];       % cost of each cost nodes
        subAcc = [];        % accuracy of cost classification nodes
    end

    methods
        function obj = GraphCost(nCost, nEntry)
            if nargin<2; nEntry = 1; end
            obj.totalCost = zeros(1,nEntry);
            obj.taskCost = zeros(nCost,nEntry);
            obj.reguCost = zeros(1,nEntry);
            obj.subCost = zeros(nCost,nEntry);
            obj.subAcc = zeros(nCost,nEntry);
        end

        % copy the cost to an outside GraphCost object
        % index specifies the location of the destination.
        function costDestination = copyCost(obj, costDestination, index)
            if length(obj.totalCost)~=length(index)
                fprintf('GraphCost: Error: the number of destination index is not equal to the number of cost in this object\n');
            end
            costDestination.totalCost(index) = obj.totalCost;
            costDestination.taskCost(:, index) = obj.taskCost;
            costDestination.reguCost(index) = obj.reguCost;
            costDestination.subCost(:, index) = obj.subCost;
            costDestination.subAcc(:, index) = obj.subAcc;
        end

        % append the cost of this object to an outside GraphCost object.
        function costDestination = appendCost(obj, costDestination)
            nEntry = length(obj.totalCost);
            nEntryInDestination = length(costDestination.totalCost);
            index = nEntryInDestination+1 : nEntryInDestination+nEntry;
            costDestination = obj.copyCost(costDestination, index);
        end

        function costDestination = appendMeanCost(obj, costDestination)
            costDestination.totalCost(end+1) = mean(obj.totalCost);
            costDestination.taskCost(:, end+1) = mean(obj.taskCost,2);
            costDestination.reguCost(end+1) = mean(obj.reguCost);
            costDestination.subCost(:, end+1) = mean(obj.subCost,2);
            costDestination.subAcc(:, end+1) = mean(obj.subAcc,2);
        end
    end

end

author: singaxiong

repo: SignalGraph


text:
% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file
% Coupled nonlinear PDE's
% Buckley Leverett equation
% dependent variables: pressure and water saturation
% Prepared for educational purposes by ** AAE **
% Written by Ali A. Eftekhari
% Last checked: June 2021
% This code works fine (two phase flow + Pc) and shows the effect of Pc
% in the simulation of a heterogeneous 2D field.
clc
%% define the geometry
Nx = 100; % number of cells in x direction
Ny = 30; % number of cells in y direction
W = 300; % [m] length of the domain in x direction
H = 30; % [m] length of the domain in y direction
% m = createMesh1D(Nx, W);
m = createMesh2D(Nx, Ny, W, H); % creates a 2D mesh
%% define the physical parametrs
krw0 = 1.0;
kro0 = 0.76;
nw = 2.4;
no = 2.0;
sor=0.12;
swc=0.09;
sws=@(sw)((sw>swc).*(sw<1-sor).*(sw-swc)/(1-sor-swc)+(sw>=1-sor).*ones(size(sw)));
kro=@(sw)((sw>=swc).*kro0.*(1-sws(sw)).^no+(sw<swc).*(1+(kro0-1)/swc*sw));
krw=@(sw)((sw<=1-sor).*krw0.*sws(sw).^nw+(sw>1-sor).*(-(1-krw0)/sor.*(1.0-sw)+1.0));
dkrwdsw=@(sw)((sw<=1-sor).*nw.*krw0.*(1/(1-sor-swc)).*sws(sw).^(nw-1)+(sw>1-sor)*((1-krw0)/sor));
dkrodsw=@(sw)((sw>=swc).*(-kro0*no*(1-sws(sw)).^(no-1))/(-swc-sor+1)+(sw<swc).*((kro0-1)/swc));
p0 = 100e5; % [bar] pressure
pin = 150e5; % [bar] injection pressure at the left boundary
u_in= 1.0/(24*3600); % [m/s] equal to 1 m/day
sw0 = swc+0.1; % initial water saturation
sw_in = 1;
mu_oil = 2e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 2e-12; % [m^2] average reservoir permeability
phi0 = 0.2; % average porosity
teta_ow=deg2rad(30);
gama_ow=0.03; % N/m
labda=10.0;
eps1=1e-6;
clx=0.2;
cly=0.2;
V_dp=0.7; % Dykstra-Parsons coef.
perm_val= field2d(Nx,Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, perm_val);
phi=createCellVariable(m, phi0);
pce=gama_ow*cos(teta_ow)*(phi./k).^0.5;
pce.value(isnan(pce.value)) = 1.0;
pce_face=gama_ow*cos(teta_ow)*(arithmeticMean(phi)./geometricMean(k)).^0.5;
pc=@(sw)(pce.*(sws(sw)+eps).^(-1/labda)); % it can also be defined for each block
dpc=@(sw)((-1/labda)*(1/(1-sor-swc)).*pce_face.*(sws(sw)+eps).^(-1/labda-1));
dpcdk=@(sw)(0.5*gama_ow*cos(teta_ow)*(geometricMean(k)./arithmeticMean(phi)).^0.5.*(sws(sw)+eps).^(-1/labda));
grad_phik=gradientTerm(phi./k);
% sw_plot=linspace(0,1, 10000);
% plot(sw_plot, pc(sw_plot))
lw = geometricMean(k)/mu_water;
lo = geometricMean(k)/mu_oil;
Lw = @(sw)(krw(sw));
Lo = @(sw)(k/mu_oil*kro(sw));
dLwdsw = @(sw)(k/mu_water*dkrwdsw(sw));
dLodsw = @(sw)(k/mu_oil*dkrodsw(sw));
%% Define the boundaries
BCp = createBC(m); % Neumann BC for pressure
BCs = createBC(m); % Neumann BC for saturation
% left boundary pressure gradient
BCp.left.a(:)=(krw(sw_in)*lw.xvalue(1,:)+kro(sw_in)*lo.xvalue(1,:)); BCp.left.b(:)=0; BCp.left.c(:)=-u_in;
% change the right boandary to constant pressure (Dirichlet)
% BCp.left.a(:)=0; BCp.left.b(:)=1; BCp.left.c(:)=pin;
BCp.right.a(:)=0; BCp.right.b(:)=1; BCp.right.c(:)=p0;
% BCs.right.a(:)=0; BCs.right.b(:)=1; BCs.right.c(:)=0.8; % end effect
% change the left boundary to constant saturation (Dirichlet)
BCs.left.a(:)=0; BCs.left.b(:)=1; BCs.left.c(:)=1;
%% define the time step and solver properties
% dt = 1000; % [s] time step
dt=(W/Nx)/u_in/20; % [s]
t_end = 1000*dt; % [s] final time
eps_p = 1e-5; % pressure accuracy
eps_sw = 1e-5; % saturation accuracy
%% define the variables
sw_old = createCellVariable(m, sw0, BCs);
p_old = createCellVariable(m, p0, BCp);
sw = sw_old;
p = p_old;
uw = -gradientTerm(p_old); % an estimation of the water velocity
%% start the main loop
% generate intial pressure profile (necessary to initialize the fully
% implicit solver)

t = 0;
while (t<t_end)
    error_p = 1e5;
    error_sw = 1e5;
    % Implicit loop
    while ((error_p>eps_p) || (error_sw>eps_sw))
        % calculate parameters
        pgrad = gradientTerm(p);
%         pcgrad=gradientTerm(pc(sw));
        sw_face = upwindMean(sw, -pgrad); % average value of water saturation
        sw_grad=gradientTerm(sw);
        sw_ave=arithmeticMean(sw);
        pcgrad=dpc(sw_ave).*sw_grad+dpcdk(sw_ave).*grad_phik;
        labdao = lo.*funceval(kro, sw_face);
        labdaw = lw.*funceval(krw, sw_face);
        dlabdaodsw = lo.*funceval(dkrodsw, sw_face);
        dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
        labda = labdao+labdaw;
        dlabdadsw = dlabdaodsw+dlabdawdsw;
        % compute [Jacobian] matrices
        Mdiffp1 = diffusionTerm(-labda);
        Mdiffp2 = diffusionTerm(-labdaw);
%         Mconvsw1 = convectionUpwindTerm(-dlabdadsw.*pgrad); % without capillary
        Mconvsw1 = convectionUpwindTerm(-dlabdawdsw.*pgrad-dlabdaodsw.*(pgrad+pcgrad)); % with capillary
        Mconvsw2 = convectionUpwindTerm(-dlabdawdsw.*pgrad);
        [Mtranssw2, RHStrans2] = transientTerm(sw_old, dt, phi);
        % Compute RHS values
        RHSpc1=divergenceTerm(labdao.*pcgrad);
        RHS1 = divergenceTerm((-dlabdawdsw.*pgrad-dlabdaodsw.*(pgrad+pcgrad)).*sw_face); % with capillary
        RHS2 = divergenceTerm(-dlabdawdsw.*sw_face.*pgrad);
        % include boundary conditions
        [Mbcp, RHSbcp] = boundaryCondition(BCp);
        [Mbcsw, RHSbcsw] = boundaryCondition(BCs);
        % Couple the equations; BC goes into the block on the main diagonal
        M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw];
        RHS = [RHS1+RHSpc1+RHSbcp; RHS2+RHStrans2+RHSbcsw];
        % solve the linear system of equations
        x = M\RHS;
        % x = agmg(M, RHS, [], 1e-10, 500, [], [p.value(:); sw.value(:)]);
        % separate the variables from the solution
        p_new = reshapeCell(m,full(x(1:(Nx+2)*(Ny+2))));
        sw_new = reshapeCell(m,full(x((Nx+2)*(Ny+2)+1:end)));
        % calculate error values
        error_p = max(abs((p_new(:)-p.value(:))./p_new(:)));
        error_sw = max(abs(sw_new(:)-sw.value(:)));
        % assign new values of p and sw
        p.value = p_new;
        sw.value = sw_new;
    end
    t=t+dt;
    p_old = p;
    sw_old = sw;
    figure(1);visualizeCells(sw); drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function Hv = hessvec_fd(v,FUN,x,gx,s)
%HESSVEC_FD   Hessian vector product finite difference approximation.
%
%   HV = HESSVEC_FD(V,FUN,X) computes a forward finite difference
%   approximation of the Hessian vector product, H(X)*V, where V is given
%   and H(X) is the Hessian of the function FUN at the point X. The
%   approximation is given by
%
%              G(X+S*V) - G(X)
%     H(X)*V = ---------------
%                     S
%
%   where G(X) is the gradient of the function FUN at the point X and
%   h is the difference step (default = 1e-8*(1+norm(X)).
%
%   HV = HESSVEC(V,FUN,X,GX) uses GX for the value of G(X), in the
%   case that is has already been computed outside of this method.
%
%   HV = HESSVEC(V,FUN,X,GX,S) uses S as the difference step.
%
%   This method should not be called directly, but only by Poblano Toolbox
%   optimization algorithms.
%
%   The number of calls to FUN is kept track of in the global variable
%   NFEV_HESSVEC_FD to facilitate using this method as a callback function
%   (e.g., from an iterative linear solver such as Matlab's SYMMLQ).
%
%   See also TN.
%
%MATLAB Poblano Toolbox.
%Copyright 2009-2012, Sandia Corporation.

%% Number of calls to FUN
global nfev_hessvec_fd;

%% Check input
if nargin < 3
    error('HESSVEC_FD => at least two input arguments are required (V,FUN,X)');
end

% Compute gradient if not given
if nargin < 4
    [f,gx] = feval(FUN,x);
    nfev_hessvec_fd = nfev_hessvec_fd + 1;
end

% Compute a difference step if not given
if nargin < 5
    s = 1e-8*(1+norm(x));
end

% Compute the gradient at the new point
[f,gxsv] = feval(FUN,x+s*v);
nfev_hessvec_fd = nfev_hessvec_fd + 1;

%% Hessian vector product approximation
Hv = (gxsv-gx)/s;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function sub_sst = extract_sst(sst,t1,t2,edge_mode)

%EXTRACT_SST: Extract event SST (Start/Stop Times) between t1 and t2
%
%USAGE: sub_sst = extract_sst(sst,t1,t2,edge_mode)
%
%INPUTS: sst - Event Start/Stop Times [Nx2 array of numeric time values]
%        t1  - time value, extract all events in sst between t1 (start)
%              and t2 (end)
%        t2  - see 't1'
%        edge_mode  - If t1 or t2 falls between an event start time and
%                stop time in sst edge_mode determines if the event is
%                deleted, retained, or partially deleted at time t1 or t2.
%           'delete' - whole event is removed from array sub_sst
%           'keep'   - whole event is retained in array sub_sst
%           'part'   - event in sub_sst is split - i.e. deleted
%                       before t1 or after t2
%
%OUTPUTS: sub_sst - Subset of origianl sst (Event Start/Stop Times)
%                   [Nx2 array of numeric time values]
%
% See also ADD_SST, CHK_T, COMPARE_SST, DELETE_SST, EXTRACT_SST, IS_SST,
%          ISTEQUAL, MERGE_SST, SEARCH_SST, SORT_SST, NAN2SST, SSD2SST,
%          SST2NAN, SST2SSD, SST2VAL, SST2WFA, WFA2SST
%
% Author: Dane Ketner, Alaska Volcano Observatory
% $Date$
% $Revision$

if nargin < 4
   error('EXTRACT_SST: Too few input arguments')
elseif nargin > 4
   error('EXTRACT_SST: Too many input arguments')
end

if is_sst(sst)
   if iscell(sst)
      for m = 1:numel(sst)
         sub_sst{m} = EXTRACT_SST(sst{m},t1,t2,edge_mode);
      end
   else
      sub_sst = EXTRACT_SST(sst,t1,t2,edge_mode);
   end
else
   error('EXTRACT_SST: Not a valid Start/Stop Time Argument')
end

%%
function sub_sst = EXTRACT_SST(sst,t1,t2,edge_mode)

[N1 P1] = search_sst(t1,sst);
[N2 P2] = search_sst(t2,sst);

% 'first' and 'last' refer to the first and last events within the time
%  span defined by t1 and t2. N1 is the event number within sst
%  corresponding to t1, and P1 is the corresponding event position with
%  relation to t1. P1 = 1 if t1 falls inside the event time. P1 = 0 if t1
%  is before the start of the event.

if P1 == 1
   if strcmpi(edge_mode,'part')
      first = [t1 sst(N1,2)];
   elseif strcmpi(edge_mode,'delete')
      first = [];
   elseif strcmpi(edge_mode,'keep')
      first = sst(N1,:);
   end
elseif P1 == 0
   first = sst(N1,:);
end

if P2 == 1
   if strcmpi(edge_mode,'part')
      last = [sst(N2,1) t2];
   elseif strcmpi(edge_mode,'delete')
      last = [];
   elseif strcmpi(edge_mode,'keep')
      last = sst(N2,:);
   end
elseif P2 == 0
   last = [];
end

sub_sst = [first; sst(N1+1:N2-1,:); last];

author: geoscience-community-codes

repo: GISMO


text:
function [H11, H12, H21, H22] = d2Abr_dV2(d2F_dV2, dF_dV1, dF_dV2, F, V, mu)
%D2ABR_DV2   Computes 2nd derivatives of |branch flow|^2 w.r.t. V.
%
%   The derivatives can be take with respect to polar or cartesian coordinates
%   of voltage, depending on the first 3 arguments. Flows could be complex
%   current or complex or real power. Notation below is based on complex power.
%
%   [H11, H12, H21, H22] = D2ABR_DV2(D2F_DV2, DF_DV1, DF_DV2, F, V, MU)
%
%   Returns 4 matrices containing the partial derivatives w.r.t. voltage
%   components (angle, magnitude or real, imaginary) of the product of a
%   vector MU with the 1st partial derivatives of the square of the magnitude
%   of branch flows.
%
%   Takes as inputs a handle to a function that evaluates the 2nd derivatives
%   of the flows (with args V and mu only), sparse first derivative matrices
%   of flow, flow vector, voltage vector V and nl x 1 vector of multipliers
%   MU. Output matrices are sparse.
%
%   Example:
%       f = branch(:, F_BUS);
%       Cf =  sparse(1:nl, f, ones(nl, 1), nl, nb);
%       [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
%       [dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = ...
%               dSbr_dV(branch, Yf, Yt, V);
%       dF_dV1 = dSf_dV1;
%       dF_dV2 = dSf_dV2;
%       F = Sf;
%       d2F_dV2 = @(V, mu)d2Sbr_dV2(Cf, Yf, V, mu, 0);
%       [H11, H12, H21, H22] = ...
%             d2Abr_dV2(d2F_dV2, dF_dV1, dF_dV2, F, V, mu);
%
%   Here the output matrices correspond to:
%     H11 = d/dV1 (dAF_dV1.' * mu)
%     H12 = d/dV2 (dAF_dV1.' * mu)
%     H21 = d/dV1 (dAF_dV2.' * mu)
%     H22 = d/dV2 (dAF_dV2.' * mu)
%
%   See also DABR_DV, DIBR_DV, DSBR_DV.
%
%   For more details on the derivations behind the derivative code used
%   in MATPOWER information, see:
%
%   [TN2]  R. D. Zimmerman, "AC Power Flows, Generalized OPF Costs and
%          their Derivatives using Complex Matrix Notation", MATPOWER
%          Technical Note 2, February 2010. [Online]. Available:
%          https://matpower.org/docs/TN2-OPF-Derivatives.pdf
%          doi: 10.5281/zenodo.3237866

%   MATPOWER
%   Copyright (c) 2008-2019, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define
nl = length(mu);

diagmu = sparse(1:nl, 1:nl, mu, nl, nl);

[F11, F12, F21, F22] = d2F_dV2(V, conj(F) .* mu);
H11 = 2 * real( F11 + dF_dV1.' * diagmu * conj(dF_dV1) );
H21 = 2 * real( F21 + dF_dV2.' * diagmu * conj(dF_dV1) );
H12 = 2 * real( F12 + dF_dV1.' * diagmu * conj(dF_dV2) );
H22 = 2 * real( F22 + dF_dV2.' * diagmu * conj(dF_dV2) );

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [major, minor, numCommits, sha] = parseVersion(ver)
    if ~isempty(strfind(lower(ver), 'fatal'))
        error('BNT:arg', 'Incorrect argument ''ver'' (type ''help helpers.parseVersion'' for details).');
    end
    hyphenInd = strfind(ver, '-'); % should be two
    if isempty(hyphenInd)
        nums = regexpi(ver, '\d+', 'match'); % check that we have some numbers (should be two of them)
        if isempty(nums)
            error('BNT:arg', 'Incorrect argument ''ver'' (type ''help helpers.parseVersion'' for details).');
        end
        major = str2double(nums{1});
        minor = str2double(nums{2});
        numCommits = '';
        sha = '';
    else
        if length(hyphenInd) ~= 2
            error('BNT:arg', 'Incorrect argument ''ver'' (type ''help helpers.parseVersion'' for details).');
        end
        verId = ver(1:hyphenInd(1)-1);
        numCommits = ver(hyphenInd(1)+1:hyphenInd(2)-1);
        sha = ver(hyphenInd(2)+2:end); % remove 'g' character
        nums = regexp(verId, '\d*', 'match');
        major = str2double(nums{1});
        minor = str2double(nums{2});
    end
end

author: kavli-ntnu


text:
function outputSummary(summaryTable, i, k, variation, shadRed)
% Internal function to print out a summaryTable
%
% USAGE:
%
%    outputSummary(summaryTable, i, k, variation, shadRed)
%
% INPUT:
%    summaryTable:   Cell containing the summary of shadow prices and reduced costs
%    i:              Row index in the summaryTable (solver)
%    k:              Column index in the summaryTable
%    variation:      String (can be either 'increase' or 'decrease')
%    shadRed:        String indicating whether the evaluation is a shadow price (SP) or reduced cost (RC)
%
% OUTPUT:
%    Prints an output to the console

    if summaryTable{i + 1, k} > 0
        fprintf([' + ', shadRed, ' is positive for metabolites that ', variation, ' OF flux\n']);
    elseif summaryTable{i + 1, k} < 0
        fprintf([' - ', shadRed, ' is negative for metabolites that ', variation, ' OF flux\n']);
    elseif summaryTable{i + 1, k} == 0
        fprintf([' * ', shadRed, ' is zero for metabolites that ', variation, ' OF flux\n']);
    end
end

author: opencobra

repo: cobratoolbox


text:
function [bankAng,angOfAttack,angOfSideslip,totalAoA] = computeAeroAnglesFromFrameBodyAxes(rVectFrame, vVectFrame, bodyXFrame, bodyYFrame, bodyZFrame)
    %Source: http://www.dept.aoe.vt.edu/~cdhall/courses/aoe5204/AircraftMotion.pdf

    [R_wind_2_frame, ~, ~, ~] = computeWindFrame(rVectFrame, vVectFrame);
    R_body_2_frame = horzcat(bodyXFrame, bodyYFrame, bodyZFrame);

    angles = real(rotm2eulARH(R_wind_2_frame' * R_body_2_frame, 'zyx'));

    bankAng = angles(3);
	angOfAttack = angles(2);
	angOfSideslip = angles(1);

    [x,y,z] = sph2cart(angleNegPiToPi(angOfSideslip),angleNegPiToPi(angOfAttack),1);
    v1 = [1;0;0];
    v2 = [x;y;z];
    totalAoA = dang(v1,v2);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = drawEllipse(varargin)
%DRAWELLIPSE Draw an ellipse on the current axis.
%
%   drawEllipse(ELLI);
%   Draws the ellipse ELLI in the form [XC YC RA RB THETA], with center
%   (XC, YC), with main axis of half-length RA and RB, and orientation
%   THETA in degrees counted counter-clockwise.
%
%   drawEllipse(XC, YC, RA, RB);
%   drawEllipse(XC, YC, RA, RB, THETA);
%   Specifies ellipse parameters as separate arguments (old syntax).
%
%   drawEllipse(..., NAME, VALUE);
%   Specifies drawing style of ellipse, see the help of plot function.
%
%   H = drawEllipse(...);
%   Also returns handles to the created line objects.
%
%   -> Parameters can also be arrays. In this case, all arrays are supposed
%   to have the same size.
%
%   Example:
%   % Draw an ellipse centered in [50 50], with semi major axis length of
%   % 40, semi minor axis length of 20, and rotated by 30 degrees.
%     figure(1); clf; hold on;
%     drawEllipse([50 50 40 20 30]);
%     axis equal; axis([0 100 10 90])
%
%   % add another ellipse with different orientation and style
%     drawEllipse([50 50 40 20 -10], 'LineWidth', 2, 'Color', 'g');
%
%
%   See also
%     ellipses2d, drawCircle, drawEllipseArc, drawEllipseAxes
%     fitEllipse, ellipseToPolygon, ellipsePoint, transformEllipse
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-12-11
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

%% Extract input arguments

% extract handle of axis to draw on
if isAxisHandle(varargin{1})
    ax = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% extract dawing style strings
styles = {};
for i = 1:length(varargin)
    if ischar(varargin{i})
        styles = varargin(i:end);
        varargin(i:end) = [];
        break;
    end
end

% extract ellipse parameters
if length(varargin) == 1
    % ellipse is given in a single array
    ellipse = varargin{1};
    x0 = ellipse(:, 1);
    y0 = ellipse(:, 2);
    a  = ellipse(:, 3);
    b  = ellipse(:, 4);
    if length(ellipse) > 4
        theta = ellipse(:, 5);
    else
        theta = zeros(size(x0));
    end

elseif length(varargin) >= 4
    % ellipse parameters given as separate arrays
    x0 = varargin{1};
    y0 = varargin{2};
    a  = varargin{3};
    b  = varargin{4};
    if length(varargin) > 4
        theta = varargin{5};
    else
        theta = zeros(size(x0));
    end

else
    error('drawEllipse: incorrect input arguments');
end

%% Process drawing of a set of ellipses

% angular positions of vertices
t = linspace(0, 2*pi, 145);

% compute position of points to draw each ellipse
h = zeros(length(x0), 1);
for i = 1:length(x0)
    % pre-compute rotation angles (given in degrees)
    cot = cosd(theta(i));
    sit = sind(theta(i));

    % compute position of points used to draw current ellipse
    xt = x0(i) + a(i) * cos(t) * cot - b(i) * sin(t) * sit;
    yt = y0(i) + a(i) * cos(t) * sit + b(i) * sin(t) * cot;

    % stores handle to graphic object
    h(i) = plot(ax, xt, yt, styles{:});
end

% return handles if required
if nargout > 0
    varargout = {h};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
clf; figure(1);

subplot('position',[0 0 1 1]);
imshow(0);

subplot('position',[0.5 0 1 1]);
imshow(1);

repo: OpenTLD


text:
%kImageIcon 'ImageIcon Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ImageIcon.pane file
%
% Parameters:
% Integer: imageiconSize 'Image Icon Size', default: 1: ' '
% InputFile: imageImagefile ' InputFile', required: 'input file'
%
% Example: kImageIcon(imageImagefile, {'imageiconSize',1;'imageImagefile',''})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The Image Icon Visual Object"
% The image icon visual object displays the image as an icon, or "miniature"
% version of the image.
%
% .item "Image Icon Size"
% This integer value specifies the size of the icon in pixels.
% The image icon size should always be used rather than setting width and
% height separately; otherwise, the image may not be proportionally correct.
%
% .item "Input File"
% The file containing the image to be dsiplayed as an image icon.
%
% .end tagged

function varargout = kImageIcon(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kImageIcon(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'imageiconSize', 1;'imageImagefile', '__input'};
maxval={512,0};
minval={-1,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'Integer','InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'ImageIcon" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function A = UpdateArchive(A,N)
% Update the external archive

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Select only the non-dominated solutions in the archive
    A = A(NDSort(A.objs,1)==1);

    %% Sort the solutions in A according to their crowding distances
    [~,rank] = sort(CrowdingDistance(A.objs),'descend');
    A        = A(rank(1:min(N,end)));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function test_cnn_gradients_are_numerically_correct
batch_x = rand(28,28,5);
batch_y = rand(10,5);
cnn.layers = {
    struct('type', 'i') %input layer
    struct('type', 'c', 'outputmaps', 2, 'kernelsize', 5) %convolution layer
    struct('type', 's', 'scale', 2) %sub sampling layer
    struct('type', 'c', 'outputmaps', 2, 'kernelsize', 5) %convolution layer
    struct('type', 's', 'scale', 2) %subsampling layer
};
cnn = cnnsetup(cnn, batch_x, batch_y);

cnn = cnnff(cnn, batch_x);
cnn = cnnbp(cnn, batch_y);
cnnnumgradcheck(cnn, batch_x, batch_y);

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function test_suite = test_expandPolygon
%TESTEXPANDPOLYGON  One-line description here, please.
%   output = testExpandPolygon(input)
%
%   Example
%   testExpandPolygon
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-06-17,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.
% Licensed under the terms of the LGPL, see the file "license.txt"

test_suite = functiontests(localfunctions);

function testSquare(testCase) %#ok<*DEFNU>
p1 = [10 10];
p2 = [20 10];
p3 = [20 20];
p4 = [10 20];
square = [p1;p2;p3;p4];

expanded5 = [5 5;25 5;25 25;5 25];
expanded = expandPolygon(square, 5);
testCase.assertTrue(length(expanded)==1);
testCase.assertTrue(distancePolygons(expanded{1}, expanded5)<1e-14);

function testOverlap(testCase)

% a polygon whose expansion overlaps at critical distance 10
poly = [10 10;190 10;190 140;110 140;110 120;150 120;150 50;50 50;50 100;90 100;90 160;10 160];

% small value: only one outline
expanded = expandPolygon(poly, 5);
testCase.assertTrue(length(expanded)==1);

% value>10: two outlines
expanded = expandPolygon(poly, 20, 'cleanupLoops', true);
testCase.assertTrue(length(expanded)==2);

% value>30: the inner outline disappear
expanded = expandPolygon(poly, 35, 'cleanupLoops', true);
testCase.assertTrue(length(expanded)==1);


text:
function TestVAD(wavRoot, output, step, savestep, DEBUG)
if nargin<3
    step = 1;
end
if nargin<4
    savestep=50;
end
if nargin<5
    DEBUG = 1;
end

files = my_dir([wavRoot '\*.wav']);

for i=1:length(files)
    PrintProgress(i, length(files), step, files{i});
    [wav, Fs] = wavread([wavRoot '\' files{i}]);
    [comb_vad_flag{i}, energy_val{i}, periodicVal{i}] = CombVAD(wav, Fs, DEBUG);
    if mod(i,savestep)==0
        save(output, 'comb_vad_flag', 'energy_val', 'periodicVal', 'files');
    end
end
save(output, 'comb_vad_flag', 'energy_val', 'periodicVal', 'files');

% example of post processing

% double frame rate
currVAD = comb_vad_flag{1};
currVAD = interp1([1:length(currVAD)]*2, currVAD, 1:(2*length(currVAD)), 'nearest', 'extrap');

% Post process VAD with a buffer length of 0.4s
[currVAD_merged, currVAD_extended] = PostProcessVAD(currVAD, 40);

end


text:
function [gx] = g_metaToM(x,P,u,inG)
% observation function for meta-learner (k-ToM vs sequence)
% function [gx] = g_metaToM(x,P,u,inF)
% The meta-learner bases her decision (a=1 or a=0) upon her prediction of
% her opponent's next move, given the game payoff table. The specific
% difficulty of the meta-learner is that she does not know with absolute
% certainty whether she faces an intentional agent (k-ToM) or an inanimate
% patter (sequence). In fact, she holds a belief about this, which is
% quantified in terms of the probability Pi that she is facing an
% intentional agent (the probability that the agent faces an inanimate
% pattern is 1-Pi).
% IN:
%   - x: hidden states (see indexing in inF.indlev)
%   - P: observation params: phi(1) = log-temperature and phi(2) = bias
%   - u: sequence of past actions:
%       u(1)= opponent's last move
%       u(2)= learner's last move
%       u(3:K+2) = sequence of K past opponent's moves
%   - inG: input structure (see prepare_metaToM.m)
% OUT:
%   - gx: updated hidden states
% [see RecToMfunction.m and f_BSL.m]

% 1- get k-ToM probabilistic decision P(a=1|k-ToM)
xktom = x(inG.ktom.indx);
inG_ktom = inG.ktom.inG;
gx_ktom = ObsRecGen(xktom,P,u,inG_ktom);

% 2- get BSL probabilistic decision P(a=1|seq)
xseq = x(inG.seq.indx);
inG_seq = inG.seq.inG;
gx_bsl = g_BSLinGame(xseq,P,u,inG_seq);

% 3- derive meta-ToM probabilistic decision P(a=1)
Pi = VBA_sigmoid(x(inG.meta.indx)); % prior P(agent=kToM)
gx = Pi*gx_ktom + (1-Pi)*gx_bsl;

author: MBB-team

repo: VBA-toolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	Implemetation of the tracker described in paper
%	"MEEM: Robust Tracking via Multiple Experts using Entropy Minimization",
%   Jianming Zhang, Shugao Ma, Stan Sclaroff, ECCV, 2014
%
%	Copyright (C) 2014 Jianming Zhang
%
%	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/>.
%
%	If you have problems about this software, please contact: jmzhang@bu.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function config = makeConfig(frame,selected_rect,use_color,use_experts,use_iif, show_img)
warning('off','MATLAB:maxNumCompThreads:Deprecated');
maxNumCompThreads(1);% The rounding error due to using differen number of threads
                     % could cause different tracking results.  On most sequences,
                     % the difference is very small, while on some challenging
                     % sequences, the difference can be substantial due to
                     % "butterfly effects". Therefore, we suggest using
                     % Spatial Robustness Evaluation (SRE) to benchmark
                     % trackers.
if nargin < 5
    use_iif = true;
end

%% default setting up

config.search_roi = 2; % ratio of the search roi to tracking window
config.padding = 40; % for object out of border

config.debug = false;
config.verbose = false;
config.display = show_img; % show tracking result at runtime
config.use_experts = use_experts;
config.use_raw_feat = false; % raw intensity feature value
config.use_iif = use_iif; % use illumination invariant feature

config.svm_thresh = -0.7; % for detecting the tracking failure
config.max_expert_sz = 4;
config.expert_update_interval = 50;
config.update_count_thresh = 1;
config.entropy_score_winsize = 5;
config.expert_lambda = 10;
config.label_prior_sigma = 15;

config.hist_nbin = 32; % histogram bins for iif computation

config.thresh_p = 0.1; % IOU threshold for positive training samples
config.thresh_n = 0.5; % IOU threshold for negative ones

%%
%  automatic setting up for determining grid sample step and feature channels
%  (do not change the following)

% check if the frame is in RGB format
config.use_color = false;
if (size(frame,3) == 3 && ~isequal(frame(:,:,1),frame(:,:,2),frame(:,:,3))) && use_color
    config.use_color = true;
end

% decide feature channel number
if config.use_color
    thr_n = 5;
else
    thr_n = 9;
end
config.thr = (1/thr_n:1/thr_n:1-1/thr_n)*255;
config.fd = numel(config.thr);

% decide image scale and pixel step for sampling feature
% rescale raw input frames propoerly would save much computation
frame_min_width = 320;
trackwin_max_dimension = 64;
template_max_numel = 144;
frame_sz = size(frame);

if max(selected_rect(3:4)) <= trackwin_max_dimension ||...
        frame_sz(2) <= frame_min_width
    config.image_scale = 1;
else
    min_scale = frame_min_width/frame_sz(2);
    config.image_scale = max(trackwin_max_dimension/max(selected_rect(3:4)),min_scale);
end
wh_rescale = selected_rect(3:4)*config.image_scale;
win_area = prod(wh_rescale);
config.ratio = (sqrt(template_max_numel/win_area));
template_sz = round(wh_rescale*config.ratio);
config.template_sz = template_sz([2 1]);

author: flyers

repo: drone-tracking


text:
function [new_mu, new_Sigma, new_Sigma2] = collapse_mog(mu, Sigma, coefs)
% COLLAPSE_MOG Collapse a mixture of Gaussians to a single Gaussian by moment matching
% [new_mu, new_Sigma] = collapse_mog(mu, Sigma, coefs)
%
% coefs(i) - weight of i'th mixture component
% mu(:,i), Sigma(:,:,i) - params of i'th mixture component

% S = sum_c w_c (S_c + m_c m_c' + m m' - 2 m_c m')
%   = sum_c w_c (S_c + m_c m_c') + m m' - 2 (sum_c m_c) m'
%   = sum_c w_c (S_c + m_c m_c') - m m'

new_mu = sum(mu * diag(coefs), 2); % weighted sum of columns

n = length(new_mu);
new_Sigma = zeros(n,n);
new_Sigma2 = zeros(n,n);
for j=1:length(coefs)
  m = mu(:,j) - new_mu;
  new_Sigma = new_Sigma + coefs(j) * (Sigma(:,:,j) + m*m');
  new_Sigma2 = new_Sigma2 + coefs(j) * (Sigma(:,:,j) + mu(:,j)*mu(:,j)');
end
%assert(approxeq(new_Sigma, new_Sigma2 - new_mu*new_mu'))

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [A,b,Aeq,beq]=addBounds(A0,b0,Aeq0,beq0,lb,ub)
%addBounds - given a sextuplet {A, b, Aeq, beq, lb, ub} of matrices and vectors
%of the kind used by the Optimization Toolbox to express linear
%constraints, the routine will convert the sextuplet to a quadruplet {A, b,
%Aeq, beq}.
%
%      [A,b,Aeq,beq]=addBounds(A,b,Aeq,beq,lb,ub)
%
%In other words, it will re-write the bounds given by vectors lb, ub as equalities
%and inequalities C*x<=c, D*x=d and append these as new rows to the matrices
%A, b, Aeq, beq.
%
%
%EXAMPLE:
%
%   >> [A,b,Aeq,beq]=addBounds([1 1 1],1,[],[],[0;0;1],[1;1;1])
%
%         A =
%
%              1     1     1
%             -1     0     0
%              0    -1     0
%              1     0     0
%              0     1     0
%
%
%         b =
%
%              1
%              0
%              0
%              1
%              1
%
%
%         Aeq =
%
%              0     0     1
%
%
%         beq =
%
%              1

  %%%%begin parsing

   if nargin<5, error 'At least 5 arguments required'; end

   if size(A0,1)~=length(b0)
       error 'Incompatible inequality matrix data sizes: size(A,1) ~= length(b)'
   end

   if size(Aeq0,1)~=length(beq0)
       error 'Incompatible equality matrix data sizes: size(Aeq,1) ~= length(beq)'
   end

     if ~exist('lb','var'), lb=[]; end
      if ~exist('ub','var'), ub=[]; end

      if ~isempty(lb)
         N=length(lb);
      elseif ~isempty(ub);
         N=length(ub);
      else
          error 'Either lb or ub must be nonempty.'
      end

      if isempty(lb)
        lb=-inf(N,1);
      end

      if isempty(ub)
        ub=+inf(N,1);
      end

      if length(ub)~=length(lb)
         error 'Arguments lb and ub must be either [] or the same lengths'
      end

      lb=lb(:); ub=ub(:);

      %%%end parsing

      lisinf=~isfinite(lb);
      uisinf=~isfinite(ub);
      idx_eq=(ub==lb)& ~lisinf & ~uisinf;

      Au=[eye(N),ub];
      Al=[-eye(N),-lb];
      Aeq=Au(idx_eq,:);
      Au(idx_eq|uisinf,:)=[];
      Al(idx_eq|lisinf,:)=[];

      A=[Al;Au]; b=A(:,end); A(:,end)=[];
                 beq=Aeq(:,end); Aeq(:,end)=[];

      A=[A0;A];
      b=[b0;b];
      Aeq=[Aeq0;Aeq];
      beq=[beq0;beq];

      if any(lb==inf | ub==-inf)
         A(end+1,end)=0; b(end+1)=-1;
      end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [nn, L]  = nntrain(nn, train_x, train_y, opts, val_x, val_y)
%NNTRAIN trains a neural net
% [nn, L] = nnff(nn, x, y, opts) trains the neural network nn with input x and
% output y for opts.numepochs epochs, with minibatches of size
% opts.batchsize. Returns a neural network nn with updated activations,
% errors, weights and biases, (nn.a, nn.e, nn.W, nn.b) and L, the sum
% squared error for each training minibatch.

assert(isfloat(train_x), 'train_x must be a float');
assert(nargin == 4 || nargin == 6,'number ofinput arguments must be 4 or 6')

loss.train.e               = [];
loss.train.e_frac          = [];
loss.val.e                 = [];
loss.val.e_frac            = [];
opts.validation = 0;
if nargin == 6
    opts.validation = 1;
end

fhandle = [];
if isfield(opts,'plot') && opts.plot == 1
    fhandle = figure();
end

m = size(train_x, 1);

batchsize = opts.batchsize;
numepochs = opts.numepochs;

numbatches = m / batchsize;
numbatches

assert(rem(numbatches, 1) == 0, 'numbatches must be a integer');

L = zeros(numepochs*numbatches,1);
n = 1;
for i = 1 : numepochs
    tic;

    kk = randperm(m);
    for l = 1 : numbatches
        batch_x = train_x(kk((l - 1) * batchsize + 1 : l * batchsize), :);

        %Add noise to input (for use in denoising autoencoder)
        if(nn.inputZeroMaskedFraction ~= 0)
            batch_x = batch_x.*(rand(size(batch_x))>nn.inputZeroMaskedFraction);
        end

        batch_y = train_y(kk((l - 1) * batchsize + 1 : l * batchsize), :);

        nn = nnff(nn, batch_x, batch_y);
        nn = nnbp(nn);
        nn = nnapplygrads(nn);

        L(n) = nn.L;

        n = n + 1;
    end

    t = toc;

    if opts.validation == 1
        loss = nneval(nn, loss, train_x, train_y, val_x, val_y);
        str_perf = sprintf('; Full-batch train mse = %f, val mse = %f', loss.train.e(end), loss.val.e(end));
    else
        loss = nneval(nn, loss, train_x, train_y);
        str_perf = sprintf('; Full-batch train err = %f', loss.train.e(end));
    end
    if ishandle(fhandle)
        nnupdatefigures(nn, fhandle, loss, opts, i);
    end

    disp(['epoch ' num2str(i) '/' num2str(opts.numepochs) '. Took ' num2str(t) ' seconds' '. Mini-batch mean squared error on training set is ' num2str(mean(L((n-numbatches):(n-1)))) str_perf]);
    nn.learningRate = nn.learningRate * nn.scaling_learningRate;
end
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment2B_K1_N30
close all;
global parameters

%STOMP PARAMETERS
%number of particles
K = 1;

%repeat experiment number of times
parameters.n_repeat = 30;
parameters.experiment_name = 'experiment2B_K1_N30.mat';
parameters.animate = 0;

%repeat the experiment E times
random_manips=[];
Gout = [];
for i=1:parameters.n_repeat
    close all
    [pk, final_manip] = path_planning_SCO_4DOF(K);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    save(parameters.experiment_name)
end

'ended'
parameters.experiment_name

repo: ARTE


text:
function [I,NextSigma,NextdeltaMu] = VarVolatility(x,mu,va,a,b)
% OTO: varational energy (and curvature) of volatile observer

if ~isinf(va)
    iva = va.^-1;
else
    iva = 0;
end
ex = exp(x(:));
v = a + ex;
I = -0.5.*iva*(x(:)-mu).^2 ...
    -0.5*log(v) ...
    -0.5*b./v;
dIdx = iva*(mu-x(:)) - 0.5.*ex./v + 0.5.*b.*ex./v.^2;
dI2dx2 = -iva -0.5.*a*ex./v.^2 -0.5.*b.*ex.*(ex-a)./v.^3;
NextSigma = -dI2dx2.^-1;
NextdeltaMu = NextSigma.*dIdx;


text:
function f = p04_f ( n, x )

%*****************************************************************************80
%
%% P04_F returns the integrand for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the arguments.
%
%    Output, real F(N), the values of the integrand.
%
  f = 0.75 * x + 1.25 * sin ( 3.0 * pi * x + 0.4 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_mrvConvertVAnat
%Validate that mrVista conversopm of vAanatomy.data to nifti works
%
%  test_mrvConvertVAnat()
%
% Tests: mrAnatConvertVAnatToT1Nifti
%
% INPUTS
%  No inputs
%
% RETURNS
%  No returns
%
% Example: test_mrvConvertVAnat()
%
% See also MRVTEST
%
% Copyright NYU team, mrVista, 2017

%% Convert vAnatomy.dat to nifti

% Get the vAnatomy.dat sample data file
vAnatFileName = mrtInstallSampleData('anatomy/anatomyV','vAnatomy',  [], [], 'filetype', 'dat');

% Convert it to nifti
newNiftiFileName = fullfile(tempdir,'vAnatomy.nii.gz');
mrAnatConvertVAnatToT1Nifti(vAnatFileName,newNiftiFileName);

% Get a stored nifti file that corresponds to the vAnatomy
oldNifitFileName = mrtInstallSampleData('anatomy/anatomyV','t1.nii', [], [], 'filetype', 'gz');

%% Compare the nifti to some stored values
n0  = niftiRead(oldNifitFileName);
n1 = niftiRead(newNiftiFileName);

n0  = niftiApplyCannonicalXform(n0);
n1 = niftiApplyCannonicalXform(n1);

assertEqual(niftiGet(n0, 'qto_xyz'), niftiGet(n1, 'qto_xyz'));
assertEqual(niftiGet(n0, 'data'), niftiGet(n1, 'data'));
assertEqual(niftiGet(n0, 'dim'), niftiGet(n1, 'dim'));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [SEF, TP] = spectral_edge(PSD,freq,f1,f2,percent)
i1 = find(freq==f1);    % f1 hertz index
i2 = find(freq==f2);    % f2 Hz index

P = PSD(i1:i2,:);       % Power in range f1-f2

edgeTP = zeros(1,size(P,2));
SEF = zeros(1,size(P,2));
TP = sum(P,1);

for j=1:size(P,2)
    i=1;
    edgeTP(j) = percent*sum(P(:,j));

    while((sum(P([1:i],j)) < edgeTP(j)) & (i <= (i2-i1)))
        i=i+1;
    end

    SEF(j) = freq(i1+i);       %spectral edge freq
end

author: alexandrebarachant


text:
function y = blkdiag( varargin )

%Disciplined convex/geometric programming information for BLKDIAG:
%   BLKDIAG imposes no convexity restrictions on its arguments.

nv = 0;
nz = 0;
sz = [ 0, 0 ];
for k = 1 : nargin,
    x  = cvx( varargin{k} );
    sx = x.size_;
    if length( sx ) > 2,
        error( 'N-D matrices not supported.' );
    end
    b  = x.basis_;
    sz = sz + sx;
    nv = max( nv, size( b, 1 ) );
    nz = nz + nnz( b );
    varargin{k} = x;
end
bz = sparse( [], [], [], prod( sz ), nz, nv );
roff = 0;
coff = 0;
for k = 1 : nargin,
    x  = varargin{k};
    b  = x.basis_;
    sx = x.size_;
    ndxr = ( roff : roff + sx( 1 ) - 1 )';
    ndxr = ndxr( :, ones( 1, sx( 2 ) ) );
    ndxc = ( coff : coff + sx( 2 ) - 1 );
    ndxc = ndxc( ones( 1, sx( 1 ) ), : );
    bz( 1 : size( b, 1 ), ndxc( : ) * sz( 1 ) + ndxr( : ) + 1 ) = b; %#ok
    roff = roff + sx( 1 );
    coff = coff + sx( 2 );
end
y = cvx( sz, bz );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function x = p11_start ( option, nvar )

%*****************************************************************************80
%
%% P11_START returns a starting point for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Output, real X(NVAR), the starting point.
%
  x = zeros ( nvar, 1 );

  return
end


text:
function z = zernpol(n,m,r,nflag)
%ZERNPOL Radial Zernike polynomials of order N and frequency M.
%   Z = ZERNPOL(N,M,R) returns the radial Zernike polynomials of
%   order N and frequency M, evaluated at R.  N is a vector of
%   positive integers (including 0), and M is a vector with the
%   same number of elements as N.  Each element k of M must be a
%   positive integer, with possible values M(k) = 0,2,4,...,N(k)
%   for N(k) even, and M(k) = 1,3,5,...,N(k) for N(k) odd.  R is
%   a vector of numbers between 0 and 1.  The output Z is a matrix
%   with one column for every (N,M) pair, and one row for every
%   element in R.
%
%   Z = ZERNPOL(N,M,R,'norm') returns the normalized Zernike poly-
%   nomials.  The normalization factor Nnm = sqrt(2*(n+1)) is
%   chosen so that the integral of (r * [Znm(r)]^2) from r=0 to
%   r=1 is unity.  For the non-normalized polynomials, Znm(r=1)=1
%   for all [n,m].
%
%   The radial Zernike polynomials are the radial portion of the
%   Zernike functions, which are an orthogonal basis on the unit
%   circle.  The series representation of the radial Zernike
%   polynomials is
%
%          (n-m)/2
%            __
%    m      \       s                                          n-2s
%   Z(r) =  /__ (-1)  [(n-s)!/(s!((n-m)/2-s)!((n+m)/2-s)!)] * r
%    n      s=0
%
%   The following table shows the first 12 polynomials.
%
%       n    m    Zernike polynomial    Normalization
%       ---------------------------------------------
%       0    0    1                        sqrt(2)
%       1    1    r                           2
%       2    0    2*r^2 - 1                sqrt(6)
%       2    2    r^2                      sqrt(6)
%       3    1    3*r^3 - 2*r              sqrt(8)
%       3    3    r^3                      sqrt(8)
%       4    0    6*r^4 - 6*r^2 + 1        sqrt(10)
%       4    2    4*r^4 - 3*r^2            sqrt(10)
%       4    4    r^4                      sqrt(10)
%       5    1    10*r^5 - 12*r^3 + 3*r    sqrt(12)
%       5    3    5*r^5 - 4*r^3            sqrt(12)
%       5    5    r^5                      sqrt(12)
%       ---------------------------------------------
%
%   Example:
%
%       % Display three example Zernike radial polynomials
%       r = 0:0.01:1;
%       n = [3 2 5];
%       m = [1 2 1];
%       z = zernpol(n,m,r);
%       figure
%       plot(r,z)
%       grid on
%       legend('Z_3^1(r)','Z_2^2(r)','Z_5^1(r)','Location','NorthWest')
%
%   See also ZERNFUN, ZERNFUN2.

% A note on the algorithm.
% ------------------------
% The radial Zernike polynomials are computed using the series
% representation shown in the Help section above. For many special
% functions, direct evaluation using the series representation can
% produce poor numerical results (floating point errors), because
% the summation often involves computing small differences between
% large successive terms in the series. (In such cases, the functions
% are often evaluated using alternative methods such as recurrence
% relations: see the Legendre functions, for example). For the Zernike
% polynomials, however, this problem does not arise, because the
% polynomials are evaluated over the finite domain r = (0,1), and
% because the coefficients for a given polynomial are generally all
% of similar magnitude.
%
% ZERNPOL has been written using a vectorized implementation: multiple
% Zernike polynomials can be computed (i.e., multiple sets of [N,M]
% values can be passed as inputs) for a vector of points R.  To achieve
% this vectorization most efficiently, the algorithm in ZERNPOL
% involves pre-determining all the powers p of R that are required to
% compute the outputs, and then compiling the {R^p} into a single
% matrix.  This avoids any redundant computation of the R^p, and
% minimizes the sizes of certain intermediate variables.
%
%   Paul Fricker 11/13/2006

% Check and prepare the inputs:
% -----------------------------
if ( ~any(size(n)==1) ) || ( ~any(size(m)==1) )
    error('zernpol:NMvectors','N and M must be vectors.')
end

if length(n)~=length(m)
    error('zernpol:NMlength','N and M must be the same length.')
end

n = n(:);
m = m(:);
length_n = length(n);

if any(mod(n-m,2))
    error('zernpol:NMmultiplesof2','All N and M must differ by multiples of 2 (including 0).')
end

if any(m<0)
    error('zernpol:Mpositive','All M must be positive.')
end

if any(m>n)
    error('zernpol:MlessthanN','Each M must be less than or equal to its corresponding N.')
end

if any( r>1 | r<0 )
    error('zernpol:Rlessthan1','All R must be between 0 and 1.')
end

if ~any(size(r)==1)
    error('zernpol:Rvector','R must be a vector.')
end

r = r(:);
length_r = length(r);

if nargin==4
    isnorm = ischar(nflag) & strcmpi(nflag,'norm');
    if ~isnorm
        error('zernpol:normalization','Unrecognized normalization flag.')
    end
else
    isnorm = false;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the Zernike Polynomials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Determine the required powers of r:
% -----------------------------------
rpowers = [];
for j = 1:length(n)
    rpowers = [rpowers m(j):2:n(j)];
end
rpowers = unique(rpowers);

% Pre-compute the values of r raised to the required powers,
% and compile them in a matrix:
% -----------------------------
if rpowers(1)==0
    rpowern = arrayfun(@(p)r.^p,rpowers(2:end),'UniformOutput',false);
    rpowern = cat(2,rpowern{:});
    rpowern = [ones(length_r,1) rpowern];
else
    rpowern = arrayfun(@(p)r.^p,rpowers,'UniformOutput',false);
    rpowern = cat(2,rpowern{:});
end

% Compute the values of the polynomials:
% --------------------------------------
z = zeros(length_r,length_n);
for j = 1:length_n
    s = 0:(n(j)-m(j))/2;
    pows = n(j):-2:m(j);
    for k = length(s):-1:1
        p = (1-2*mod(s(k),2))* ...
                   prod(2:(n(j)-s(k)))/          ...
                   prod(2:s(k))/                 ...
                   prod(2:((n(j)-m(j))/2-s(k)))/ ...
                   prod(2:((n(j)+m(j))/2-s(k)));
        idx = (pows(k)==rpowers);
        z(:,j) = z(:,j) + p*rpowern(:,idx);
    end

    if isnorm
        z(:,j) = z(:,j)*sqrt(2*(n(j)+1));
    end
end

% EOF zernpol

author: goGPS-Project


text:
function [Data_TV,avg_votenum] = Reversed_Sparse_Ball_faster(data, Sigma, Sigma_noise, method, index_mat_vote, S)

if strcmp( method, 'Std2' ) == 1 % Sparse Standard Ball Tensor Voting with my code

    disp('Std2');
    disp('faster version');
    D = size( data, 2 );
    N = size( data, 1 );
    votenum = 0;
    data_inv=(data(:, 1:D))';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    F = S<sqrt(Sigma./3);
    Data_TV(:,1:D) = data(:,1:D);
    size(F)
    for i = 1 : size( data, 1 )
       Data_TV(i,D+1) = 1;
       Tensor = zeros(D,D); % 2nd order tensor - matrix
       temp_nei_i=index_mat_vote(:,i);
       length(temp_nei_i);
       l=length(temp_nei_i);
        for j = 1:l
        idxNeig = temp_nei_i(j);
       % idxNeig = (j);
        dis = S(j,i);
           % if (dis > 0) & (F(i,idxNeig))
         if (dis > 0) & (F(j,i))
                votenum = votenum + 1;
                dir = ( data(i,1:D) - data(idxNeig,1:D) ) ;
                Tensor = Tensor + exp( -dis^2/ Sigma ) * (eye(D) - dir'*dir);

            end
        end

       temp_dim=1;
       [U1,S1,V1] = rsvd(Tensor,temp_dim);
 %Extract first temp_dim eigenvlaues and k eigenvectors:

 %disp('fast SVD computing');

 if  isempty(S1)==1
     S1 = zeros(temp_dim,D);
     U1 = zeros(D,temp_dim);
     V1 = zeros(D,temp_dim);
 end
      Lam_k=S1(1:temp_dim,:);
      Vec_k=U1*V1';

% %         size(U1)
      % [L,J] = sort( diag(Lam) );
        [L,J] = sort( diag(Lam_k) );
        for k = 1 : temp_dim
        %Data_TV(i,D+1+k) =  Lam(J(D-k+1),J(D-k+1));% here the ranking is very important, a slight error before but revised now^_^
        Data_TV(i,D+1+k) =  Lam_k(J(temp_dim-k+1),J(temp_dim-k+1));
        end
        Data_TV(i, D+1+temp_dim+1 : (2*D)+1 )= 0;
        for k = 1 : temp_dim % instead of D
        %Data_TV( i, 2*D+1+(k-1)*D+1 : 2*D+1+k*D ) = Vec(:,J(D-k+1));
        Data_TV( i, 2*D+1+(k-1)*D+1 : 2*D+1+k*D ) = Vec_k(:,J(temp_dim-k+1));
        end
      %   error_vec(i)=error;
    end
    avg_votenum=votenum/N;
    disp('Average Vote is'); disp(votenum/N);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [feat_prob, locreg_pred, next_pred, rpn_prob, rpn_bbox] = cnn_process_image(input, net, sigmoid)

    % switch width and height for Caffe
    input = permute(input, [2 1 3]);

    blob_size = [size(input), 1];
    net.blobs('data').reshape(blob_size);

    blob_size = [size(input), 1];

    batch = zeros(blob_size, 'single');
    batch(:,:,:,1) = input;
    batches = cell(1,1);
    batches{1} = batch;
    zs = net.forward(batches);

    prob_idx = find_string(net.outputs, 'prob');
    locref_idx = find_string(net.outputs, 'loc_pred');
    next_idx = find_string(net.outputs, 'next_pred');
    rpn_prob_idx = find_string(net.outputs, 'rpn_prob');
    rpn_bbox_idx = find_string(net.outputs, 'rpn_bbox_pred');

    % fc6_data = net.blobs('fc6-conv').get_data();

    locreg_pred = [];
    next_pred = [];
    rpn_prob = [];
    rpn_bbox = [];

    if prob_idx ~= -1
        feat_prob = zs{prob_idx};
        feat_prob = permute(feat_prob, [2 1 3]);
        if ~sigmoid
            % remove background class
            feat_prob = feat_prob(:,:,2:end);
        end
    end

    if locref_idx ~= -1
        fmap = zs{locref_idx};
        sz = size(fmap);
        fmap = reshape(fmap, sz(1), sz(2), 2, []);
        fmap = permute(fmap, [2 1 4 3]);
        locreg_pred = fmap;
    end

    if next_idx ~= -1
        fmap = zs{next_idx};
        sz = size(fmap);
        fmap = reshape(fmap, sz(1), sz(2), 2, []);
        fmap = permute(fmap, [2 1 4 3]);
        next_pred = fmap;
    end

    if rpn_prob_idx ~= -1
        fmap = zs{rpn_prob_idx};
        fmap = permute(fmap, [2 1 3]);
        rpn_prob = fmap;
    end

    if rpn_bbox_idx ~= -1
        fmap = zs{rpn_bbox_idx};
        sz = size(fmap);
        fmap = reshape(fmap, sz(1), sz(2), 4, []);
        fmap = permute(fmap, [2 1 4 3]);
        rpn_bbox = fmap;
    end
end

function idx = find_string(list, str)
    for k = 1:length(list)
        cand = list{k};
        if ~isempty(strfind(cand, str))
            idx = k;
            return;
        end
    end
    idx = -1;
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [fmodel,removed] = filterMetMNXIDsViaRxns(model,mnx,ignoreComp,keepAtLeastOne)
%filterMetMNXIDsViaRxns  Remove met MNXID associations based on their rxns.
%
% filterMetMNXIDsViaRxns determines which met MNXIDs (if any) should be
% removed for mets associated with multiple MNXIDs. This is accomplished by
% finding all model rxns that involve the met, obtaining their associated
% MNX rxn IDs, and retrieving those rxns from the MNX database. Any MNXIDs
% associated with the met that are not included in that set of rxns from
% the MNX database will be removed.
%
% NOTE: metabolites that occur only in reactions that have no associated
%       rxnMNXIDs will be skipped, as there is insufficient evidence to
%       properly filter their associated metMNXIDs.
%
%
% USAGE:
%
%   [fmodel,removed] = filterMetMNXIDsViaRxns(model,mnx,ignoreComp,keepAtLeastOne);
%
% INPUT:
%
%   model    A model structure containing reaction and metabolite MNX ID
%            association fields ("rxnMNXID" and "metMNXID", respectively).
%
%   mnx      (Optional) An MNX database structure, containing reaction-
%            related information retrieved from the MNX database, generated
%            using the following command: mnx = buildMNXmodel('rxn');
%            By default, the function will automatically run the above
%            command to regenerate the MNX database structure (slower).
%
%   ignoreComp   (Optional, Default FALSE) If TRUE, metabolite compartments
%                will be ignored. In this case, identical mets of different
%                compartments will be lumped together, and when searching
%                for rxns involving the metabolite, the compartment will be
%                ignored.
%
%   keepAtLeastOne  (Optional, Default FALSE) In some cases, none of the
%                   MNXIDs associated with a metabolite are found in any of
%                   the reactions, and will result in a met with zero MNXID
%                   associations in the filtered model.
%                   If keepAtLeastOne = TRUE, then in cases such as this,
%                   the MNXID associations will not be removed in the
%                   filtered model, but will be indicated in the "removed"
%                   structure, with the text "ALL SHOULD BE REMOVED!" next
%                   to the metabolite ID.
%
% OUTPUT:
%
%   fmodel   A filtered model, which has removed metabolite MNXIDs that
%            were not found in any of the MNX reactions associated with the
%            model and that metabolite.
%
%   removed  A structure containing more detailed information on the met
%            MNXID associations that were removed.
%

% handle input arguments
if nargin < 2 || isempty(mnx)
    mnx = buildMNXmodel('rxn');
end
if nargin < 3
    ignoreComp = false;
end
if nargin < 4
    keepAtLeastOne = false;
end

% initialize "removed" structure and "fmodel" outputs
fmodel = model;
removed.mets = {};
removed.metMNXID = {};

% if metMNXID and/or rxnMNXID field contains multiple columns, consolidate
% into a single column with nested cell entries
if size(model.metMNXID,2) > 1
    model.metMNXID = nestCell(model.metMNXID,true);
end
if size(model.rxnMNXID,2) > 1
    model.rxnMNXID = nestCell(model.rxnMNXID,true);
end

if ( ignoreComp )

    % check if model has already merged met compartments, or if "mets"
    % contains non-unique elements
    if length(unique(model.mets)) ~= length(model.mets)
        error('Input model "mets" field must contain unique (non-repeated) elements.');
    elseif any(regexp(model.mets{1},'\d$'))
        % Note: this test is specific to models whose met IDs end in
        % numbers (e.g., "m00001"), with compartments appended to the end
        % of the ID (e.g., "m00001c" or "m00001[c]").
        fprintf('\nIt appears that the model compartments have already been merged.\n');
        fprintf('The "ignoreComp" flag will be ignored, since it will have no effect.\n');
        S = model.S;
        ignoreComp = false;
    else

        % strip compartment label from model met ID
        if endsWith(model.mets{1},']')
            % compartment name is formatted as "m00001[c]"
            model.mets = regexprep(model.mets,'\[\w\]$','');
        else
            % compartment name is formatted as "m00001c"
            model.mets = regexprep(model.mets,'.$','');
        end

        % check if ignoring compartments will actually do anything
        if length(unique(model.mets)) == length(model.mets)
            fprintf('\nIt appears that the model compartments have already been merged.\n');
            fprintf('The "ignoreComp" flag will be ignored, since it will have no effect.\n');
            S = model.S;
            ignoreComp = false;
        else

            fprintf('\nMetabolite compartments will be ignored. NOTE: This process will merge metMNXIDs\n');
            fprintf('of mets that are identical except for their compartment. If any mets are associated\n');
            fprintf('with compartment-specific metMNXIDs, this is not recommended.\n');

            % If ignoring compartments, convert the stoich matrix into a
            % binary met-rxn association matrix (i.e., set all nonzero
            % entries = 1), and combine all associations for metabolites
            % that are identical except for their compartment. Also combine
            % their metMNXIDs.
            S = (model.S ~= 0);
            [~,uniq_ind,met_groups] = unique(model.mets);
            h = waitbar(0,'Merging model metabolites across compartments...');
            for i = 1:max(met_groups)
                ind = (met_groups == i);
                S(ind,:) = repmat(any(S(ind,:),1),sum(ind),1);
                model.metMNXID(ind) = repmat({unique(horzcat(model.metMNXID{ind}))},sum(ind),1);
                waitbar(i/max(met_groups),h);
            end
            close(h);

            % now merge mets of different compartments into single met
            model.mets_nocomp = model.mets;  % first save these for indexing later
            model.mets = model.mets(uniq_ind);
            S = S(uniq_ind,:);
            model.metMNXID = model.metMNXID(uniq_ind);

        end

    end

else
    S = model.S;
end

% % Identify all mets with two or more MNXID associations. Mets with one or
% % zero MNXIDs will be ignored.
% multi_ind = find(cellfun(@numel,model.metMNXID) > 1);

% Iterate through mets with multiple MNXIDs, and determine which MNXIDs
% should be removed.
h = waitbar(0,'Processing metabolites...');
for i = 1:length(model.mets)

    % get indices of all rxns in which the met participates
    rxn_ind = S(i,:) ~= 0;

    % obtain a unique list of all the rxn MNXIDs associated with those
    % rxns, as well as their indices in the MNX database structure
    MNX_rxnIDs = unique(horzcat(model.rxnMNXID{rxn_ind}));
    MNX_rxn_ind = ismember(mnx.rxns,MNX_rxnIDs);

    if isempty(MNX_rxnIDs)
        % If none of the rxns are associated to any rxn MNXIDs, just skip
        % this metabolite. Otherwise, all met MNXIDs associated with the
        % metabolite will be removed, which is not so helpful.
        continue
    end

    % obtain unique set of all mets participating in the MNX rxns
    MNX_metIDs = unique([mnx.rxnMets{MNX_rxn_ind}]);

    % determine which of the MNXIDs currently associated to the model met
    % are not found in any of the MNX rxns
    rem_mets = ~ismember(model.metMNXID{i},MNX_metIDs);

    if ~any(rem_mets)
        % if no met MNXID associations should be removed, continue to next met
        continue
    elseif all(rem_mets) && (keepAtLeastOne)
        % In this case, all the MNXIDs associated with this met are to be
        % removed, but the user has indicated that they do not want this to
        % happen. Therefore,  NONE will be removed, but the information
        % will be reported in the "removed" structure, so these cases can
        % be manually evaluated by the user.
        removed.mets = [removed.mets; model.mets(i)];
        removed.metMNXID = [removed.metMNXID; {'ALL SHOULD BE REMOVED!'}];
    else
        % remove one or more MNXID associations from the metabolite
        removed.mets = [removed.mets; model.mets(i)];
        removed.metMNXID = [removed.metMNXID; {model.metMNXID{i}(rem_mets)}];
        model.metMNXID{i}(rem_mets) = [];
    end

    waitbar(i/length(model.mets),h);
end
close(h);

if any(cellfun(@(x) ismember({'ALL SHOULD BE REMOVED!'},x),removed.metMNXID))
    fprintf('\n*** NOTE: There was at least one case where ALL the MNXIDs associated\n');
    fprintf('          with a metabolite were not found in any of the involved rxns,\n');
    fprintf('          and therefore NONE were removed. See "removed" structure for\n');
    fprintf('          further information.\n\n');
end

% update filtered model (fmodel) metMNXID field
if ( ignoreComp )
    % if compartment was ignored (and thus mets were merged), distribute
    % the updated metMNXID assignments to all compartment-versions of each
    % metabolite.
    for i = 1:length(model.mets)
        ind = ismember(model.mets_nocomp,model.mets(i));
        fmodel.metMNXID(ind) = model.metMNXID(i);
    end
else
    fmodel.metMNXID = model.metMNXID;
end

author: SysBioChalmers

repo: Human-GEM


text:
function [ num_int, pint ] = plane_normal_triangle_int_3d ( pp, normal, t )

%*****************************************************************************80
%
%% PLANE_NORMAL_TRIANGLE_INT_3D: intersection ( normal plane, triangle ) in 3D.
%
%  Discussion:
%
%    The normal form of a plane in 3D is:
%
%      PP is a point on the plane,
%      N is a normal vector to the plane.
%
%    There may be 0, 1, 2 or 3 points of intersection returned.
%
%    If two intersection points are returned, then the entire line
%    between them comprises points of intersection.
%
%    If three intersection points are returned, then all points of
%    the triangle intersect the plane.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PP(3), a point on the plane.
%
%    Input, real NORMAL(3), a normal vector to the plane.
%
%    Input, real T(3,3), the vertices of the triangle.
%
%    Output, integer NUM_INT, the number of intersection points returned.
%
%    Output, real PINT(3,3), the coordinates of the
%    intersection points.
%
  dim_num = 3;

  num_int = 0;
  pint = [];
%
%  Compute the signed distances between the vertices and the plane.
%
  d = - normal(1:dim_num) * pp(1:dim_num)';

  dist1 = normal(1:dim_num) * t(1:dim_num,1) + d;
  dist2 = normal(1:dim_num) * t(1:dim_num,2) + d;
  dist3 = normal(1:dim_num) * t(1:dim_num,3) + d;
%
%  Consider any zero distances.
%
  if ( dist1 == 0.0 )

    num_int = num_int + 1;
    pint(1:dim_num,num_int) = t(1:dim_num,1);

  end

  if ( dist2 == 0.0 )

    num_int = num_int + 1;
    pint(1:dim_num,num_int) = t(1:dim_num,2);

  end

  if ( dist3 == 0.0 )

    num_int = num_int + 1;
    pint(1:dim_num,num_int) = t(1:dim_num,3);

  end
%
%  If 2 or 3 of the nodes intersect, we're already done.
%
  if ( 2 <= num_int )
    return
  end
%
%  If one node intersects, then we're done unless the other two
%  are of opposite signs.
%
  if ( num_int == 1 )

    if ( dist1 == 0.0 )

      [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,2), ...
        t(1:dim_num,3), dist2, dist3, num_int, pint );

    elseif ( dist2 == 0.0 )

      [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), ...
        t(1:dim_num,3), dist1, dist3, num_int, pint );

    elseif ( dist3 == 0.0 )

      [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), ...
        t(1:dim_num,2), dist1, dist2, num_int, pint );

    end

    return

  end
%
%  All nodal distances are nonzero, and there is at least one
%  positive and one negative.
%
  if ( dist1 * dist2 < 0.0 & dist1 * dist3 < 0.0 )

    [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), ...
      t(1:dim_num,2), dist1, dist2, num_int, pint );

    [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,1), ...
      t(1:dim_num,3), dist1, dist3, num_int, pint );

  elseif ( dist2 * dist1 < 0.0 & dist2 * dist3 < 0.0 )

    [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,2), ...
      t(1:dim_num,1), dist2, dist1, num_int, pint );

    [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,2), ...
      t(1:dim_num,3), dist2, dist3, num_int, pint );

  elseif ( dist3 * dist1 < 0.0 & dist3 * dist2 < 0.0 )

    [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,3), ...
      t(1:dim_num,1), dist3, dist1, num_int, pint );

    [ num_int, pint ] = plane_imp_triangle_int_add_3d ( t(1:dim_num,3), ...
      t(1:dim_num,2), dist3, dist2, num_int, pint );

  end

  return
end


text:
function [cl2,classes] = image_histogram1d(varargin)
% Visualize a change_point map stored in hewma_cp.img
% (output from hewma2)
% and classify voxels into groupings based on CP
%
% :Usage:
% ::
%
%     [cl2,classes] = image_histogram1d([image name],[overlay])
%

name = 'hewma_cp.img';
%name = 'hewma_runlen.img';
overlay = [];

if length(varargin) > 0, name = varargin{1};,end
if length(varargin) > 1, overlay = varargin{2};,end

% ----------------------------------------------------
% initial viewing of CP map
% ----------------------------------------------------

[CLU.XYZ,CLU.XYZmm,CLU.Z,CLU.V] = img2voxel(name);
CLU.XYZ = CLU.XYZ(1:3,:);
%cl = mask2clusters(name);

%cluster_orthviews(cl,'overlay',overlay);
v2 = CLU.Z;

% ----------------------------------------------------
% load data
% ----------------------------------------------------

%v = spm_read_vols(spm_vol(name));
%v2 = v(:);

%discretize for convenience
v2 = round(v2.*100);
%v2(abs(v2) < eps | isnan(v2)) = [];  % shoul

% ----------------------------------------------------
% make histogram and get clusters (classifications) of CPs
% ----------------------------------------------------

nbins = input('Number of bins: '); %unique(v2);
tor_fig; f1 = gcf; [h,x] = hist(v2./100,nbins); hh = bar(x,h); set(hh,'FaceColor',[.5 .5 .5])

nclasses = input('How many classes?');
str = sprintf('%3.0f Classes: color image overlay starting at which class (e.g., 1 for all classes): ',nclasses)
startclass = input(str);

err = 1; indx = 1;
while err
    try
        classes = kmeans(v2, nclasses);   % ,'start','uniform');
        err = 0;
    catch
    end
    indx = indx + 1;
    if indx == 11, disp('kmeans: tried 10 times.  No solution.'); err = 0;, return, end
end

% For each x, find the modal class
% this is used to color the histogram
x100 = x*100;
binwid = diff(x100); binwid = binwid(1)./2;

for i =1:length(x100)
    xrange = [x100(i)-binwid x100(i)+binwid];
    xclass = classes(v2>xrange(1) & v2<=xrange(2));
    binclass(i) = mode(xclass);
end

%v2 = round(v2.*100);
% % classmap's elements are the range of values in v2
% mincp = min(v2);
% maxcp = max(v2);
% for i = mincp:maxcp, tmp = unique(classes(find(v2==i)));
%     if isempty(tmp), classmap(i) = 0;,
%     else,classmap(i) = tmp(1);, end
% end

% CLU = clusters2CLU(cl);
% CLU.Z = v2';
% CLU.Z(abs(CLU.Z)<eps) = NaN;

CLU.cp = v2';
v2(v2<=0) = 1;   %make sure no negative numbers or zeros for some reason
CLU.Z = classes';  %classmap(v2);

% ----------------------------------------------------
% define colors and sort by class size
% ----------------------------------------------------

%colors = {[1 0 0] [0 1 0] [0 0 1] [1 1 0] [1 0 1] [0 1 1]};
%while length(colors) < nclasses, colors = [colors colors];,end

colors = hot(nclasses+2 - startclass + 1);
colors = colors(2:end-1,:);

%un-plotted classes
colors = [repmat(.5,startclass-1,3); colors];

for i = 1:length(colors), col{i} = colors(i,:);,end
colors = col;
%colors = col(2:end);    %(floor(length(col)./3):end);

% sort by change point, ascending
for i = 1:nclasses,
    meancp(i) = mean(v2(classes==i));,
    nvox(i) = sum(classes==i);
    indx(i) = i;
end
[meancp,i] = sort(meancp,2,'ascend');
%nvox(i) = nvox(1:length(i));
%indx(i) = indx(1:length(i));
nvox  = nvox(i);
indx = indx(i);

meancp = meancp ./ 100; % convert back to original units.  resolution = 2 dec. places.

%table
fprintf(1,'\nClass\tNum. Voxels\tMean Value\tColor\t\t\n');
for j =1:nclasses
    fprintf(1,'%3.0f\t%3.0f\t%3.2f\t%3.2f\t%3.2f\t%3.2f\t\n',indx(j),nvox(j),meancp(j),colors{j}(1),colors{j}(2),colors{j}(3));
end
fprintf(1,'\n');

% sort by class size
%for i = 1:nclasses, nvox(i) = sum(classes==i);,end
%[nvox,i] = sort(nvox,2,'descend');

%colors(i) = colors(1:length(i));

% ----------------------------------------------------
% re-plot histogram with color codes
% ----------------------------------------------------
figure(f1);

clear i
for i = startclass:nclasses
    %     wh = find(classmap == indx(i)); range = [min(wh) max(wh)+1];
    %     wh = find(x*100 <= range(2) & x*100 >= range(1));
    %     hh = bar(x(wh),h(wh)); set(hh,'FaceColor',colors{i});

    wh = find(binclass == indx(i));
    hh = bar(x(wh),h(wh)); set(hh,'FaceColor',colors{i});

end
xlabel('Image value')
ylabel('Number of voxels')

% ----------------------------------------------------
% re-make separate clusters for each class
% and plot on brain
% ----------------------------------------------------

clear cl2
for i = startclass:nclasses
    CLUtmp = CLU;
    wh = find(CLUtmp.Z == indx(i));
    CLUtmp.XYZmm = CLUtmp.XYZmm(:,wh);
    CLUtmp.XYZ = CLUtmp.XYZ(:,wh);
    CLUtmp.Z = CLUtmp.Z(:,wh);
    CLUtmp.cp = CLUtmp.cp(:,wh);
    cl2{i} = tor_extract_rois([],CLUtmp,CLUtmp);

    if i == startclass
        cluster_orthviews(cl2{i},colors(i),'overlay',overlay);
    else
        cluster_orthviews(cl2{i},colors(i),'add','overlay',overlay);
    end
end

% write mask of sig. results
domask = input('Write mask of colormapped voxels? (1 or 0): ');
if domask
    xyz = [];
    for i=1:length(cl2)
        if ~isempty(cl2{i}), xyz = [xyz cl2{i}.XYZ];, end
    end
    V = spm_vol(name);
    V.fname = 'img_hist_mask.img';
    mask = voxel2mask(xyz', V.dim(1:3));
    spm_write_vol(V,mask);
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function g = gpsimLogLikeGradients(model)

% GPSIMLOGLIKEGRADIENTS Compute the gradients of the log likelihood of a GPSIM model.
% FORMAT
% DESC computes the gradients of the log likelihood of the given
% Gaussian process for use in a single input motif protein network.
% ARG model : the model for which the log likelihood is computed.
% RETURN g : the gradients of the parameters of the model.
%
% SEEALSO : gpsimCreate, gpsimLogLikelihood, gpsimGradient
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% MODIFIED : Pei Gao, 2008

% SHEFFIELDML

covGrad = -model.invK + model.invK*model.m*model.m'*model.invK;
covGrad = 0.5*covGrad;
if isfield(model, 'proteinPrior') && ~isempty(model.proteinPrior)
  g = kernGradient(model.kern, model.timesCell, covGrad);
else
  g = kernGradient(model.kern, model.t, covGrad);
end

% In case we need priors in.
% Add contribution of any priors
if isfield(model, 'bprior'),
  g = g + kernPriorGradient(model.kern);
end

gmuFull = model.m'*model.invK;

if isfield(model, 'proteinPrior') && ~isempty(model.proteinPrior)
  if model.includeNoise
    ind = model.kern.comp{1}.diagBlockDim{1} + (1:model.kern.comp{1}.diagBlockDim{2});
    gmu = zeros(size(1, model.numGenes));

    for i = 1:model.numGenes
      gmu(i) = sum(gmuFull(ind));
      ind = ind + model.kern.comp{1}.diagBlockDim{i+1};
    end
  else
    ind = model.kern.diagBlockDim{1} + (1:model.kern.diagBlockDim{2});
    gmu = zeros(size(1, model.numGenes));

    for i = 1:model.numGenes
      gmu(i) = sum(gmuFull(ind));
      ind = ind + model.kern.diagBlockDim{i+1};
    end
  end

else
  numData = size(model.t, 1);
  ind = 1:numData;
  gmu = zeros(size(1, model.numGenes));
  for i = 1:model.numGenes
    gmu(i) = sum(gmuFull(ind));
    ind = ind + numData;
  end
end

gb = gmu./model.D;
fhandle = str2func([model.bTransform 'Transform']);
% In case we need priors in.
% Add prior on B if it exists.
if isfield(model, 'bprior');
  gb = gb + priorGradient(model.bprior, model.B);
end

gb = gb.*fhandle(model.B, 'gradfact');

% This is a nasty hack to add the influence of the D in the mean to
% the gradient already computed for the kernel. This is all very
% clunky and sensitive to changes that take place elsewhere in the
% code ...
gd = -gmu.*model.B./(model.D.*model.D);
if model.kern.numBlocks == 1
  decayIndices = 1;
elseif model.kern.numBlocks>1
  if isfield(model, 'proteinPrior') && ~isempty(model.proteinPrior)
    decayIndices = [3];
    for i = 3:model.kern.numBlocks
      decayIndices(end+1) = decayIndices(end) + 2;
    end
  else
    decayIndices = [1 4];
    for i = 3:model.kern.numBlocks
      decayIndices(end+1) = decayIndices(end) + 2;
    end
  end
end

% Account for decay in mean.
g(decayIndices) = g(decayIndices) ...
    + gd.*expTransform(model.D, 'gradfact');

g = [g gb];

if isfield(model, 'fix')
  for i = 1:length(model.fix)
    g(model.fix(i).index) = 0;
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Comparing Gaussian Blur and Guided Image Filtering
% See http://dsp.stackexchange.com/questions/29041

%% General Parameters and Initialization

clear();
close('all');

% set(0, 'DefaultFigureWindowStyle', 'docked');
defaultLooseInset = get(0, 'DefaultAxesLooseInset');
% set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

titleFontSize   = 14;
axisFotnSize    = 12;
stringFontSize  = 12;

thinLineWidth   = 2;
normalLineWidth = 3;
thickLineWidth  = 4;

smallSizeData   = 36;
mediumSizeData  = 48;
bigSizeData     = 60;

randomNumberStream  = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber     = round(sum(clock()));
% subStreamNumber    = 57162;
% subStreamNumber    = 2143;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

%% Setting Constants

FALSE   = 0;
TRUE    = 1;

OFF = 0;
ON  = 1;

RADIUS_TO_STD_FACTOR    = 5;

%% Setting Parameters

numRows = 400;
numCols = 400;

grayLevel1 = 245 / 255;
grayLevel2 = 10 / 255;

noiseStd = 10 / 255;

refLine1Row     = 201;
vRefLine1ColIdx = [1:numCols];

gaussianFilterRadius    = 7;
vGuidedFilterRadius     = [7, 7];
guidedFilterSmoothing   = 0.01;

%% Creating Data

borderColIdx = round(numCols / 2);

% Refrence Image
mRefImage = zeros([numRows, numCols]);
mRefImage(:, 1:borderColIdx)            = grayLevel1;
mRefImage(:, (borderColIdx + 1):end)    = grayLevel2;

% Noisy Image
mNoisyImage = mRefImage + (noiseStd * randn([numRows, numCols]));

% Gaussian Filtered Image
mGaussFilteredImage = ApplyGaussianBlur(mNoisyImage, gaussianFilterRadius, RADIUS_TO_STD_FACTOR);

% Guided Filtered Image
mGuidedFilteredImage = imguidedfilter(mNoisyImage, 'NeighborhoodSize', vGuidedFilterRadius, 'DegreeOfSmoothing', guidedFilterSmoothing);

%% Displaying Results

% Displaying Reference Image
mImgDisplay                                     = repmat(mRefImage, [1, 1, 3]);
mImgDisplay(refLine1Row, vRefLine1ColIdx, 1)    = 1;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 2)    = 0;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 3)    = 0;
for iRow = 1:numRows
    mImgDisplay(iRow, iRow, 1) = 0;
    mImgDisplay(iRow, iRow, 2) = 1;
    mImgDisplay(iRow, iRow, 3) = 0;
end

hFigure         = figure();
set(hFigure, 'Units', 'pixels', 'Position', [100, 100, 500, 500]);
hAxes           = axes('Units', 'pixels', 'Position', [50, 50, numCols, numRows]);
hImageObject    = image(mImgDisplay);
set(get(hAxes, 'Title'), 'String', 'Reference Image', 'FontSize', titleFontSize);

% Extract Line 01
vLine1 = mRefImage(refLine1Row, vRefLine1ColIdx);
% Extract Line 02
vLine2 = zeros([numRows, 1]);
for iRow = 1:numRows
    vLine2(iRow) = mRefImage(iRow, iRow);
end

hFigure = figure();
hAxes   = axes();
hLineSeries = plot([1:numCols], vLine1, [1:numRows], vLine2);
set(hLineSeries(1), 'LineWidth', normalLineWidth, 'Color', 'r');
set(hLineSeries(2), 'LineWidth', normalLineWidth, 'Color', 'g');
set(get(hAxes, 'Title'), 'String', 'Values Across Lines - Reference Image', 'FontSize', titleFontSize);
hLegend = legend({['Line #01'], ['Line #02']});
set(hLegend, 'FontSize', axisFotnSize);

% Displaying Noisy Image
mImgDisplay                                     = repmat(mNoisyImage, [1, 1, 3]);
mImgDisplay(refLine1Row, vRefLine1ColIdx, 1)    = 1;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 2)    = 0;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 3)    = 0;
for iRow = 1:numRows
    mImgDisplay(iRow, iRow, 1) = 0;
    mImgDisplay(iRow, iRow, 2) = 1;
    mImgDisplay(iRow, iRow, 3) = 0;
end

hFigure         = figure();
set(hFigure, 'Units', 'pixels', 'Position', [100, 100, 500, 500]);
hAxes           = axes('Units', 'pixels', 'Position', [50, 50, numCols, numRows]);
hImageObject    = image(mImgDisplay);
set(get(hAxes, 'Title'), 'String', 'Noisy Image', 'FontSize', titleFontSize);

% Extract Line 01
vLine1 = mNoisyImage(refLine1Row, vRefLine1ColIdx);
% Extract Line 02
vLine2 = zeros([numRows, 1]);
for iRow = 1:numRows
    vLine2(iRow) = mNoisyImage(iRow, iRow);
end

hFigure = figure();
hAxes   = axes();
hLineSeries = plot([1:numCols], vLine1, [1:numRows], vLine2);
set(hLineSeries(1), 'LineWidth', normalLineWidth, 'Color', 'r');
set(hLineSeries(2), 'LineWidth', normalLineWidth, 'Color', 'g');
set(get(hAxes, 'Title'), 'String', 'Values Across Lines - Noisy Image', 'FontSize', titleFontSize);
hLegend = legend({['Line #01'], ['Line #02']});
set(hLegend, 'FontSize', axisFotnSize);

% Displaying Gaussian Filtered Image
mImgDisplay                                     = repmat(mGaussFilteredImage, [1, 1, 3]);
mImgDisplay(refLine1Row, vRefLine1ColIdx, 1)    = 1;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 2)    = 0;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 3)    = 0;
for iRow = 1:numRows
    mImgDisplay(iRow, iRow, 1) = 0;
    mImgDisplay(iRow, iRow, 2) = 1;
    mImgDisplay(iRow, iRow, 3) = 0;
end

hFigure         = figure();
set(hFigure, 'Units', 'pixels', 'Position', [100, 100, 500, 500]);
hAxes           = axes('Units', 'pixels', 'Position', [50, 50, numCols, numRows]);
hImageObject    = image(mImgDisplay);
set(get(hAxes, 'Title'), 'String', 'Gaussian Filtered Image', 'FontSize', titleFontSize);

% Extract Line 01
vLine1 = mGaussFilteredImage(refLine1Row, vRefLine1ColIdx);
% Extract Line 02
vLine2 = zeros([numRows, 1]);
for iRow = 1:numRows
    vLine2(iRow) = mGaussFilteredImage(iRow, iRow);
end

hFigure = figure();
hAxes   = axes();
hLineSeries = plot([1:numCols], vLine1, [1:numRows], vLine2);
set(hLineSeries(1), 'LineWidth', normalLineWidth, 'Color', 'r');
set(hLineSeries(2), 'LineWidth', normalLineWidth, 'Color', 'g');
set(get(hAxes, 'Title'), 'String', 'Values Across Lines - Gaussian Filtered Image', 'FontSize', titleFontSize);
hLegend = legend({['Line #01'], ['Line #02']});
set(hLegend, 'FontSize', axisFotnSize);

% Displaying Guided Filtered Image
mImgDisplay                                     = repmat(mGuidedFilteredImage, [1, 1, 3]);
mImgDisplay(refLine1Row, vRefLine1ColIdx, 1)    = 1;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 2)    = 0;
mImgDisplay(refLine1Row, vRefLine1ColIdx, 3)    = 0;
for iRow = 1:numRows
    mImgDisplay(iRow, iRow, 1) = 0;
    mImgDisplay(iRow, iRow, 2) = 1;
    mImgDisplay(iRow, iRow, 3) = 0;
end

hFigure         = figure();
set(hFigure, 'Units', 'pixels', 'Position', [100, 100, 500, 500]);
hAxes           = axes('Units', 'pixels', 'Position', [50, 50, numCols, numRows]);
hImageObject    = image(mImgDisplay);
set(get(hAxes, 'Title'), 'String', 'Guided Filtered Image', 'FontSize', titleFontSize);

% Extract Line 01
vLine1 = mGuidedFilteredImage(refLine1Row, vRefLine1ColIdx);
% Extract Line 02
vLine2 = zeros([numRows, 1]);
for iRow = 1:numRows
    vLine2(iRow) = mGuidedFilteredImage(iRow, iRow);
end

hFigure = figure();
hAxes   = axes();
hLineSeries = plot([1:numCols], vLine1, [1:numRows], vLine2);
set(hLineSeries(1), 'LineWidth', normalLineWidth, 'Color', 'r');
set(hLineSeries(2), 'LineWidth', normalLineWidth, 'Color', 'g');
set(get(hAxes, 'Title'), 'String', 'Values Across Lines - Guided Filtered Image', 'FontSize', titleFontSize);
hLegend = legend({['Line #01'], ['Line #02']});
set(hLegend, 'FontSize', axisFotnSize);

%% Restore Defaults
set(0, 'DefaultFigureWindowStyle', 'normal');
set(0, 'DefaultAxesLooseInset', defaultLooseInset);

author: RoyiAvital

repo: StackExchangeCodes